SEQanswers Overlapping and non-Overlapping pair-end reads with Tophat
 User Name Remember Me? Password
 Register FAQ Members List Calendar Search Today's Posts Mark Forums Read

 Similar Threads Thread Thread Starter Forum Replies Last Post wenhuang Bioinformatics 25 04-29-2011 10:51 AM chariko Bioinformatics 2 04-28-2011 11:52 PM FredOnSeq Illumina/Solexa 6 04-18-2011 05:19 PM FredOnSeq Bioinformatics 2 09-09-2010 01:27 AM krobison Bioinformatics 1 04-30-2010 11:58 AM

 10-15-2011, 12: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, 04: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 - 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.
10-16-2011, 10:26 AM   #3
upendra_35
Senior Member

Location: USA

Join Date: Apr 2010
Posts: 102

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:
mean(sizes, na.rm=TRUE)
median(sizes, na.rm=TRUE)
sd(sizes, na.rm=TRUE)

Hope this helps
Attached Files
 get_insert_sizes_from_bowtie_aln.pl (1.7 KB, 126 views)

 10-16-2011, 06: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, 06: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.