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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
|
import os
import subprocess
import pysam
from TestUtils import BAM_DATADIR, force_str
def build_pileup_with_samtoolsshell(fn):
os.system("samtools mpileup {} 2> /dev/null | wc -l > /dev/null".format(fn))
return 2998
def build_pileup_with_samtoolspipe(fn):
FNULL = open(os.devnull, 'w')
with subprocess.Popen(["samtools", "mpileup", fn],
stdin=subprocess.PIPE,
stdout=subprocess.PIPE,
stderr=FNULL) as proc:
return len(proc.stdout.readlines())
def build_pileup_with_pysam(*args, **kwargs):
with pysam.AlignmentFile(*args, **kwargs) as inf:
return len(list(inf.pileup(stepper="samtools")))
def build_depth_with_samtoolsshell(fn):
os.system(
"samtools mpileup {} 2> /dev/null | awk '{{a += $4}} END {{print a}}' > /dev/null".format(fn))
return 107241
def build_depth_with_samtoolspipe(fn):
FNULL = open(os.devnull, 'w')
with subprocess.Popen(["samtools", "mpileup", fn],
stdin=subprocess.PIPE,
stdout=subprocess.PIPE,
stderr=FNULL) as proc:
data = [x.split() for x in proc.stdout.readlines()]
return [int(x[3]) for x in data]
def build_depth_with_filter_with_pysam(*args, **kwargs):
with pysam.AlignmentFile(*args, **kwargs) as inf:
return [x.get_num_aligned() for x in inf.pileup(stepper="samtools")]
def build_depth_with_pysam(*args, **kwargs):
with pysam.AlignmentFile(*args, **kwargs) as inf:
return [x.nsegments for x in inf.pileup(stepper="samtools")]
def build_query_bases_with_samtoolsshell(fn):
os.system("samtools mpileup {} 2> /dev/null | awk '{{a = a $5}} END {{print a}}' | wc -c > /dev/null".format(fn))
return 116308
def build_query_bases_with_samtoolspipe(fn, *args, **kwargs):
FNULL = open(os.devnull, 'w')
with subprocess.Popen(["samtools", "mpileup", fn] + list(args),
stdin=subprocess.PIPE,
stdout=subprocess.PIPE,
stderr=FNULL) as proc:
stdout = proc.stdout.read().decode()
return [x.split()[4] for x in stdout.splitlines()]
def build_query_bases_with_samtoolspysam(fn, *args):
return [x.split()[4] for x in pysam.samtools.mpileup(fn, *args).splitlines()]
def build_query_bases_with_pysam_pileups(*args, **kwargs):
total_pileup = []
with pysam.AlignmentFile(*args, **kwargs) as inf:
total_pileup = [
[r.alignment.query_sequence[r.query_position_or_next]
for r in column.pileups if r.query_position_or_next is not None]
for column in inf.pileup(stepper="samtools")]
return total_pileup
def build_query_qualities_with_pysam_pileups(*args, **kwargs):
total_pileup = []
with pysam.AlignmentFile(*args, **kwargs) as inf:
total_pileup = [
[r.alignment.query_qualities[r.query_position_or_next]
for r in column.pileups if r.query_position_or_next is not None]
for column in inf.pileup(stepper="samtools")]
return total_pileup
def build_query_bases_with_pysam(fn, *args, **kwargs):
total_pileup = []
with pysam.AlignmentFile(fn) as inf:
total_pileup = [column.get_query_sequences(
mark_ends=True, add_indels=True, mark_matches=True) for column in
inf.pileup(*args, **kwargs)]
return total_pileup
def build_query_names_with_pysam(*args, **kwargs):
total_pileup = []
with pysam.AlignmentFile(*args, **kwargs) as inf:
total_pileup = [column.get_query_names() for column in
inf.pileup(stepper="samtools")]
return total_pileup
def build_query_qualities_with_samtoolspipe(fn):
FNULL = open(os.devnull, 'w')
with subprocess.Popen(["samtools", "mpileup", fn],
stdin=subprocess.PIPE,
stdout=subprocess.PIPE,
stderr=FNULL) as proc:
data = [force_str(x).split()[5] for x in proc.stdout.readlines()]
return data
def build_query_qualities_with_pysam(*args, **kwargs):
total_pileup = []
with pysam.AlignmentFile(*args, **kwargs) as inf:
total_pileup = [column.get_query_qualities() for column in
inf.pileup(stepper="samtools")]
return total_pileup
def build_mapping_qualities_with_samtoolspipe(fn):
FNULL = open(os.devnull, 'w')
with subprocess.Popen(["samtools", "mpileup", "-s", fn],
stdin=subprocess.PIPE,
stdout=subprocess.PIPE,
stderr=FNULL) as proc:
data = [force_str(x).split()[6] for x in proc.stdout.readlines()]
return data
def build_mapping_qualities_with_pysam(*args, **kwargs):
total_pileup = []
with pysam.AlignmentFile(*args, **kwargs) as inf:
total_pileup = [column.get_mapping_qualities() for column in
inf.pileup(stepper="samtools")]
return total_pileup
def build_query_positions_with_samtoolspipe(fn):
FNULL = open(os.devnull, 'w')
with subprocess.Popen(["samtools", "mpileup", "-O", fn],
stdin=subprocess.PIPE,
stdout=subprocess.PIPE,
stderr=FNULL) as proc:
data = [list(map(int, force_str(x).split()[6].split(",")))
for x in proc.stdout.readlines()]
return data
def build_query_positions_with_pysam(*args, **kwargs):
total_pileup = []
with pysam.AlignmentFile(*args, **kwargs) as inf:
total_pileup = [column.get_query_positions() for column in
inf.pileup(stepper="samtools")]
return total_pileup
|
