indelratio.py

import pysam
import sys
import matplotlib.pyplot as plt
import numpy as np

def main(bam):
	samfile = pysam.AlignmentFile(bam)
	ratios = np.array([indelratio(read.cigartuples) for read in samfile.fetch()])
	n, bins, patches = plt.hist(ratios, bins=[i/10 for i in range(1,50, 1)])
	plt.xlabel('indel ratio')
	plt.ylabel('reads')
	plt.savefig("Indelratio.png")
	print(np.median(ratios))
	

def indelratio(cigartuples):
	ins = sum([oplen for op, oplen in cigartuples if op == 1])
	dels = sum([oplen for op, oplen in cigartuples if op == 2]) or 1
	return ins/dels


if __name__ == "__main__":
	main(sys.argv[1])