/usr/share/doc/python-pysam/tests/example.py is in python-pysam-tests 0.7.7-1ubuntu1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 | ## This script contains some example code
## illustrating ways to to use the pysam
## interface to samtools.
##
## The unit tests in the script pysam_test.py
## contain more examples.
##
import pysam
samfile = pysam.Samfile( "ex1.bam", "rb" )
print "###################"
# check different ways to iterate
print len(list(samfile.fetch()))
print len(list(samfile.fetch( "chr1", 10, 200 )))
print len(list(samfile.fetch( region="chr1:10-200" )))
print len(list(samfile.fetch( "chr1" )))
print len(list(samfile.fetch( region="chr1")))
print len(list(samfile.fetch( "chr2" )))
print len(list(samfile.fetch( region="chr2")))
print len(list(samfile.fetch()))
print len(list(samfile.fetch( "chr1" )))
print len(list(samfile.fetch( region="chr1")))
print len(list(samfile.fetch()))
print len(list(samfile.pileup( "chr1", 10, 200 )))
print len(list(samfile.pileup( region="chr1:10-200" )))
print len(list(samfile.pileup( "chr1" )))
print len(list(samfile.pileup( region="chr1")))
print len(list(samfile.pileup( "chr2" )))
print len(list(samfile.pileup( region="chr2")))
print len(list(samfile.pileup()))
print len(list(samfile.pileup()))
print "########### fetch with callback ################"
def my_fetch_callback( alignment ): print str(alignment)
samfile.fetch( region="chr1:10-200", callback=my_fetch_callback )
print "########## pileup with callback ################"
def my_pileup_callback( column ): print str(column)
samfile.pileup( region="chr1:10-200", callback=my_pileup_callback )
print "########### Using a callback object ###### "
class Counter:
mCounts = 0
def __call__(self, alignment):
self.mCounts += 1
c = Counter()
samfile.fetch( region = "chr1:10-200", callback = c )
print "counts=", c.mCounts
print "########### Calling a samtools command line function ############"
for p in pysam.mpileup( "-c", "ex1.bam" ):
print str(p)
print pysam.mpileup.getMessages()
print "########### Investigating headers #######################"
# playing arount with headers
samfile = pysam.Samfile( "ex3.sam", "r" )
print samfile.references
print samfile.lengths
print samfile.text
print samfile.header
header = samfile.header
samfile.close()
header["HD"]["SO"] = "unsorted"
outfile = pysam.Samfile( "out.sam", "wh",
header = header )
outfile.close()
|