Aligning reads to a reference genome
Next we align reads to the genome. This example uses
Rsubread but any aligner that support RNA-seq alignment and
gives standard BAM output can be used here. The Rsubread
package is not available on Windows, so the following alignment step
will only run on Mac OSX or Linux.
Rsubread::buildindex(basename=file.path(data_dir, "ERCC_index"), reference=ERCCfa_fn)
##
## ========== _____ _ _ ____ _____ ______ _____
## ===== / ____| | | | _ \| __ \| ____| /\ | __ \
## ===== | (___ | | | | |_) | |__) | |__ / \ | | | |
## ==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
## ==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
## ========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
## Rsubread 2.21.0
##
## //================================= setting ==================================\\
## || ||
## || Index name : ERCC_index ||
## || Index space : base space ||
## || Index split : no-split ||
## || Repeat threshold : 100 repeats ||
## || Gapped index : no ||
## || ||
## || Free / total memory : 48.0GB / 125.4GB ||
## || ||
## || Input files : 1 file in total ||
## || o ERCC92.fa ||
## || ||
## \\============================================================================//
##
## //================================= Running ==================================\\
## || ||
## || Check the integrity of provided reference sequences ... ||
## || There were 1 notes for reference sequences. ||
## || The notes can be found in the log file, '/tmp/RtmpF20gvE/ERCC_index.log'. ||
## || Scan uninformative subreads in reference sequences ... ||
## || 1 uninformative subreads were found. ||
## || These subreads were excluded from index building. ||
## || Estimate the index size... ||
## || 8%, 0 mins elapsed, rate=121.2k bps/s ||
## || 49%, 0 mins elapsed, rate=702.7k bps/s ||
## || 66%, 0 mins elapsed, rate=924.0k bps/s ||
## || 3.0 GB of memory is needed for index building. ||
## || Build the index... ||
## || 8%, 0 mins elapsed, rate=12.5k bps/s ||
## || 49%, 0 mins elapsed, rate=74.6k bps/s ||
## || 66%, 0 mins elapsed, rate=99.1k bps/s ||
## || Save current index block... ||
## || [ 0.0% finished ] ||
## || [ 10.0% finished ] ||
## || [ 20.0% finished ] ||
## || [ 30.0% finished ] ||
## || [ 40.0% finished ] ||
## || [ 50.0% finished ] ||
## || [ 60.0% finished ] ||
## || [ 70.0% finished ] ||
## || [ 80.0% finished ] ||
## || [ 90.0% finished ] ||
## || [ 100.0% finished ] ||
## || ||
## || Total running time: 0.2 minutes. ||
## || Index /tmp/RtmpF20gvE/ERCC_index was successfully built. ||
## || ||
## \\============================================================================//
Rsubread::align(index=file.path(data_dir, "ERCC_index"),
readfile1=file.path(data_dir, "combined.fastq.gz"),
output_file=file.path(data_dir, "out.aln.bam"), phredOffset=64)
##
## ========== _____ _ _ ____ _____ ______ _____
## ===== / ____| | | | _ \| __ \| ____| /\ | __ \
## ===== | (___ | | | | |_) | |__) | |__ / \ | | | |
## ==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
## ==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
## ========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
## Rsubread 2.21.0
##
## //================================= setting ==================================\\
## || ||
## || Function : Read alignment (RNA-Seq) ||
## || Input file : combined.fastq.gz ||
## || Output file : out.aln.bam (BAM) ||
## || Index name : ERCC_index ||
## || ||
## || ------------------------------------ ||
## || ||
## || Threads : 1 ||
## || Phred offset : 64 ||
## || Min votes : 3 / 10 ||
## || Max mismatches : 3 ||
## || Max indel length : 5 ||
## || Report multi-mapping reads : yes ||
## || Max alignments per multi-mapping read : 1 ||
## || ||
## \\============================================================================//
##
## //=============== Running (29-Oct-2024 19:28:13, pid=1772257) ================\\
## || ||
## || Check the input reads. ||
## || The input file contains base space reads. ||
## || Initialise the memory objects. ||
## || Estimate the mean read length. ||
## || The range of Phred scores observed in the data is [1,1] ||
## || Create the output BAM file. ||
## || Check the index. ||
## || Init the voting space. ||
## || Global environment is initialised. ||
## || Load the 1-th index block... ||
## || The index block has been loaded. ||
## || Start read mapping in chunk. ||
## || 23% completed, 0.6 mins elapsed, rate=137.1k reads per second ||
## || 43% completed, 0.6 mins elapsed, rate=169.5k reads per second ||
## || 62% completed, 0.6 mins elapsed, rate=185.8k reads per second ||
## || 82% completed, 0.6 mins elapsed, rate=194.1k reads per second ||
## || 72% completed, 0.6 mins elapsed, rate=0.6k reads per second ||
## || 78% completed, 0.6 mins elapsed, rate=0.7k reads per second ||
## || 83% completed, 0.6 mins elapsed, rate=0.7k reads per second ||
## || 88% completed, 0.6 mins elapsed, rate=0.8k reads per second ||
## || 92% completed, 0.6 mins elapsed, rate=0.8k reads per second ||
## || 97% completed, 0.6 mins elapsed, rate=0.9k reads per second ||
## || 102% completed, 0.6 mins elapsed, rate=0.9k reads per second ||
## || 107% completed, 0.6 mins elapsed, rate=1.0k reads per second ||
## || 113% completed, 0.6 mins elapsed, rate=1.0k reads per second ||
## || ||
## || Completed successfully. ||
## || ||
## \\==================================== ====================================//
##
## //================================ Summary =================================\\
## || ||
## || Total reads : 46,009 ||
## || Mapped : 46,009 (100.0%) ||
## || Uniquely mapped : 46,009 ||
## || Multi-mapping : 0 ||
## || ||
## || Unmapped : 0 ||
## || ||
## || Indels : 0 ||
## || ||
## || Running time : 0.6 minutes ||
## || ||
## \\============================================================================//
## out.aln.bam
## Total_reads 46009
## Mapped_reads 46009
## Uniquely_mapped_reads 46009
## Multi_mapping_reads 0
## Unmapped_reads 0
## Indels 0
Assigning reads to annotated exons
After the read alignment, we assign reads to exons based on the
annotation provided using the sc_exon_mapping function.
Currently scPipe only supports annotation in
gff3 or bed format.
sc_exon_mapping(file.path(data_dir, "out.aln.bam"),
file.path(data_dir, "out.map.bam"),
ERCCanno_fn)
## adding annotation files...
## time elapsed: 0 milliseconds
##
## annotating exon features...
## updating progress every 3 minutes...
## 46009 reads processed, 46k reads/sec
## number of read processed: 46009
## unique map to exon: 46009 (100.00%)
## ambiguous map to multiple exon: 0 (0.00%)
## map to intron: 0 (0.00%)
## not mapped: 0 (0.00%)
## unaligned: 0 (0.00%)
De-multiplexing data and counting genes
Next we use the sc_demultiplex function to split the
reads into separate .csv files for each cell in the
/count subfolder. Each file contains three columns and each
row corresponds to a distinct read. The first column contains the gene
ID that the read maps to, the second column gives the UMI sequence for
that read and the third column is the distance to the transcript end
position (TES) in bp. This file will be used for UMI deduplication and
to generate a gene count matrix by calling the
sc_gene_counting function.
sc_demultiplex(file.path(data_dir, "out.map.bam"), data_dir, barcode_annotation_fn, has_UMI=FALSE)
## demultiplexing reads by barcode...
## Warning: mitochondrial chromosome not found using chromosome name `MT`.
## time elapsed: 90 milliseconds
sc_gene_counting(data_dir, barcode_annotation_fn)
## summarising gene counts...
## time elapsed: 61 milliseconds
We have now completed the preprocessing steps. The gene count matrix
is available as a .csv file in
data_dir/gene_count.csv and quality control statistics are
saved in the data_dir/stat folder. These data are useful
for later quality control (QC).
