Overlapping and non-Overlapping pair-end reads with Tophat

senpeng · 10-15-2011, 01:43 PM

we recently tried to use Tophat to align our illumina 2*100bp pair end data.

It seems that part of our data has a overlap from 30bp to 80bp, and some of the data doesn't have overlap.

if I have a mix of pair-ended reads such as some read pairs are overlapping some are not;
---> May I set up the inner distance to(-r option) -10 and then the std deviation(--mate-std-dev ) for the inner distance at 70 for example?
---> In that case, will Tophat look for inner distance from -80 to +60? (-80 meaning an overlap of 80

Thanks so much for your answer.

cedance · 10-16-2011, 05:19 AM

The way I go about figuring these parameters is this:
I use BWA and align the paired end reads and obtain PE.sam file. From that, I use picard tools to obtain only those reads that are aligned to the reference genome. From this sam file, I calculate the inner distance of all pairs. From the vector of all such inner distances, I compute the 1st and 3rd quantiles, Q1 and Q3 and then calculate the inter quartile range IQ = Q3-Q1. From here, I take all the values that fall between Q1 - 2*IQ to Q1 + 2*IQ. I compute the mean and standard deviation of all these values (this is similar to what BWA does) and provide it to tophat. It seems to work great. Most of my reads are mapped under the given inner distance and standard deviation and the other (few) reads with larger inner distance. There'll always be a few reads that'll have larger inner distance and I consider them outliers and discard in the computation of mean and SD.

Hope this helps.

upendra_35 · 10-16-2011, 11:26 AM

I normally do this way.....

I map the paired end reads using bowtie to obtain SAM files (Bowtie was chosen because of its speed). From the SAM files i extract the optimized insert_length and insert_length_sd. For obtaining the insert_length information from the SAM file i use the following perl script (Attached. The usage is quite simple). The output from this script is then put into R to calculate the mean and sd of the insert_length. You can use then use this information for tophat...

perl:
get_insert_sizes_from_sam.pl 300bp_pe_def.filt.sam > 300bp_pe_def.sizes

R:
sizes <- as.numeric(readLines("300bp_pe_def.sizes"))
mean(sizes, na.rm=TRUE)
median(sizes, na.rm=TRUE)
sd(sizes, na.rm=TRUE)

Hope this helps

senpeng · 10-16-2011, 07:33 PM

Thanks for your answers, cedance & upendra_35. Does that mean I have to run Bowtie or BWA first to get these numbers? I tried to apply different inner distances to test data(small of course) and didn't see much difference. How on earth will this parameter affect the final result?
For example, if I set the -r as 50 and SDV as 0, does that mean the program will not check paired of distance greater than 50?
Thanks again for your attention.

upendra_35 · 10-16-2011, 07:43 PM

Unfortunately you have to run either BWA or Bowtie to get the insert_length and insert_length SD info for tophat. If i am correct the parameters are quite important for correct mapping of the reads onto the reference. If you had tried a few more insert_lengths you would have seen the difference. The best way to compare the different lengths is to apply samtools flagstat on the BAM files and see how many reads mapped vs non mapped.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Overlapping paired end - tophat	wenhuang	Bioinformatics	25	04-29-2011 11:51 AM
Determine paired end overlapping	chariko	Bioinformatics	2	04-29-2011 12:52 AM
Questions about overlapping paired-end reads...	FredOnSeq	Illumina/Solexa	6	04-18-2011 06:19 PM
How to manage overlapping paired-end reads?	FredOnSeq	Bioinformatics	2	09-09-2010 02:27 AM
How do variant callers deal with overlapping paired end reads?	krobison	Bioinformatics	1	04-30-2010 12:58 PM

10-15-2011, 01:43 PM	#1
senpeng Member Location: phoenix Join Date: Sep 2011 Posts: 10	Overlapping and non-Overlapping pair-end reads with Tophat we recently tried to use Tophat to align our illumina 2*100bp pair end data. It seems that part of our data has a overlap from 30bp to 80bp, and some of the data doesn't have overlap. if I have a mix of pair-ended reads such as some read pairs are overlapping some are not; ---> May I set up the inner distance to(-r option) -10 and then the std deviation(--mate-std-dev ) for the inner distance at 70 for example? ---> In that case, will Tophat look for inner distance from -80 to +60? (-80 meaning an overlap of 80 Thanks so much for your answer.

10-16-2011, 05:19 AM	#2
cedance Senior Member Location: Germany Join Date: Feb 2011 Posts: 108	The way I go about figuring these parameters is this: I use BWA and align the paired end reads and obtain PE.sam file. From that, I use picard tools to obtain only those reads that are aligned to the reference genome. From this sam file, I calculate the inner distance of all pairs. From the vector of all such inner distances, I compute the 1st and 3rd quantiles, Q1 and Q3 and then calculate the inter quartile range IQ = Q3-Q1. From here, I take all the values that fall between Q1 - 2IQ to Q1 + 2IQ. I compute the mean and standard deviation of all these values (this is similar to what BWA does) and provide it to tophat. It seems to work great. Most of my reads are mapped under the given inner distance and standard deviation and the other (few) reads with larger inner distance. There'll always be a few reads that'll have larger inner distance and I consider them outliers and discard in the computation of mean and SD. Hope this helps.

10-16-2011, 07:33 PM	#4
senpeng Member Location: phoenix Join Date: Sep 2011 Posts: 10	Thanks for your answers, cedance & upendra_35. Does that mean I have to run Bowtie or BWA first to get these numbers? I tried to apply different inner distances to test data(small of course) and didn't see much difference. How on earth will this parameter affect the final result? For example, if I set the -r as 50 and SDV as 0, does that mean the program will not check paired of distance greater than 50? Thanks again for your attention.

10-16-2011, 07:43 PM	#5
upendra_35 Senior Member Location: USA Join Date: Apr 2010 Posts: 102	Unfortunately you have to run either BWA or Bowtie to get the insert_length and insert_length SD info for tophat. If i am correct the parameters are quite important for correct mapping of the reads onto the reference. If you had tried a few more insert_lengths you would have seen the difference. The best way to compare the different lengths is to apply samtools flagstat on the BAM files and see how many reads mapped vs non mapped.