Repository navigation
Unix V: Comparing two quantification methods
Course: HT17 Unix with applications to next generation sequencing data (SC00025)
Here we will compare two ways of generating expression gene levels from RNAseq data. For this exercise you need to know that the sequencing was paired-end and directional. The organism is Drosophila melanogaster, the fruit fly.
Copy Fly_R1.fastq.gz and Fly_2.fastq.gz from /home/marcela/unix/files, if you don't have them already.
First, we need to map our reads to the fly genome. Bowtie2 is one of many mappers that can be used to align sequencing reads to long reference sequences. If you type bowtie2 you will get a pretty big list of different parameters!.
Q1. What would be a minimum command line to run bowtie2 using Fly_R1.fastq.gz and Fly_R2.fastq.gz as input read files? You can use Bowtie2Index/genome as the Index filename prefix
Instead, copy the alignment file Fly.sam.gz from /home/marcela/unix/files/.
The SAM format harbors information on where and how each read is mapped in relation to the genome. Here is more complete explanation of the SAM format. Have look at it!
Q2. How many reads are mapped to the chromosome 3R?
The next step is to compress this data type into a BAM file, the binary format of SAM. For that we will be using samtools, a set of tools to manipulate alignments:
samtools view -Sb --threads 8 Fly.sam | samtools sort -o Fly.sort.bam - --threads 8
Q3. What does this cmd line mean?
Now we need to count how many reads lie within each gene, so we can quantify their expression level. There are of course, several ways/algorithms that count the number of reads mapped to a specific gene. For this exercise we provide you with the output files (counts) from two different counting approaches: A and B. How these algorithms work is out of the scope of this course, but ask if you want to know more!.
Copy the output file from the two approaches: Fly_A.counts and Fly_B.counts from /home/marcela/unix/files/ and inspect them.
Create two lists, one for each approach, listing genes with more than 9999 of coverage, in other words, with a 5 digit coverage. Compare them.
Q4. How many highly covered genes are shared by the two methods?
Previous: Unix IV: NGS data quality
Created by Marcela Dávila