tests/example.py

   1 ## This script contains some example code
   2 ## illustrating ways to to use the pysam
   3 ## interface to samtools.
   4 ##
   5 ## The unit tests in the script pysam_test.py
   6 ## contain more examples.
   7 ##
   8
   9 import pysam
  10
  11 samfile = pysam.Samfile( "ex1.bam", "rb" )
  12
  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()))
  26
  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()))
  35
  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 )
  39
  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 )
  43
  44
  45 print "########### Using a callback object ###### "
  46
  47 class Counter:
  48     mCounts = 0
  49     def __call__(self, alignment):
  50         self.mCounts += 1
  51
  52 c = Counter()
  53 samfile.fetch( region = "chr1:10-200", callback = c )
  54 print "counts=", c.mCounts
  55
  56 print "########### Calling a samtools command line function ############"
  57
  58 for p in pysam.mpileup( "-c", "ex1.bam" ):
  59     print str(p)
  60
  61 print pysam.mpileup.getMessages()
  62
  63 print "########### Investigating headers #######################"
  64
  65 # playing arount with headers
  66 samfile = pysam.Samfile( "ex3.sam", "r" )
  67 print samfile.references
  68 print samfile.lengths
  69 print samfile.text
  70 print samfile.header
  71 header = samfile.header
  72 samfile.close()
  73
  74 header["HD"]["SO"] = "unsorted"
  75 outfile = pysam.Samfile( "out.sam", "wh",
  76                          header = header )
  77
  78 outfile.close()
  79