scRNA-seq: RNA quantification

1 RNA quantification

Like for bulk RNA sequencing, we start the analysis of scRNASeq data with a list of fastq files.

For the first part of this practical, we are going to work on a 1k cells 1:1 Mixture of Fresh Frozen Human (HEK293T) and Mouse (NIH3T3) Cells available from the 10x website. First, we download the fastq files (6.34 GB).

mkdir data/hgmm_1k
wget http://cf.10xgenomics.com/samples/cell-exp/2.1.0/hgmm_1k/hgmm_1k_fastqs.tar -O data/hgmm_1k/hgmm_1k_fastqs.tar
cd data/hgmm_1k/
tar -xvf hgmm_1k_fastqs.tar
cd ../..

On pedago-ngs.univ-lyon1.fr, you can find the file here

mkdir data
ln -s /data/share/MADT/TP_lmodolo/data/hgmm_1k/ data/hgmm_1k

This is a toy dataset, which is not representative of the size of the data generated in scRNASeq.

Check the size of the data generated for the Single-cell RNA-seq of 1.3 million brain cells from E18 mice.

The processing of scRNASeq fastq file can be a tedious and time-consuming process. In addition to the size of the data to handle, for each read, we have to:

• locate the read on the reference genome or transcritome
• identify the cell molecular barcode
• identify the UMI
• identify the sample molecular barcode (if any)

To ease these steps, the barcode are often at the start or the end of the reads, where the read quality is the lowest. Most scRNASeq pipeline will clip the barcode / UMI information from the read and add them to the read name. Therefore, the read can be mapped and we can assign it to a given cell and remove PCR duplicate afterward.

Luckily for you, the kb pipeline (Kallisto-Bustools) can handle most of these steps in few commands (Melsted et al.).

2 Reference indexing

