1 ## This script contains some example code
2 ## illustrating ways to to use the pysam
3 ## interface to samtools.
5 ## The unit tests in the script pysam_test.py
6 ## contain more examples.
11 samfile = pysam.Samfile( "ex1.bam", "rb" )
13 print "###################"
14 # check different ways to iterate
15 print len(list(samfile.fetch()))
16 print len(list(samfile.fetch( "chr1", 10, 200 )))
17 print len(list(samfile.fetch( region="chr1:10-200" )))
18 print len(list(samfile.fetch( "chr1" )))
19 print len(list(samfile.fetch( region="chr1")))
20 print len(list(samfile.fetch( "chr2" )))
21 print len(list(samfile.fetch( region="chr2")))
22 print len(list(samfile.fetch()))
23 print len(list(samfile.fetch( "chr1" )))
24 print len(list(samfile.fetch( region="chr1")))
25 print len(list(samfile.fetch()))
27 print len(list(samfile.pileup( "chr1", 10, 200 )))
28 print len(list(samfile.pileup( region="chr1:10-200" )))
29 print len(list(samfile.pileup( "chr1" )))
30 print len(list(samfile.pileup( region="chr1")))
31 print len(list(samfile.pileup( "chr2" )))
32 print len(list(samfile.pileup( region="chr2")))
33 print len(list(samfile.pileup()))
34 print len(list(samfile.pileup()))
36 print "########### fetch with callback ################"
37 def my_fetch_callback( alignment ): print str(alignment)
38 samfile.fetch( region="chr1:10-200", callback=my_fetch_callback )
40 print "########## pileup with callback ################"
41 def my_pileup_callback( column ): print str(column)
42 samfile.pileup( region="chr1:10-200", callback=my_pileup_callback )
45 print "########### Using a callback object ###### "
49 def __call__(self, alignment):
53 samfile.fetch( region = "chr1:10-200", callback = c )
54 print "counts=", c.mCounts
56 print "########### Calling a samtools command line function ############"
58 for p in pysam.mpileup( "-c", "ex1.bam" ):
61 print pysam.mpileup.getMessages()
63 print "########### Investigating headers #######################"
65 # playing arount with headers
66 samfile = pysam.Samfile( "ex3.sam", "r" )
67 print samfile.references
71 header = samfile.header
74 header["HD"]["SO"] = "unsorted"
75 outfile = pysam.Samfile( "out.sam", "wh",