AssemblyAnalysis/CompareToStrandSeq.py at master · ChaissonLab/AssemblyAnalysis · GitHub

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
#!/usr/bin/env python
import sys

title=["correct-het", "correct-hom", "missed-het", "missed-hom", "fp-het", "fp-hom", "het-to-hom-alt", "hom-alt-to-het"]
correct_het=0
correct_hom=1
missed_het=2
missed_hom=3
fp_het=4
fp_hom=5
het_to_hom_alt=6
hom_alt_to_het=7
n=8
tally = [[0]*n, [0]*n]


for line in sys.stdin:
    if line[0] == "#":
        continue
    vals=line.split()
    gts=[vals[-3].split(":")[0], vals[-2].split(":")[0]]
    sts=vals[-1].split(":")[0].replace("|","/")
#    print((gts,sts))
    if sts == "1/0":
        sts = "0/1"

    for i in [0,1]:
        if sts == "0/0":
            if gts[i] == "0/1":
                tally[i][fp_het]+=1
            elif gts[i] == "1/1":
                tally[i][fp_hom]+=1
        if sts == "0/1":
            if gts[i] == "0/1":
                tally[i][correct_het]+=1
            elif gts[i] == "1/1":
                tally[i][het_to_hom_alt]+=1
            elif gts[i] == "0/0":
                tally[i][missed_het]+=1
        if sts == "1/1":
            if gts[i] == "0/0":
                tally[i][missed_hom]+=1
            if gts[i] == "0/1":
                tally[i][hom_alt_to_het] +=1
            if gts[i] == "1/1":
                tally[i][correct_hom] +=1

print("\t".join(title))
print("\t".join([str(i) for  i in tally[0]]))
print("\t".join([str(i) for i in tally[1]]))