Entering edit mode
7.2 years ago
ammarsabir15
▴
70
Using the code below I am trying to plot the values for coverage of a BAM file,
#!/usr/bin/python3
import matplotlib.pyplot as plt
import pysam,sys
bamfile = sys.argv[1] # Input bam file
samfile = pysam.Samfile(bamfile, "rb" )
samfile.pileup(max_depth = 1000000)
X = []
for pileupcolumn in samfile.pileup():
X.append(pileupcolumn.n)
Y = range(len(X))
plt.title(bamfile)
plt.xlabel('Bases')
plt.ylabel('Coverage')
plt.plot(Y, X,c = 'g')
plt.grid(True, lw = 2, ls = '-', c = '.75')
plt.xlim(0, 16569)
plt.axhline(y=230, color='r', linestyle='--')
plt.savefig('Coverage_plot11c.png', c = 'k')
But the code is not working right and giving following plot for which coverage seems to be limited to 8000 which can be seen here Coverage plot
But the exact values for coverage are much higher than those plotted by using pysam, I have calculted coverage for BAM file using bamtools and part of that file is shown below,
gi|251831106|ref|NC_012920.1| 0 13130
gi|251831106|ref|NC_012920.1| 1 13207
gi|251831106|ref|NC_012920.1| 2 13344
gi|251831106|ref|NC_012920.1| 3 13366
gi|251831106|ref|NC_012920.1| 4 13377
gi|251831106|ref|NC_012920.1| 5 13465
gi|251831106|ref|NC_012920.1| 6 13470
gi|251831106|ref|NC_012920.1| 7 13601
gi|251831106|ref|NC_012920.1| 8 13886
gi|251831106|ref|NC_012920.1| 9 13912
gi|251831106|ref|NC_012920.1| 10 14013
gi|251831106|ref|NC_012920.1| 11 14015
gi|251831106|ref|NC_012920.1| 12 14048
gi|251831106|ref|NC_012920.1| 13 14057
gi|251831106|ref|NC_012920.1| 14 14061
gi|251831106|ref|NC_012920.1| 15 14076
gi|251831106|ref|NC_012920.1| 16 14275
gi|251831106|ref|NC_012920.1| 17 14311
gi|251831106|ref|NC_012920.1| 18 14576
gi|251831106|ref|NC_012920.1| 19 14579
gi|251831106|ref|NC_012920.1| 20 14763
gi|251831106|ref|NC_012920.1| 21 14772
gi|251831106|ref|NC_012920.1| 22 14813
gi|251831106|ref|NC_012920.1| 23 14819
gi|251831106|ref|NC_012920.1| 24 14829
gi|251831106|ref|NC_012920.1| 25 14874
gi|251831106|ref|NC_012920.1| 26 14904
gi|251831106|ref|NC_012920.1| 27 15016
gi|251831106|ref|NC_012920.1| 28 15061
gi|251831106|ref|NC_012920.1| 29 15067
gi|251831106|ref|NC_012920.1| 30 15151
gi|251831106|ref|NC_012920.1| 31 15152
gi|251831106|ref|NC_012920.1| 32 15199
gi|251831106|ref|NC_012920.1| 33 15192
gi|251831106|ref|NC_012920.1| 34 15149
gi|251831106|ref|NC_012920.1| 35 15216
gi|251831106|ref|NC_012920.1| 36 14942
gi|251831106|ref|NC_012920.1| 37 14920
gi|251831106|ref|NC_012920.1| 38 14922
gi|251831106|ref|NC_012920.1| 39 14842
gi|251831106|ref|NC_012920.1| 40 15065
gi|251831106|ref|NC_012920.1| 41 14990
gi|251831106|ref|NC_012920.1| 42 14998
gi|251831106|ref|NC_012920.1| 43 14915
gi|251831106|ref|NC_012920.1| 44 14894
gi|251831106|ref|NC_012920.1| 45 14884
gi|251831106|ref|NC_012920.1| 46 14857
gi|251831106|ref|NC_012920.1| 47 14840
Is there a way that I can change this code to plot the exact values for coverage of BAM file ??
Thanks a lot, now it plots exact values.