Single-cell_read_alignment_and_DropEst_library_quantification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## __Section 1:__ Single-cell_read_alignment_and_DropEst_library_quantification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this section is to set up the containerized environment for running the DropEst pipeline, running the the process of barcode demultiplexing, read alignment, and droplet count estimation. The inputs to this step are the raw sequencing files and the outputs to this step are an estimated droplet read count matrix and its corresponding quality control reports."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Before You Begin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This protocol requires three types of input files, which will have to be in an accessible directory:\n",
    "    \n",
    "A primary assembly .fa file corresponding to the reference genome for read alignment. Here we use release 85 of the Ensembl human reference genome.\n",
    "\n",
    "    - i.e.: Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa\n",
    "    - Located at ftp://ftp.ensembl.org/pub/release-85/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa.gz \n",
    "    \n",
    "A gene annotation file corresponding to the same reference genome and release version\n",
    "\n",
    "    - i.e.: Homo_sapiens.GRCh38.85.annotated.gtf\n",
    "    - Located at https://www.mc.vanderbilt.edu/vumcdept/cellbio/laulab/data/Homo_sapiens.GRCh38.85.annotated.gtf.gz\n",
    "    \n",
    "Two .fastq files, corresponding to the paired R1 and R2 sequencing reads of the target sequenced single-cell RNA-seq library. \n",
    "\n",
    "    - i.e.: 3907-AS-1-CCGCGGTT-AGCGCTAG_S343_sub_R1_001.fastq and 3907-AS-1-CCGCGGTT-AGCGCTAG_S343_sub_R2_001.fastq\n",
    "    - Located at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM5068493\n",
    "\n",
    "__Note that the above is a subsampled sequencing library for the sake of testing Step 1 of this pipeline.__\n",
    "\n",
    "For the full version of the 3907 dataset, which may take several hours to run depending on the available hardware, see the following information:\n",
    "\n",
    "    - i.e.: 3907-AS-1-CCGCGGTT-AGCGCTAG_S343_R1_001.fastq and 3907-AS-1-CCGCGGTT-AGCGCTAG_S343_R2_001.fastq\n",
    "    - Located at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM4804820\n",
    "    \n",
    "Further exploration of dropEst configurations can be done after cloning the dropEst repository to a local directory with the following command:\n",
    "\n",
    "    - git clone https://github.com/kharchenkolab/dropEst.git"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Singularity Installation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use the provided containerized environment, we need to install the required golang libraries and Singularity itself. For Ubuntu, a Debian-based linux distribution, we first update our software libraries through apt-get"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#If you do not have administrative privleges, please coordinate with your system administrator to update your Ubuntu packages\n",
    "!sudo apt-get update && sudo apt-get install -y \\build-essential \\libssl-dev \\uuid-dev \\libgpgme11-dev \\squashfs-tools \\libseccomp-dev \\pkg-config"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next install golang, by obtaining the installation files, extracting them, and updating paths for usage. In this case version 1.15.6 is used. Further instructions can be found at: https://golang.org/doc/install"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Download to current directory\n",
    "!wget https://dl.google.com/go/go1.15.6.linux-amd64.tar.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Extract files to local directory\n",
    "!sudo tar -C /usr/local -xzf go1.15.6.linux-amd64.tar.gz\n",
    "#Alternatively extract to a user directory\n",
    "!tar -C /home/$USER -xzf go1.15.6.linux-amd64.tar.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Update path\n",
    "!export PATH=$PATH:/usr/local/go/bin\n",
    "#If saved to a user directory\n",
    "!export PATH=$PATH:/home/$USER/go/bin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Double check installation by checking go version\n",
    "!go version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can install Singularity from a release, in this case version 3.7.0 is used. Further instructions can be found at: https://sylabs.io/guides/3.0/user-guide/installation.html#download-and-install-singularity-from-a-release"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Download and extract the necessary Singularity libraries to a user directory, bypassing requirement for admin access\n",
    "!export VERSION=3.7.0 && # adjust this as necessary \\\n",
    "    mkdir -p /home/$USER/go/src/github.com/sylabs && \\\n",
    "    cd /home/$USER/go/src/github.com/sylabs && \\\n",
    "    wget https://github.com/sylabs/singularity/releases/download/v${VERSION}/singularity-${VERSION}.tar.gz && \\\n",
    "    tar -xzf singularity-${VERSION}.tar.gz && \\\n",
    "    cd ./singularity && \\\n",
    "    ./mconfig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Navigate to the installation directory and compile\n",
    "!cd /home/$USER/go/src/github.com/sylabs/singularity/builddir && \\\n",
    "    make && \\\n",
    "    make install"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, pull the Singularity image from dockerhub, this may take a while depending on download speeds. This container comes preconfigured with DropEst, STAR, and scRNABatchQC, along with their supporting libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!singularity pull docker://ramiremars/star_dropest:star_protocols_pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ensure that your working directory contains the following, if following the example detailed in our manuscript in STAR Protocols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "3907-AS-1-CCGCGGTT-AGCGCTAG_S343_sub_R1_001.fastq #Read 1 \n",
    "3907-AS-1-CCGCGGTT-AGCGCTAG_S343_sub_R2_001.fastq #Read 2\n",
    "Homo_sapiens.GRCh38.85.annotated.gtf #GTF file\n",
    "Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa #Reference genome\n",
    "star_dropest_star_protocols_pipeline.sif #Singularity image file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### __1.__ Running DropTag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DropTag is first run to demultiplex the .fastq read files. This is the first step of the DropEst pipeline, and further documentation can be found at: https://github.com/hms-dbmi/dropEst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!singularity exec -e star_dropest_star_protocols_pipeline.sif droptag -c /usr/share/dropEst/configs/indrop_v1_2.xml \\\n",
    "    3907-AS-1-CCGCGGTT-AGCGCTAG_S343_sub_R1_001.fastq 3907-AS-1-CCGCGGTT-AGCGCTAG_S343_sub_R2_001.fastq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The primary output from this step is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "3907-AS-1-CCGCGGTT-AGCGCTAG_S343_sub_R2_001.fastq.tagged.1.fastq.gz #Tagged fastq for demultiplexing and droplet estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### __2.__ Generating the STAR index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The STAR aligner used for read alignment requires an index file to be generated, we generate it with the following command, which creates a directory to operate within. Further documentation can be found at: https://github.com/alexdobin/STAR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir STAR_index && singularity exec -e star_dropest_star_protocols_pipeline.sif STAR --runThreadN 16 --runMode genomeGenerate \\\n",
    "    --genomeDir STAR_index --genomeFastaFiles Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa --sjdbGTFfile Homo_sapiens.GRCh38.85.annotated.gtf --sjdbOverhang 99"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main output from this step is A directory, named __STAR_index__, which contains several files related to the generated genomic index."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### __3.__ Read alignment with the STAR aligner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "STAR aligner is a splice-aware aligner which we use to align sequenced reads to the reference genome. In this example, there is only one tagged fastq file, this number will increase with read depth, the --readFilesIn parameter will need to be adjusted accordingly to include all such tagged fastq files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!singularity exec -e star_dropest_star_protocols_pipeline.sif STAR --genomeDir STAR_index \\\n",
    "    --readFilesIn 3907-AS-1-CCGCGGTT-AGCGCTAG_S343_sub_R2_001.fastq.tagged.1.fastq.gz  --outSAMmultNmax 1 \\\n",
    "    --runThreadN 16 --readNameSeparator space --outSAMunmapped Within --outSAMtype BAM SortedByCoordinate \\\n",
    "    --outFileNamePrefix 3907-AS-1-CCGCGGTT-AGCGCTAG_S343_sub  --readFilesCommand gunzip -c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main output from this step is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "3907-AS-1-CCGCGGTT-AGCGCTAG_S343_subAligned.sortedByCoord.out.bam #A sorted bam file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### __4.__ Droplet count estimation with DropEst"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "STAR aligner is a splice-aware aligner which we use to align sequenced reads to the reference genome. This step may take several hours depending on the read depth of the library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!singularity exec -e star_dropest_star_protocols_pipeline.sif dropest -m -V -b -F -o sample_name \\\n",
    "    -g Homo_sapiens.GRCh38.85.annotated.gtf -L eiEIBA -c /usr/share/dropEst/configs/indrop_v1_2.xml 3907-AS-1-CCGCGGTT-AGCGCTAG_S343_subAligned.sortedByCoord.out.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main outputs from this step are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "3907-AS-1-CCGCGGTT-AGCGCTAG_S343_subAligned.sortedByCoord.out.tagged.bam #A sorted and tagged .bam file\n",
    "3907-AS-1-CCGCGGTT-AGCGCTAG_S343_subAligned.sortedByCoord.out.filtered.bam #A sorted, tagged, and filtered .bam file\n",
    "sample_name.rds #A .rds file which contains aggregated count matrix information"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### __5.__ RDS conversion to sparse droplet count matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of DropEst used in this step contains the droplet matrix that will go into downstream scRNA-seq analyses and further droplet matrix quality control. Our script, RDS_to_sparsematrix.r generates a .mtx along with its corresponding column and row names so that analysis can be performed in non-R environments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!singularity exec -e star_dropest_star_protocols_pipeline.sif R --vanilla --slave -f /R_scripts/RDS_to_sparesematrix.r --args sample_name.rds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main outputs from this step are two directories named __cmraw_sparse__ and __cm_sparse__ which contain sparse matrices, the former being completely unfiltered and the latter with a baseline level of filtering. Further a .csv based on the cm_sparse is also generated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_name.rds_cm.csv #A dense matrix of the droplet counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### __6.__ scRNABatchQC report generation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, initial droplet estimation QC metrics can be generated from the output R file from step 1.5. Further documentation can be found at: https://github.com/liuqivandy/scRNABatchQC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!singularity exec -e star_dropest_star_protocols_pipeline.sif R --vanilla --slave -f /R_scripts/scRNABatchQC.r --args hsapiens sample_name.rds_cm.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main output from this step is a directory named __scRNABatchQC_report_files__, which contains a report regarding the quality of the input dense matrix "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}