RNA sequencing data analysis - alignment and reads counting using cufflinks

This analysis was performed using R (ver. 3.1.0).


RNA-Seq is a valuable experiment for quantifying both the types and the amount of RNA molecules in a sample. The aims of this article is to show you :

  • How to align the reads using tophat2?
  • How to count the number of reads per features using cufflinks?
  • How to visualize bam file using IGV ?

Download data

We will work with the Hammer et al dataset, as prepared by the ReCount website.

The Hammer et al paper: http://www.ncbi.nlm.nih.gov/pubmed?term=20452967 The GEO page: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE20895 The sample I will align: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM539553 which points to the SRA: ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByExp/sra/SRX%2FSRX020%2FSRX020088/SRR042499/

download the file

#download the file
wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByExp/sra/SRX%2FSRX020%2FSRX020088/SRR042499/SRR042499.sra
#convert to FastQ
fastq-dump SRR042499.sra

Download reference genome from tophat website

wget wget [url=ftp://igenome:G3nom3s4u@ussd-ftp.illumina.com/Rattus_norvegicus/Ensembl/RGSC3.4/Rattus_norvegicus_Ensembl_RGSC3.4.tar.gz]ftp://igenome:G3nom3s4u@ussd-ftp.illumina.com/Rattus_norvegicus/Ensembl/RGSC3.4/Rattus_norvegicus_Ensembl_RGSC3.4.tar.gz[/url]
tar zxvf Rattus_norvegicus_Ensembl_RGSC3.4.tar.gz

The genome was downloaded from Illumina iGenomes

Align the reads

The tophat call to align the reads (this take account for splice read) :

tophat2 -o tophat_out -p 10 /path/to/Rattus_norvegicus/Ensembl/RGSC3.4/Sequence/Bowtie2Index/genome SRR042499.fastq

-o tophat_out : Output directory -p 10 : number of threads to use

Tophat2 is an aligner which takes account splice reads.

Tophat2 output

Tophat2 output

Use IGV to visualize

I’m going to just show you what these look like in IGV. So I’ve downloaded the accepted hits.bam file and the transcripts .GTF file to my local directory. So I then use this line here to– it’s a grep line which removes any lines which match this. I wanted to remove lines which had FPKM 0.

grep -v ‘FPKM “0.0000000000”’ transcripts.gtf | less

And then I downloaded the chromosome 1 FASTA file and the rat GTF file so that I could look at these in IGV. So now here’s my IGV window where I’ve loaded in the coverage. And this line here is the ensemble GTF file and below it is the transcripts.gtf file from Cufflinks.

So they look very similar but we can zoom in and try to find some differences.

For visualizing:

ftp://ftp.ensembl.org/pub/release-69/fasta/rattus_norvegicus/dna/Rattus_norvegicus.RGSC3.4.69.dna.chromosome.1.fa.gz ftp://ftp.ensembl.org/pub/release-69/gtf/rattus_norvegicus/Rattus_norvegicus.RGSC3.4.69.gtf.gz

Licence: Licence References: https://github.com/genomicsclass

Enjoyed this article? I’d be very grateful if you’d help it spread by emailing it to a friend, or sharing it on Twitter, Facebook or Linked In.

Show me some love with the like buttons below... Thank you and please don't forget to share and comment below!!
Avez vous aimé cet article? Je vous serais très reconnaissant si vous aidiez à sa diffusion en l'envoyant par courriel à un ami ou en le partageant sur Twitter, Facebook ou Linked In.

Montrez-moi un peu d'amour avec les like ci-dessous ... Merci et n'oubliez pas, s'il vous plaît, de partager et de commenter ci-dessous!

This page has been seen 14689 times