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?

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.

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.

Calculate the mean of the distribution of TLEN numbers generated in step 2.

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:

Sort both the SAM files using queryname.

Remove secondary alignments for a read so that read order and read indexes become same in both the SAM files.

Now for every read, calculate absolute difference in their mapping position (Ignore reads that are mapped to different chromosomes

Take the mean of the distribution and subtract the read length.

Picard tools CollectInsertSizeMetrics takes only one bam file as input. If we merge the two files also..it 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.

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.

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