Skip to content

Unix V: Comparing two quantification methods

Jump to bottom
Bioinformatics Core Facility Gothenburg edited this page Aug 31, 2017 · 12 revisions

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.

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?

Please don't run! Since it can take a long time!

Instead, copy the alignment file Fly.sam from /home/marcela/unix/files/.

The SAM format harbours 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 expression mean?

Now we need to count how many reads lie within each gene, so we know their expression level. There are of course, several ways of counting these reads, so let's compare two different approaches: htseq-count and bedtools coverage.

If you inspect the help pages for both tools (htseq-count is not installed in the server, you can have a look here), you will see that both need a feature file, a file where the positions of the genes we would like to quantify are specified. Copy the genes.gtf file from /home/marcela/unix/files/ and inspect it.

Since the program is not installed, copy the Fly_gene.htseq from /home/marcela/unix/files/

Q4. If you were to run htseq-count, would this expression be correct? if not, what would you change?

htseq-count bam -r pos -s yes -a 20 -q Fly.sort.bam genes.gtf > Fly_gene.htseq

Now, we will modify genes.gtf. While inspecting the file, you noticed that there is information on coding regions, exons, stop codons, etc. htseq-count was developed to quantify the exons in a gene and report back one single value: the gene count. bedtools coverage is not thought in the same way, so it will give you the counts of each line in your feature file, even though they belong to the same gene. Besides, you also noticed that there are genes with different transcripts, and for each transcript the gene is specified several times. To make this simpler, do the following:

  • extract only the exons (look at the 3rd column)
  • keep the columns containing chromosome, start, end, strand and the gene name (you may do some editing)
  • keep unique entries

your file should look like:

chr2L   21823   22687   -       Ir21a
chr2L   22743   22935   -       Ir21a
chr2L   22994   23873   -       Ir21a
chr2L   22998   24237   +       CR43609
chr2L   26766   26964   -       Cda5
chr2L   27053   27490   -       Cda5

Name your file exon_uniq.gtf.

Run genome coverage using the flag -counts and requiring strandedness. Save your results in a file named Fly_exon_uniq.counts

Q5. What is your command line?

Once you have inspected your count files (Fly_gene.htseq and Fly_exon_uniq.counts) create two lists, one for each gene, listing genes with more than 9999 of coverage, in other words, with a 5 digit coverage.

Q6. How many highly covered genes are shared by the two methods?

Next: Unix VI: Disease causing mutations

Previous: Unix IV: NGS data quality

Clone this wiki locally