SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Forcing paired-end data mapped as single-end in SAM puggie Bioinformatics 1 03-16-2013 10:50 AM
Help with Illumina Paired-End Data adamba Bioinformatics 5 04-16-2012 12:36 PM
PRINSEQ and paired-end data Rockx Bioinformatics 1 03-10-2012 10:02 AM
Does Cufflinks support single-end and paired end data together ? ersenkavak Bioinformatics 1 10-22-2010 07:26 AM
Paired-end Illumina data mchaisso Bioinformatics 7 07-17-2008 11:52 AM

Reply
 
Thread Tools
Old 03-24-2013, 08:06 AM   #1
Amative
Member
 
Location: USA

Join Date: Dec 2011
Posts: 45
Default RPKM for paired-end data

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!
Amative is offline   Reply With Quote
Old 03-25-2013, 03:41 AM   #2
davidblaney
Member
 
Location: Oxford, UK

Join Date: Nov 2011
Posts: 17
Default

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.
davidblaney is offline   Reply With Quote
Old 03-25-2013, 07:09 AM   #3
Amative
Member
 
Location: USA

Join Date: Dec 2011
Posts: 45
Default

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
Amative is offline   Reply With Quote
Old 03-28-2013, 06:37 AM   #4
davidblaney
Member
 
Location: Oxford, UK

Join Date: Nov 2011
Posts: 17
Default

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

Quote:
CollectMultipleMetrics
CollectTargetedPcrMetrics
CollectRnaSeqMetrics
http://picard.sourceforge.net/comman...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

Quote:
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
davidblaney is offline   Reply With Quote
Old 03-28-2013, 07:47 AM   #5
Amative
Member
 
Location: USA

Join Date: Dec 2011
Posts: 45
Default

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.
Amative is offline   Reply With Quote
Old 03-31-2013, 04:24 AM   #6
davidblaney
Member
 
Location: Oxford, UK

Join Date: Nov 2011
Posts: 17
Default

No worries,

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

http://www.nature.com/nprot/journal/....2012.016.html
davidblaney is offline   Reply With Quote
Old 04-03-2013, 09:05 AM   #7
Amative
Member
 
Location: USA

Join Date: Dec 2011
Posts: 45
Default

Excellent, I will look at it as well
Amative is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 03:04 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO