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
79
|
## 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()
|
