Estimating paired-end read insert length from SAM/BAM files

I wrote a single Python script to estimate the paired-end read insert length (or fragment length) from read mapping information (i.e., SAM/BAM files). The algorithm is simple: check the TLEN field in the SAM format, throw out pair-end reads whose pairs are too far away, and use them to estimate the mean and variance of the insert length.

This script is also able to provide a detailed distribution of read length and read span for your convenience. Please refer to the detailed usage below.

This script is distributed in GitHub now.

Usage:

getinsertsize.py [ SAM file | -]

or

samtools view [ BAM file ] | getinsertsize.py - 


Detailed Usage:


usage: getinsertsize.py [-h] [--span-distribution-file SPAN_DISTRIBUTION_FILE]
                        [--read-distribution-file READ_DISTRIBUTION_FILE]
                        SAMFILE

Automatically estimate the insert size of the paired-end reads for a given
SAM/BAM file.

positional arguments:
  SAMFILE               Input SAM file (use - from standard input)

optional arguments:
  -h, --help            show this help message and exit
  --span-distribution-file SPAN_DISTRIBUTION_FILE, -s SPAN_DISTRIBUTION_FILE
                        Write the distribution of the paired-end read span
                        into a text file with name SPAN_DISTRIBUTION_FILE.
                        This text file is tab-delimited, each line containing
                        two numbers: the span and the number of such paired-
                        end reads.
  --read-distribution-file READ_DISTRIBUTION_FILE, -r READ_DISTRIBUTION_FILE
                        Write the distribution of the paired-end read length
                        into a text file with name READ_DISTRIBUTION_FILE.
                        This text file is tab-delimited, each line containing
                        two numbers: the read length and the number of such
                        paired-end reads.

Sample output:

Read length: mean 90.6697303194, STD=15.9446036414

Possible read length and their counts:

{108: 43070882, 76: 50882326}

Read span: mean 165.217445903, STD=32.8914834802

Note: If the SAM/BAM file size is too large, it is accurate enough to estimate based on a few reads (like 1 millioin). In this case, you can run the script as follows:

head -n 1000000 [ SAM file ] |  getinsertsize.py -

or

samtools view [ BAM file ] | head -n 1000000 | getinsertsize.py -

Note: According to the SAM definition, the read span "equals the number of bases from the leftmost mapped base to the rightmost mapped base". This span is the distance between two reads in a paired-end read PLUS 2 times read length. Read span is different from the "mate-inner-distance" in Tophat (-r option), which measures only the distance between two reads in a paired-end read.

RNA-Seq Read Simulator

RNA-Seq is now a common protocol to study the expression of genes or transcripts. For research purposes, there are many simulators to simulate RNA-Seq reads, like Flux Simulator. But many times I found it hard to use because: 1) there are complicated parameters, 2) it requires large memory and 3) it crushes frequently. So I wrote a few scripts to generate simulated RNA-Seq reads, and publish them in a package "RNASeqReadSimulator".

 RNASeqReadSimulator is a set of scripts generating simulated RNA-Seq reads. It provides users a simple tool to generate RNA-Seq reads for research purposes, and a framework to allow experienced users to expand functions. RNASeqReadSimulator offers the following features:

1. It allows users to randomly assign expression levels of transcripts and generate simulated single-end or paired-end RNA-Seq reads. 
2. It is able to generate RNA-Seq reads that have a specified positional bias profile. 
3. It is able to simulate random read errors from sequencing platforms. 
4. The simulator consists of a few simple Python scripts. All scripts are command line driven, allowing users to invoke and design more functions.

The webpage of RNASeqReadSimulator is here. You can find the source code in GitHub.