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.

	#!/usr/bin/env python
	'''
	Automatically estimate insert size of the paired-end reads for a given SAM/BAM file.
	Usage: getinsertsize.py <SAM file> or samtools view <BAM file> | getinsertsize.py -
	Author: Wei Li
	Copyright (c) <2015> <Wei Li>
	Permission is hereby granted, free of charge, to any person obtaining a copy
	of this software and associated documentation files (the "Software"), to deal
	in the Software without restriction, including without limitation the rights
	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
	copies of the Software, and to permit persons to whom the Software is
	furnished to do so, subject to the following conditions:
	The above copyright notice and this permission notice shall be included in
	all copies or substantial portions of the Software.
	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
	THE SOFTWARE.
	'''
	from __future__ import print_function
	import sys;
	import pydoc;
	import os;
	import re;
	import fileinput;
	import math;
	import argparse;
	parser=argparse.ArgumentParser(description='Automatically estimate the insert size of the paired-end reads for a given SAM/BAM file.');
	parser.add_argument('SAMFILE',type=argparse.FileType('r'),help='Input SAM file (use - from standard input)');
	parser.add_argument('--span-distribution-file','-s',type=argparse.FileType('w'),help='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.');
	parser.add_argument('--read-distribution-file','-r',type=argparse.FileType('w'),help='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.');
	args=parser.parse_args();
	plrdlen={};
	plrdspan={};
	def getmeanval(dic,maxbound=-1):
	nsum=0; n=0;
	for (k,v) in dic.items():
	if maxbound!=-1 and k>maxbound:
	continue;
	nsum=nsum+k*v;
	n=n+v;
	meanv=nsum*1.0/n;
	nsum=0; n=0;
	for (k,v) in dic.items():
	if maxbound!=-1 and k>maxbound:
	continue;
	nsum=nsum+(k-meanv)*(k-meanv)*v;
	n=n+v;
	varv=math.sqrt(nsum*1.0/(n-1));
	return (meanv,varv);
	objmrl=re.compile('([0-9]+)M$');
	objmtj=re.compile('NH:i:(\d+)');
	nline=0;
	for lines in args.SAMFILE:
	field=lines.strip().split();
	nline=nline+1;
	if nline%1000000==0:
	print(str(nline/1000000)+'M...',file=sys.stderr);
	if len(field)<12:
	continue;
	try:
	mrl=objmrl.match(field[5]);
	if mrl==None: # ignore non-perfect reads
	continue;
	readlen=int(mrl.group(1));
	if readlen in plrdlen.keys():
	plrdlen[readlen]=plrdlen[readlen]+1;
	else:
	plrdlen[readlen]=1;
	if field[6]!='=':
	continue;
	dist=int(field[8]);
	if dist<=0: # ignore neg dist
	continue;
	mtj=objmtj.search(lines);
	# if mtj==None:
	# continue;
	# if int(mtj.group(1))!=1:
	# continue;
	#print(field[0]+' '+str(dist));
	if dist in plrdspan.keys():
	plrdspan[dist]=plrdspan[dist]+1;
	else:
	plrdspan[dist]=1;
	except ValueError:
	continue;
	#print(str(plrdlen));
	#print(str(plrdspan));
	# get the maximum value
	readlenval=getmeanval(plrdlen);
	print('Read length: mean '+str(readlenval[0])+', STD='+str(readlenval[1]));
	# print('Possible read lengths and their counts:');
	# print(str(plrdlen));
	if args.span_distribution_file is not None:
	for k in sorted(plrdspan.keys()):
	print(str(k)+'\t'+str(plrdspan[k]),file=args.span_distribution_file);
	if args.read_distribution_file is not None:
	for k in sorted(plrdlen.keys()):
	print(str(k)+'\t'+str(plrdlen[k]),file=args.read_distribution_file);
	if len(plrdspan)==0:
	print('No qualified paired-end reads found. Are they single-end reads?');
	else:
	maxv=max(plrdspan,key=plrdspan.get);
	spanval=getmeanval(plrdspan,maxbound=maxv*3);
	print('Read span: mean '+str(spanval[0])+', STD='+str(spanval[1]));
	# print('maxv:'+str(maxv));

view raw getinsertsize.py hosted with ❤ by GitHub

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.