-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathvcf_exome_metrics.py
115 lines (90 loc) · 4.03 KB
/
vcf_exome_metrics.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
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from subprocess import call
import os
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("-i", "--input", help="GVCF.GZ file (can be the location on S3)")
parser.add_argument("-t", "--target", help="Target File")
args = parser.parse_args()
vcf_file = args.input
target_file = args.target
print(vcf_file, target_file)
base=os.path.basename(vcf_file)
# base_name = os.path.splitext(base)[0]
base_name = base.split('.')[0]
print(base_name)
#create one folder per sample
output_folder = "/home/ubuntu/projects/output/reports/vcf/%s" % (base_name)
if not os.path.exists(output_folder):
os.makedirs(output_folder)
output_base = "%s/%s" % (output_folder, base_name)
memory_use = "15g"
gvcftools_path = "/home/ubuntu/projects/programs/gvcftools-0.16/bin"
vcftools_path = "/home/ubuntu/projects/programs/vcftools/vcftools-0.1.14/src"
bcftools_path = "/home/ubuntu/projects/programs/bcftools/bcftools-1.3.1"
snpeff_path = "/home/ubuntu/projects/programs/snpeff/snpEff"
#/home/ubuntu/projects/programs/vcftools/vcftools-0.1.14/src/cpp/
#/home/ubuntu/projects/programs/vcftools/vcftools-0.1.14/src/perl/
#extract vcf from gvcf
print('extract vcf from gvcf')
#gzip -dc ../../input/WGC081270U.g.vcf.gz | ../../programs/gvcftools-0.16/bin/extract_variants | bgzip -c > WGC081270U.vcf.gz
command = """gzip -dc %s | %s/extract_variants | bgzip -c > %s.vcf.gz
""" % (vcf_file, gvcftools_path, output_base)
output = call(command, shell=True)
print(output)
#index with tabix
command = """tabix -p vcf %s.vcf.gz""" % (output_base)
output = call(command, shell=True)
print(output)
print('vcftools stats')
#vcftools metrics
command = """%s/perl/vcf-stats %s.vcf.gz > %s.vcftools.stats.txt
""" % (vcftools_path, output_base, output_base)
output = call(command, shell=True)
print(output)
#*coverage
#bcftools metrics
print('bcftools stats')
command = "%s/bcftools stats %s.vcf.gz > %s.bcftools.stats.txt" % (bcftools_path, output_base, output_base)
output = call(command, shell=True)
print(output)
#snpeff
#extract vcf
print('extract vcf')
command = "bgzip -d -c %s.vcf.gz > %s.vcf" % (output_base, output_base)
output = call(command, shell=True)
print(output)
print('snpeff')
command = "java -Xmx5g -jar %s/snpEff.jar eff -stats %s.snpeff.full.html -i vcf GRCh37.75 %s.vcf > %s.snpeff.full.vcf" % (snpeff_path, output_base, output_base, output_base)
output = call(command, shell=True)
print(output)
print('filtering')
#filtering VCF
command = "%s/bcftools filter -T %s %s.vcf > %s.filtered.exons.vcf" % (bcftools_path, target_file, output_base, output_base)
output = call(command, shell=True)
print(output)
#filtered.exons.q50.dp50.vcf
command = "%s/bcftools filter -T %s -i'QUAL>50 && FMT/DP>50' %s.vcf > %s.filtered.exons.q50.dp50.vcf" % (bcftools_path, target_file, output_base, output_base)
output = call(command, shell=True)
print(output)
#WGC081270U.filtered.exons.q100.dp100.vcf
command = "%s/bcftools filter -T %s -i'QUAL>100 && FMT/DP>100' %s.vcf > %s.filtered.exons.q100.dp100.vcf" % (bcftools_path, target_file, output_base, output_base)
output = call(command, shell=True)
print(output)
print('snpeff')
command = "java -Xmx5g -jar %s/snpEff.jar eff -stats %s.snpeff.exons.html -i vcf GRCh37.75 %s.filtered.exons.vcf > %s.snpeff.exons.vcf" % (snpeff_path, output_base, output_base, output_base)
output = call(command, shell=True)
print(output)
print('snpeff')
command = "java -Xmx5g -jar %s/snpEff.jar eff -stats %s.snpeff.exons.q100.dp100.html -i vcf GRCh37.75 %s.filtered.exons.q100.dp100.vcf > %s.snpeff.exons.q100.dp100.vcf" % (snpeff_path, output_base, output_base, output_base)
output = call(command, shell=True)
print(output)
print('bcftools stats')
command = "%s/bcftools stats %s.filtered.exons.q100.dp100.vcf > %s.bcftools.stats.exons.q100.dp100.txt" % (bcftools_path, output_base, output_base)
output = call(command, shell=True)
print(output)
print('plot vcf-stats')
command = "%s/plot-vcfstats %s.bcftools.stats.exons.q100.dp100.txt -p %s" % (bcftools_path, output_base, output_base)
output = call(command, shell=True)
print(output)