The kb pipeline uses Kallisto. Therefore, we need a reference index before processing the data. The dataset is a mixture of Human and Mouse cells; thus, we want to quantify the transcript of a Human-Mouse hybrid reference. We can run the following command to build the index with the kallistobustools container (docker://lbmc/tp_scrna:0.1.0)

wget ftp://ftp.ensembl.org/pub/release-99/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz -O data/hgmm_1k/hs_cdna.fa.gz
wget ftp://ftp.ensembl.org/pub/release-99/fasta/mus_musculus/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz -O data/hgmm_1k/mm_cdna.fa.gz

kallisto index -i data/hgmm_1k/hs_mm_tr_index.idx data/hgmm_1k/hs_cdna.fa.gz data/hgmm_1k/mm_cdna.fa.gz

[build] loading fasta file data/hgmm_1k/hs_cdna.fa.gz
[build] loading fasta file data/hgmm_1k/mm_cdna.fa.gz
[build] k-mer length: 31
[build] warning: clipped off poly-A tail (longer than 10)
        from 2130 target sequences
[build] warning: replaced 8 non-ACGUT characters in the input sequence
        with pseudorandom nucleotides
[build] counting k-mers ... done.
[build] building target de Bruijn graph ...  done
[build] creating equivalence classes ...  done
[build] target de Bruijn graph has 2192295 contigs and contains 209728067 k-mers

3 RNA quantification

Here we will generate the bus file. bus stands for Barbode, UMI, Set. In text form, it is a table whose first column is the barcode. The second column is the UMI that are associated with the barcode. The third column is the index of the equivalence class reads with the UMI maps. The fourth column is the count of reads with this barcode, UMI, and equivalence class combination, which is ignored as one UMI should stand for one molecule. See this paper for more detail.

Also, since the reads were generated with the 10x Genomics Chromium Single Cell v2 Chemistry, the -x 10xv2 argument is used. To view others supported technologies, and see the information kb needs, run kallisto bus --list.

mkdir -p results/hgmm_1k/
kallisto bus \
  -i data/hgmm_1k/hs_mm_tr_index.idx \
  -o results/hgmm_1k/ -x 10xv2 -t 10 \
  data/hgmm_1k/fastqs/hgmm_1k_S1_L001_R1_001.fastq.gz data/hgmm_1k/fastqs/hgmm_1k_S1_L001_R2_001.fastq.gz \
  data/hgmm_1k/fastqs/hgmm_1k_S1_L002_R1_001.fastq.gz data/hgmm_1k/fastqs/hgmm_1k_S1_L002_R2_001.fastq.gz \
  data/hgmm_1k/fastqs/hgmm_1k_S1_L003_R1_001.fastq.gz data/hgmm_1k/fastqs/hgmm_1k_S1_L003_R2_001.fastq.gz \
  data/hgmm_1k/fastqs/hgmm_1k_S1_L004_R1_001.fastq.gz data/hgmm_1k/fastqs/hgmm_1k_S1_L004_R2_001.fastq.gz \
  data/hgmm_1k/fastqs/hgmm_1k_S1_L005_R1_001.fastq.gz data/hgmm_1k/fastqs/hgmm_1k_S1_L005_R2_001.fastq.gz \
  data/hgmm_1k/fastqs/hgmm_1k_S1_L006_R1_001.fastq.gz data/hgmm_1k/fastqs/hgmm_1k_S1_L006_R2_001.fastq.gz \
  data/hgmm_1k/fastqs/hgmm_1k_S1_L007_R1_001.fastq.gz data/hgmm_1k/fastqs/hgmm_1k_S1_L007_R2_001.fastq.gz \
  data/hgmm_1k/fastqs/hgmm_1k_S1_L008_R1_001.fastq.gz data/hgmm_1k/fastqs/hgmm_1k_S1_L008_R2_001.fastq.gz

[index] k-mer length: 31
[index] number of targets: 309,785
[index] number of k-mers: 209,728,067
[index] number of equivalence classes: 1,288,400
[quant] will process sample 1: data/hgmm_1k/fastqs/hgmm_1k_S1_L001_R1_001.fastq.gz
                               data/hgmm_1k/fastqs/hgmm_1k_S1_L001_R2_001.fastq.gz
[quant] will process sample 2: data/hgmm_1k/fastqs/hgmm_1k_S1_L002_R1_001.fastq.gz
                               data/hgmm_1k/fastqs/hgmm_1k_S1_L002_R2_001.fastq.gz
[quant] will process sample 3: data/hgmm_1k/fastqs/hgmm_1k_S1_L003_R1_001.fastq.gz
                               data/hgmm_1k/fastqs/hgmm_1k_S1_L003_R2_001.fastq.gz
[quant] will process sample 4: data/hgmm_1k/fastqs/hgmm_1k_S1_L004_R1_001.fastq.gz
                               data/hgmm_1k/fastqs/hgmm_1k_S1_L004_R2_001.fastq.gz
[quant] will process sample 5: data/hgmm_1k/fastqs/hgmm_1k_S1_L005_R1_001.fastq.gz
                               data/hgmm_1k/fastqs/hgmm_1k_S1_L005_R2_001.fastq.gz
[quant] will process sample 6: data/hgmm_1k/fastqs/hgmm_1k_S1_L006_R1_001.fastq.gz
                               data/hgmm_1k/fastqs/hgmm_1k_S1_L006_R2_001.fastq.gz
[quant] will process sample 7: data/hgmm_1k/fastqs/hgmm_1k_S1_L007_R1_001.fastq.gz
                               data/hgmm_1k/fastqs/hgmm_1k_S1_L007_R2_001.fastq.gz
[quant] will process sample 8: data/hgmm_1k/fastqs/hgmm_1k_S1_L008_R1_001.fastq.gz
                               data/hgmm_1k/fastqs/hgmm_1k_S1_L008_R2_001.fastq.gz
[quant] finding pseudoalignments for the reads ...
 done
[quant] processed 63,252,296 reads, 52,281,921 reads pseudoaligned

The quantification run in around 2 min on 10 threads and creates 4 files:

ls results/hgmm_1k/

• matrix.ec: A text file with two columns. The first column is the 0 based index of equivalence classes. The second column is the set of transcripts (denoted by 0 based index based on order of appearance in the transcriptome fasta file) present in the corresponding equivalence class.
• output.bus: The data represented in bus format. This is a binary file, so don’t use something like read.table to read it into R.
• run_info.json: Information about the call to kallisto bus, including the command used, number and percentage of reads pseudoaligned, version of kallisto used, etc.
• transcript.txt: A text file with one column, which is the transcripts present in the data, in the same order as in the transcriptome fasta file.

With the current 3’ sequencing protocols, most people are only interested in how many UMIs per gene per cell we can quantify. We will quantify this from the bus output, and to do so, we need to find which gene corresponds to each transcript. The t2g.txt file is a transcript-to-gene mapping, which allows us to have access to genes counts.

dir.create("results/hgmm_1k", recursive = T)
tr2g <- transcript2gene(
  fasta_file = c("data/hgmm_1k/hs_cdna.fa.gz", "./data/hgmm_1k/mm_cdna.fa.gz"),
  kallisto_out_path = "results/hgmm_1k/"
)
head(tr2g)
save_tr2g_bustools(tr2g, "results/hgmm_1k/tr2g_hgmm.tsv")

The /data/share/MADT/TP_lmodolo/10xv2_whitelist.txt contains all the barcodes known to be present in the kit v2 is provided by 10x. We are ready to make the gene count matrix. First, bustools runs barcode error correction on the bus file. Then, the corrected bus file is sorted by barcode, UMI, and equivalence classes. Then the UMIs are counted and the counts are collapsed into gene level.

mkdir -p results/hgmm_1k/genecount tmp
bustools correct -w data/whitelist_v2.txt -p results/hgmm_1k/output.bus | \
  bustools sort -T tmp/ -t 4 -p - | \
  bustools count -o results/hgmm_1k/genecount/genes -g results/hgmm_1k/tr2g_hgmm.tsv \
  -e results/hgmm_1k/matrix.ec -t results/hgmm_1k/transcripts.txt --genecounts

Found 737280 barcodes in the whitelist
Processed 52281921 BUS records
In whitelist = 50825261
Corrected    = 348752
Uncorrected  = 1107908

The run logs information is accessible in json format in the results/hmm_1k/ directory.

You can check how to load the data in R with the single-cell experiment or direclty go to the quality control of the data