Overlapping and non-Overlapping pair-end reads with Tophat

[senpeng]

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]

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]

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]

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]

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.