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.
if(.Platform$OS.type != "windows"){
Rsubread::buildindex(basename=file.path(data_dir, "ERCC_index"), reference=ERCCfa_fn)
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.12.3
##
## //================================= setting ==================================\\
## || ||
## || Index name : ERCC_index ||
## || Index space : base space ||
## || Index split : no-split ||
## || Repeat threshold : 100 repeats ||
## || Gapped index : no ||
## || ||
## || Free / total memory : 69.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/RtmpuUvI9J/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=134.1k bps/s ||
## || 49%, 0 mins elapsed, rate=774.2k bps/s ||
## || 66%, 0 mins elapsed, rate=1015.2k bps/s ||
## || 3.0 GB of memory is needed for index building. ||
## || Build the index... ||
## || 8%, 0 mins elapsed, rate=13.4k bps/s ||
## || 49%, 0 mins elapsed, rate=79.9k bps/s ||
## || 66%, 0 mins elapsed, rate=106.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/RtmpuUvI9J/ERCC_index was successfully built. ||
## || ||
## \\============================================================================//
##
##
## ========== _____ _ _ ____ _____ ______ _____
## ===== / ____| | | | _ \| __ \| ____| /\ | __ \
## ===== | (___ | | | | |_) | |__) | |__ / \ | | | |
## ==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
## ==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
## ========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
## Rsubread 2.12.3
##
## //================================= 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 (07-Mar-2023 18:02:10, pid=2647699) ================\\
## || ||
## || 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.5 mins elapsed, rate=131.9k reads per second ||
## || 43% completed, 0.5 mins elapsed, rate=156.6k reads per second ||
## || 62% completed, 0.5 mins elapsed, rate=170.6k reads per second ||
## || 82% completed, 0.5 mins elapsed, rate=176.4k reads per second ||
## || 72% completed, 0.5 mins elapsed, rate=0.8k reads per second ||
## || 78% completed, 0.5 mins elapsed, rate=0.9k reads per second ||
## || 83% completed, 0.5 mins elapsed, rate=0.9k reads per second ||
## || 88% completed, 0.5 mins elapsed, rate=1.0k reads per second ||
## || 92% completed, 0.5 mins elapsed, rate=1.0k reads per second ||
## || 97% completed, 0.5 mins elapsed, rate=1.1k reads per second ||
## || 102% completed, 0.5 mins elapsed, rate=1.1k reads per second ||
## || 107% completed, 0.5 mins elapsed, rate=1.2k reads per second ||
## || 113% completed, 0.5 mins elapsed, rate=1.2k 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.5 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.
if(.Platform$OS.type != "windows"){
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.
if(.Platform$OS.type != "windows"){
sc_demultiplex(file.path(data_dir, "out.map.bam"), data_dir, barcode_annotation_fn, has_UMI=FALSE)
sc_gene_counting(data_dir, barcode_annotation_fn)
}
## demultiplexing reads by barcode...
## Warning: mitochondrial chromosome not found using chromosome name `MT`.
## time elapsed: 78 milliseconds
##
## summarising gene counts...
## time elapsed: 57 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).
