exercise 2 Quantify with the RSEM program

RSEM depends on an existing annotation and will only scores transcripts that are present in the given annotation file. We will compare the alignments produced by RSEM and tophat and this will become clear.

The first step is to prepare the transcript set that we will quantify. We selected the UCSC genes which is a very comprehensive, albeit a bit noisy dataset. As with all the data in this activity we will only use the subset of the genes that map to the genome regions we are using.

2.1 Calculate expression

RSEM now is ready to align and then attempt to perform read assignment and counting for each isoform in the file provided above. You must process each one of the 6 libraries:

cd ~/RNASeqWS/transcriptomics

mkdir rsem

rsem-calculate-expression --paired-end -p 2 \
 --output-genome-bam fastq.quantification/control_rep1.1.fq \
fastq.quantification/control_rep1.2.fq genome.quantification/mm10.rsem rsem/ctrl1.rsem

And similarly for each of the other 5 libraries

/project/umw_biocore/bin/class/runRSEM.bash

You should take the time to familiarize yourself with the output

2.2 Create consolidated table

In the bin directory (we provide a simple script to take all the independent RSEM output and combine it into a single table, which is then useful for inspection and ready for differential gene expression analysis.

To find out what the script does you may type the following command in the transcriptomics directory

cd ~/RNASeqWS/transcriptomics

perl bin/rsem.to.table.pl -help

We will generate two tables with isoform and gene level expression:

perl bin/rsem.to.table.pl -out rsem.gene.summary.count.txt -indir rsem -gene_iso genes -quantType expected_count
perl bin/rsem.to.table.pl -out rsem.isoforms.summary.count.txt -indir rsem -gene_iso isoforms -quantType expected_count

You should now take the time to inspect these tables, find genes that look affected by the experiment.

2.3 Visualize the raw data:

We are now ready to look at the data. Download the IGV program from the IGV site and transfer your files from the hpcc to your laptop. In general you will not transfer the files. There are ways to access the files from your dekstop or laptop directly at the HPCC, we'll cover this later. Since the files are very small transfer time will not be a problem. Use Filezila or the scp command if you are on a mac, to copy the files onto your laptop/desktop, please ensure you note to which directory you transfer the files as later you will need to load them into IGV.

Launch the IGV browser, and use the File -> load to load the files onto the browser. Load only one or two .bam files to begin with. Don't forget to choose mm10 as a genome build from the corner

A few genes are good examples of differentially expressed genes. For example the whole region around the key Fgf21 gene is upregulated in experiment vs controls, while the gene Crebbp is downregulated in experiments vs controls. To point your browser to either gene just type or copy the name of the gene in the location box at the top.

We will revisit these genes below when we do differential gene expression.

mkdir ~/RNASeqWS/transcriptomics/tophat

cd ~/RNASeqWS/transcriptomics

tophat2 --library-type fr-firststrand --segment-length 20 -G  genome.quantification/ucsc.gtf \
-o  tophat/ctrl1 genome.quantification/mm10 fastq.quantification/control_rep1.1.fq \
fastq.quantification/control_rep1.2.fq

/project/umw_biocore/bin/class/runTophat.bash

mv tophat/ctrl1/accepted_hits.bam tophat/ctrl1.bam

/project/umw_biocore/bin/class/moveBamFiles.bash 

cd ~/RNASeqWS/transcriptomics/tophat

samtools sort ctrl1.bam ctrl1.sorted
samtools index ctrl1.sorted.bam

/project/umw_biocore/bin/class/makeSortIndex.bash

module load R/3.0.1
Rscript  bin/rsem.quant.r

colData = as.data.frame((colnames(d)));
colData[1,2] = "control";
colData[2,2] = "expr";
colData[3,2] = "expr";
colData[4,2] = "control";
colData[5,2] = "expr";
colData[6,2] = "control";