Question: Estimating Insert Size From Paired End Data.
gravatar for geek_y
6.6 years ago by
geek_y11k wrote:


I have paired end data from illumina. To estimate the insert size in silico ( from scratch ), I have aligned the reads as single end reads to the genome ( mouse ). Now I have the two alignment files ( SAM/BAM). I would like to estimate the insert size ( distribution plot ) from the two SAM/BAM files.

The existing tools will look for the "=" field in the SAM format to consider it as the corresponding mate but not read name. Any one is aware of any tool that estimates the insert size based on reads names from the two sam files ?


picard paired-end alignment • 22k views
ADD COMMENTlink modified 3.0 years ago by jmodlis0 • written 6.6 years ago by geek_y11k

I have written a script that calculates insert size from single end alignments. Thanks all for your help.

ADD REPLYlink written 6.5 years ago by geek_y11k
gravatar for Ashutosh Pandey
6.6 years ago by
Ashutosh Pandey12k wrote:

Estimating insert size the easy way:

  1. Align the paired-end data in a combined way against the reference genome. By combined, I mean you should provide both the fastq files to the aligner.
  2. Extract information from the ninth column of the SAM file (TLEN). To be more accurate, you can only select the number in the ninth column of those paired-end read pairs that are uniquely aligned against the genome.
  3. Calculate the mean of the distribution of TLEN numbers generated in step 2.
  4. Subtract length of a read (For example 75 bp or 100 bp) from the mean to get insert size.

For aligners like BWA, you need not to give the insert size. It does it for you automatically. See below paragraph from BWA manual to know how it does it:

Estimating Insert Size Distribution

BWA estimates the insert size distribution per 2561024 read pairs. It first collects pairs of reads with both ends mapped with a single-end quality 20 or higher and then calculates median (Q2), lower and higher quartile (Q1 and Q3). It estimates the mean and the variance of the insert size distribution from pairs whose insert sizes are within interval [Q1-2(Q3-Q1), Q3+2(Q3-Q1)]. The maximum distance x for a pair considered to be properly paired (SAM flag 0x2) is calculated by solving equation Phi((x-mu)/sigma)=x/Lp0, where mu is the mean, sigma is the standard error of the insert size distribution, L is the length of the genome, p0 is prior of anomalous pair and Phi() is the standard cumulative distribution function. For mapping Illumina short-insert reads to the human genome, x is about 6-7 sigma away from the mean. Quartiles, mean, variance and x will be printed to the standard error output.

Estimating insert size your way or hard way:

  1. Sort both the SAM files using queryname.
  2. Remove secondary alignments for a read so that read order and read indexes become same in both the SAM files.
  3. Now for every read, calculate absolute difference in their mapping position (Ignore reads that are mapped to different chromosomes
  4. Take the mean of the distribution and subtract the read length.
ADD COMMENTlink modified 9 months ago by RamRS30k • written 6.6 years ago by Ashutosh Pandey12k

Ashutosh Pandey.. Thanks. I have the idea how to do it. Just I am looking if there is any tool available.

ADD REPLYlink written 6.6 years ago by geek_y11k

I apologize, I thought you are a newbie. Well give it a try. It may work for you. I mean this feature may update the TLEN values for your alignments in the new bam file.

ADD REPLYlink written 6.6 years ago by Ashutosh Pandey12k

Dear Ashutosh,

When I take the two reads from both the files, the bed format looks like this:

68331263        68331364        HWI-1KL120:91:C0CBKACXX:1:1101:1420:2186_R1        42      +
68331437        68331536        HWI-1KL120:91:C0CBKACXX:1:1101:1420:2186_R2        42      -

if I take R1 position - R2 position, I end up in negative value. Is there any thing to do with strand information?

ADD REPLYlink modified 9 months ago by RamRS30k • written 6.6 years ago by geek_y11k

Why do you subtract the read length in the fourth step? Do the insert size include the read length? Insert Size And Fragment Size ?

ADD REPLYlink written 3.3 years ago by verne9120
gravatar for Irsan
6.6 years ago by
Irsan7.2k wrote:

Picard tools CollectInsertSizeMetrics

ADD COMMENTlink written 6.6 years ago by Irsan7.2k

Picard tools CollectInsertSizeMetrics takes only one bam file as input. If we merge the two files is not giving any output. The main problem is how the tool makes a decision whether the two reads are paired. Is it based on read names or the "=" tag in the SAM file.

ADD REPLYlink written 6.6 years ago by geek_y11k

I have never used this feature from Picard but I think MergeBAMAlignment should do the job. It can take seperate BAM files each representing first and second pair.

ADD REPLYlink modified 6.6 years ago • written 6.6 years ago by Ashutosh Pandey12k

You may be using MergeSamFiles I dont think it will update or calculate the value of the TLEN or ninth column in the new BAM file.

ADD REPLYlink written 6.6 years ago by Ashutosh Pandey12k
gravatar for Rayan Chikhi
6.4 years ago by
Rayan Chikhi1.4k
France, Lille, CNRS
Rayan Chikhi1.4k wrote:

This script enables to get an insert size estimation very quickly (based on BWA's intermediate alignment results)

ADD COMMENTlink modified 9 months ago by RamRS30k • written 6.4 years ago by Rayan Chikhi1.4k
gravatar for jmodlis
3.0 years ago by
Durham, NC, USA
jmodlis0 wrote:

I know this post is really old, but if anyone is still looking for the answer to this (like I was) if you use SOAPdenovo2, the average insert size will be estimated during the assembly -- the average insert size is outputted into the log file, along with the standard deviation. The insert size estimate will change slightly depending on the assembly, but it should at least give you an idea.

ADD COMMENTlink written 3.0 years ago by jmodlis0
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1287 users visited in the last hour