/
cgCaller.py
executable file
·164 lines (144 loc) · 6.94 KB
/
cgCaller.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
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
161
162
163
164
#!/usr/bin/env python3
import sys
import argparse
import gzip
import bz2
import re
import string
import utils
import io
argparser = argparse.ArgumentParser()
argparser.add_argument("chr", help="Chromosome in the masterVar file")
argparser.add_argument("sample_id", help="Sample ID to put into the VCF file")
argparser.add_argument("out_mask", help="mask-file to write to, will be gzipped, so better add .gz to the end.")
argparser.add_argument("input", help="Complete Genomics masterVarBeta file (uncompressed or compressed with gzip or bzip2)")
argparser.add_argument("--max_pos", type=int, default=0)
argparser.add_argument("--mastervar_version", default="2.0", choices=['2.0','2.4'], help="version of masterVar file")
argparser.add_argument("--legend_file", help="Impute2 reference panel legend file, can be gzipped or not")
args = argparser.parse_args()
if args.mastervar_version == "2.0":
qpass = "VQHIGH"
elif args.mastervar_version == "2.4":
qpass = "PASS"
mask_generator = utils.MaskGenerator(args.out_mask, args.chr)
sites_parser = None
if args.legend_file is not None:
sites_parser = utils.LegendParser(args.legend_file)
print("##fileformat=VCFv4.1")
print('##FORMAT=<ID=GT,Number=1,Type=String,Description="Phased Genotype">')
print("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t{}".format(args.sample_id))
input_file = None
if args.input[-3:] == ".gz":
input_file = io.TextIOWrapper(gzip.open(args.input, "r"))
elif args.input[-4:] == ".bz2":
input_file = io.TextIOWrapper(bz2.open(args.input, "r"))
else:
input_file = open(args.input, "r")
line_count = 0
chromosome_read = False
for line in input_file:
if line[0] == '#' or line[0] == '>' or line == "\n":
continue
line_count += 1
fields = line.strip().split()
chrom = fields[2]
begin = int(fields[3])
end = int(fields[4])
if line_count % 100000 == 0:
sys.stderr.write("processing chromosome {}, position {}\n".format(chrom, begin))
if chrom != args.chr:
if chromosome_read:
break
else:
continue
if args.max_pos > 0 and end > args.max_pos:
break
chromosome_read = True
zygosity = fields[5]
var_type = fields[6]
if var_type == "ref" and zygosity == "hom":
for i in range(begin + 1, end + 1):
mask_generator.addCalledPosition(i)
if sites_parser is not None:
while not sites_parser.end and sites_parser.pos < i:
sites_parser.tick()
if sites_parser.pos == i:
print("{chrom}\t{pos}\t.\t{ref_a}\t{alt_a}\t.\tPASS\t.\tGT\t0/0".format(chrom=args.chr, pos=i,
ref_a=sites_parser.ref_a,
alt_a=sites_parser.alt_a))
if var_type == "snp":
if zygosity in ["hom", "het-ref", "het-alt"]:
if args.mastervar_version == "2.0" and (end - begin) != 1:
assert (end - begin) == 1
elif args.mastervar_version == "2.4" and (end - begin) != 1:
continue
allele_ref = fields[7]
if sites_parser is not None:
while not sites_parser.end and sites_parser.pos < begin + 1:
sites_parser.tick()
if sites_parser.pos == begin + 1:
assert allele_ref == sites_parser.ref_a
allele_1 = fields[8]
allele_2 = fields[9]
allele1_qual = fields[14]
allele2_qual = fields[15]
if allele1_qual == qpass and allele2_qual == qpass:
mask_generator.addCalledPosition(begin + 1)
allele_indices = []
alt_alleles = []
if allele_1 == allele_ref:
allele_indices.append(0)
else:
alt_alleles.append(allele_1)
allele_indices.append(1)
if allele_2 == allele_ref:
allele_indices.append(0)
elif allele_2 == allele_1:
allele_indices.append(1)
else:
alt_alleles.append(allele_2)
allele_indices.append(len(alt_alleles))
print("{chrom}\t{pos}\t.\t{ref_a}\t{alt_a}\t.\tPASS\t.\tGT\t{gen1}/{gen2}".format(chrom=args.chr,
pos=begin+1, ref_a=allele_ref, alt_a=",".join(alt_alleles), gen1=allele_indices[0], gen2=allele_indices[1]))
# if var_type == "sub":
# if zygosity in ["hom", "het-ref", "het-alt"]:
# allele_ref = fields[7]
# allele_1 = fields[8]
# allele_2 = fields[9]
# allele1_qual = fields[14]
# allele2_qual = fields[15]
# if allele1_qual == "VQHIGH" and allele2_qual == "VQHIGH" and len(allele_ref)==len(allele_1) and len(allele_ref)==len(allele_2):
# for i in range(0,len(allele_1)):
# mask_generator.addCalledPosition(begin + 1)
# allele_indices = []
# alt_alleles = []
# if allele_1[i] != allele_ref[i]:
# alt_alleles.append(allele_1[i])
# allele_indices.append(1)
# elif allele_1[i] == allele_ref[i]:
# allele_indices.append(0)
# if allele_2[i] == allele_ref[i]:
# allele_indices.append(0)
# elif allele_2[i] == allele_1[i]:
# allele_indices.append(1)
# else:
# if len(alt_alleles)==0:
# alt_alleles.append(allele_2[i])
# allele_indices.append(1)
# else:
# alt_alleles.append(allele_2[i])
# allele_indices.append(2)
# if allele_indices[0]==0 and allele_indices[1]==0:
# if sites_parser is not None:
# while not sites_parser.end and sites_parser.pos < begin+i+1:
# sites_parser.tick()
# if sites_parser.pos == begin+i+1:
# assert allele_ref[i] == sites_parser.ref_a
# print("{chrom}\t{pos}\t.\t{ref_a}\t{alt_a}\t.\tPASS\t.\tGT\t0/0".format(chrom=args.chr, pos=begin+i+1,ref_a=sites_parser.ref_a,alt_a=sites_parser.alt_a))
# else:
# if sites_parser is not None:
# while not sites_parser.end and sites_parser.pos < begin+i+1:
# sites_parser.tick()
# if sites_parser.pos == begin+i+1:
# assert allele_ref[i] == sites_parser.ref_a
# print("{chrom}\t{pos}\t.\t{ref_a}\t{alt_a}\t.\tPASS\t.\tGT\t{gen1}/{gen2}".format(chrom=args.chr, pos=begin+i+1, ref_a=allele_ref[i], alt_a=",".join(alt_alleles), gen1=allele_indices[0], gen2=allele_indices[1]))