Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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])
Add Comment
Please, Sign In to add comment