-
Notifications
You must be signed in to change notification settings - Fork 0
/
report-alignments.py
56 lines (48 loc) · 1.75 KB
/
report-alignments.py
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
from __future__ import print_function
import sys
import os
def countAlignments(sam):
'''
Description: counts aligned sequences (unique accessions)
In: sam (filename)
Out: count (int)
'''
accs = dict()
try:
fh = open(sam)
except:
sys.exit("cannot open file: " + sam)
for line in fh:
line = line.rstrip()
if line.startswith("@"): # skip sam header
continue
line = line.split("\t")
if line[1] != "4": # check if alignment flag is NOT 4 (unmapped)
accs[line[0]] = accs.get(line[0], 0) + 1
count = len(accs)
return(count)
if __name__ == '__main__':
if len(sys.argv) < 2:
print("Usage:", sys.argv[0], "*.sam")
fhValign = open("report-ALIGNED-V.txt", "w")
fhJalign = open("report-ALIGNED-J.txt", "w")
fhValignAgain = open("report-ALIGNED-AGAIN-V.txt", "w")
fhJalignAgain = open("report-ALIGNED-AGAIN-J.txt", "w")
fhOther = open("report-ALIGNED-OTHER.txt", "w")
for sam in sys.argv[1:]:
sam_basename = os.path.basename(sam)
count = countAlignments(sam)
if "V_" in sam_basename and "-e-clean.sam" in sam_basename:
print(sam_basename, count, file=fhValign)
elif "J_" in sam_basename and "-e-clean.sam" in sam_basename:
print(sam_basename, count, file=fhJalign)
elif "V_" in sam_basename and "aligned-again" in sam_basename:
print(sam_basename, count, file=fhValignAgain)
elif "J_" in sam_basename and "aligned-again" in sam_basename:
print(sam_basename, count, file=fhJalignAgain)
else:
print(sam_basename, count, file=fhOther)
fhValign.close()
fhJalign.close()
fhValignAgain.close()
fhJalignAgain.close()