RPKM for paired-end data

[Amative]

Hello,

I have aligned my RNA-Seq paired-end data to a reference using Bowtie2:
Bowtie2 -x Ref -1 pair1.fq -2 pair2.fq -S out.sam

In order to calculate RPKM we need to know the total number of mappable reads in the experiment which was obvious in the case of single end data. I am confused now which number to use as "mappable reads in the experiment" in the case of paired-end data. This is Bowtie2 alignment summary

19231140 reads; of these:
19231140 (100.00%) were paired; of these:
1790875 (9.31%) aligned concordantly 0 times
4367667 (22.71%) aligned concordantly exactly 1 time
13072598 (67.98%) aligned concordantly >1 times
----
1790875 pairs aligned concordantly 0 times; of these:
30148 (1.68%) aligned discordantly 1 time
----
1760727 pairs aligned 0 times concordantly or discordantly; of these:
3521454 mates make up the pairs; of these:
2645377 (75.12%) aligned 0 times
322397 (9.16%) aligned exactly 1 time
553680 (15.72%) aligned >1 times
93.12% overall alignment rate


Thank you in advance!

[davidblaney]

So to get 'mappable reads' I would look at the first output which is the read stats.

---
(1) 19231140 reads; of these:
(2) 19231140 (100.00%) were paired; of these:
(3) 1790875 (9.31%) aligned concordantly 0 times
(4) 4367667 (22.71%) aligned concordantly exactly 1 time
(5) 13072598 (67.98%) aligned concordantly >1 times
---

So in a plain English:

Out of 19231140 reads (1) all where paired (2), 1790875 did not align (3) and 17440265 (4+5) aligned 1 or more times.

therefore 17440265 reads were mapped.

The next 2 outputs just explain why the unmapped reads (3) didn't align I guess.

Hope that helps.

[Amative]

Thanks David, I have been looking for tools to analyze the output file from bowtie2 (SAM formatted) and found a couple (samtools flagstat and RSeQC):
The output from samtools flagstat is:

in total (QC-passed reads + QC-failed reads) 38462280 + 0
duplicates 0 + 0
mapped (93.12%:-nan%) 35816903 + 0
paired in sequencing 38462280 + 0
read1 19231140 + 0
read2 19231140 + 0
properly paired (90.69%:-nan%) 34880530 + 0
with itself and mate mapped 35157518 + 0
singletons (1.71%:-nan%) 659385 + 0
with mate mapped to a different chr 174452 + 0
with mate mapped to a different chr (mapQ>=5) 45969 + 0

you can see that the overall alignment reported by bowtie2 is what samtools flagstat called "mapped (93.12%:-nan%) 35816903 + 0".

Now, I am not sure wither to add the singletons (read pairs in which only one tag is successfully mapped) to the "mapped" or not, in order to use it to calculate RPKM

[davidblaney]

Yes I have tried those packages before. I also like Picards tool set for alignment metrics such as

CollectMultipleMetrics
CollectTargetedPcrMetrics
CollectRnaSeqMetrics

http://picard.sourceforge.net/command-line-overview.shtml

I guess the question to ask would be weather you could change tack a little and use FPKM.

Take a look at what cufflinks says here

What's the difference between FPKM and RPKM?

They're almost the same thing. RPKM stands for Reads Per Kilobase of transcript per Million mapped reads. FPKM stands for Fragments Per Kilobase of transcript per Million mapped reads. In RNA-Seq, the relative expression of a transcript is proportional to the number of cDNA fragments that originate from it. Paired-end RNA-Seq experiments produce two reads per fragment, but that doesn't necessarily mean that both reads will be mappable. For example, the second read is of poor quality. If we were to count reads rather than fragments, we might double-count some fragments but not others, leading to a skewed expression value. Thus, FPKM is calculated by counting fragments, not reads.

Also look at this tread :

http://seqanswers.com/forums/showthread.php?t=10332

David

[Amative]

Thanks Daivd, That was helpful. I guess I am going to try Cufflinks and see how it goes. As I already have the SAM/BAM files, so it should be easy to try it.

[davidblaney]

No worries,

This is a good paper to help with analysis using cummbeRbund also.

http://www.nature.com/nprot/journal/v7/n3/full/nprot.2012.016.html

[Amative]

Excellent, I will look at it as well