Module 1 · From reads to counts

Videos: T1 What does RNA-seq measure? and T2 From reads to counts. Script: scripts/01_consolidate_kallisto.R

Goal. Turn FASTQ files into a gene × sample table of estimated counts, and understand every step.

Part A (optional, needs kallisto): quantify one sample yourself

You need three things:

• Reads. The FASTQ files for one sample from PRJNA412001. The samples are paired-end, so there are two files (_1 and _2) per sample, with about 40 million read pairs each. The files are large, so start the download early.
• Transcriptome. The mouse cDNA from Ensembl: Mus_musculus.GRCm38.cdna.all.fa.gz. The provided results used GRCm38; with a newer release your numbers will differ slightly.

* Transcript → gene map. From BioMart. Choose Ensembl Genes → Mouse genes (GRCm38). Under Filters → Gene, choose protein_coding for both Gene type and Transcript type. Under Attributes, choose Transcript stable ID version, Gene stable ID and Gene name. Export the result as a TSV. This is the file data/mart_export.txt.

# build the index once per transcriptome (it chops every transcript into k-mers; k = 31 by default)
kallisto index -i mm_grcm38.idx Mus_musculus.GRCm38.cdna.all.fa.gz
# quantify one paired-end sample
kallisto quant -i mm_grcm38.idx -o SRR6068402 SRR6068402_1.fastq.gz SRR6068402_2.fastq.gz

Open SRR6068402/abundance.tsv and SRR6068402/run_info.json.

Q1.1 What do the columns length, eff_length, est_counts and tpm mean? Why is est_counts not a whole number?

Q1.2 In run_info.json, what fraction of reads pseudo-aligned (p_pseudoaligned)? Compare it with the other seven samples in data/kallisto/*/run_info.json. What could explain a low value?

Part B: combine 8 samples into gene-level tables

Run scripts/01_consolidate_kallisto.R. It:

1. reads the eight abundance.tsv files into transcript × sample tables of counts and TPM;
2. keeps the protein coding transcripts listed in mart_export.txt, and sums them per gene with rowsum();
3. saves results/count_gene.txt and results/tpm_gene.txt, and checks them against the provided tables.

Q1.3 How many transcripts are in the kallisto index? How many protein coding transcripts and genes do we keep?

Q1.4 The gene Pkm has several transcripts. How many of them get almost all of its reads? The paper behind this data set is about a switch between two Pkm isoforms. Could you see such a switch in the gene-level table? Which table would you use instead?

Q1.5 (theory, video T2) A read is compatible with transcripts {tx1, tx2} by its first k-mer, {tx1} by its middle k-mers, and {tx1, tx3} by its last k-mer. Which transcripts is the read compatible with, and why do we take the intersection, not the union?

Q1.6 (theory, video T2) 30 reads fit only tx1, 10 fit only tx2, and 40 fit both. The transcripts have equal length. Do two rounds of EM by hand, starting from 50 : 50. Where does it converge?