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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
|
#!/usr/bin/python3
# Copyright (C) 2011 Ion Torrent Systems, Inc. All Rights Reserved
import sys
import re
import os
import subprocess
import time
import json
from optparse import OptionParser
def printtime(message, *args):
if args:
message = message % args
print("[ " + time.strftime('%a %Y-%m-%d %X %Z') + " ] " + message)
sys.stdout.flush()
sys.stderr.flush()
def RunCommand(command,description):
printtime(' ')
printtime('Task : ' + description)
printtime('Command : ' + command)
printtime(' ')
stat = subprocess.call(command,shell=True)
if stat != 0:
printtime('ERROR: command failed with status %d' % stat)
sys.exit(1)
def main():
parser = OptionParser()
parser.add_option('-b', '--region-bed', help='BED file specifing regions over which variant calls will be limited or filtered to', dest='bedfile')
parser.add_option('-s', '--hotspot-vcf', help='VCF file specifying exact hotspot positions. TVC will force an evaluation for the alleles specified at the hotspot positions in this VCF file. For details, please visit Hotspot Variants (optional)', dest='hotspot_vcf')
parser.add_option('-i', '--input-bam', help='Input BAM file(s) containing aligned reads. Multiple file names must be concatenated with commas (required)', dest='bamfile')
parser.add_option('-n', '--normal-bam', help='BAM file(s) containing aligned reads for normal (reference) sample (IR only, optional)', dest='normbamfile')
parser.add_option('-g', '--sample-name', help='Sample Name for Test Sample (IR only, optional)', dest='testsamplename')
parser.add_option('-r', '--reference-fasta', help='FASTA file containing reference genome (required)', dest='reference')
parser.add_option('-o', '--output-dir', help='Output directory (default: current)', dest='outdir', default='.')
parser.add_option('-p', '--parameters-file', help='JSON file containing variant calling parameters. This file can be obtained from https://ampliseq.com for your panel. If not provided, default params will be used. for more information about parameters, please visit TVC 4.x Parameters Description (optional, recommended)', dest='paramfile')
parser.add_option('-m', '--error-motifs', help='System dependent motifs file helps improve variant calling accuracy. For Hi-Q chemistry use $TVC_ROOT_DIR/share/TVC/sse/ampliseqexome_germline_p1_hiq_motifset.txt else use $TVC_ROOT_DIR/share/TVC/sse/motifset.txt (optional)', dest='errormotifsfile')
parser.add_option('-N', '--num-threads', help='Set TVC number of threads (default: 12)', dest='numthreads',default='12')
parser.add_option('-B', '--bin-dir', help='Directory path to location of variant caller programs. Defaults to the directory this script is located', dest='bindir')
parser.add_option('-t', '--tvc-root-dir', help='Directory path to TVC root directory', dest='tvcrootdir')
parser.add_option('-G', '--generate-gvcf', help='Request generation of gvcf file in addition to vcf (on/off, default off)', dest='generate_gvcf', default='off')
parser.add_option( '--primer-trim-bed', help='Perform primer trimming using provided unmerged BED file. (optional, recommended for ampliseq)', dest='ptrim_bed')
parser.add_option( '--postprocessed-bam',help='If provided, a primer trimmed BAM file will be produced for IGV viewing. This file does not contain the flow space data and should not be used as an input to TVC. Use of this option may increase TVC run time (optional, not recommended)', dest='postprocessed_bam')
(options, args) = parser.parse_args()
#parser.add_option('-B', '--build-dir', help='Directory path to location of IR Build. Defaults to the directory this script is located', dest='builddir', default=os.path.realpath(__file__))
#parser.add_option('-m', '--motifs-file', help='File containing motif set', dest='motifsfile')
#parser.add_option('-s', '--hotspot-file', help='VCF or BED file specifying exact hotspot positions', dest='hotspot_vcf')
#parser.add_option('-c', '--contig-name', help='Name of the Contig over which variant calls will be limited', dest='contigname')
# Sort out directories and locations of the binaries
if not os.path.exists(options.outdir):
os.makedirs(options.outdir)
if options.tvcrootdir:
bin_dir = os.path.join(options.tvcrootdir, 'bin')
elif options.bindir:
bin_dir = options.bindir
else:
bin_dir = os.path.dirname(os.path.realpath(__file__))
bin_dir = os.path.normpath(bin_dir)
path_to_tvc = os.path.join(bin_dir,'tvc')
if not os.path.exists(path_to_tvc):
path_to_tvc = 'tvc'
path_to_tvcassembly = os.path.join(bin_dir,'tvcassembly')
if not os.path.exists(path_to_tvcassembly):
path_to_tvcassembly = 'tvcassembly'
path_to_tvcutils = os.path.join(bin_dir,'tvcutils')
if not os.path.exists(path_to_tvcutils):
path_to_tvcutils = 'tvcutils'
# Verify that all pre-conditions are met
if not options.bamfile or not options.reference:
parser.print_help()
exit(1)
if not os.path.isdir(options.outdir):
printtime('ERROR: No output directory found at: ' + options.outdir)
sys.exit(1)
if not os.path.exists(options.reference):
printtime('ERROR: No reference file found at: ' + options.reference)
sys.exit(1)
if not os.path.exists(options.reference+'.fai'):
printtime('ERROR: No reference index file found at: ' + options.reference + '.fai')
sys.exit(1)
for bam_filename in options.bamfile.split(','):
if not os.path.exists(bam_filename):
printtime('ERROR: No bam file found at: ' + bam_filename)
sys.exit(1)
if not os.path.exists(bam_filename+'.bai'):
printtime('ERROR: No bam index file found at: ' + bam_filename + '.bai')
sys.exit(1)
if options.hotspot_vcf:
if not os.path.exists(options.hotspot_vcf):
printtime('ERROR: No hotspots vcf file found at: ' + options.hotspot_vcf)
sys.exit(1)
if options.ptrim_bed:
if not os.path.exists(options.ptrim_bed):
printtime('ERROR: No primer trim bed file found at: ' + options.ptrim_bed)
sys.exit(1)
if options.bedfile:
if not os.path.exists(options.bedfile):
printtime('ERROR: No target regions bed file found at: ' + options.bedfile)
sys.exit(1)
# write effective bed file
if options.ptrim_bed and not options.bedfile:
options.bedfile = options.ptrim_bed
if options.ptrim_bed:
tvcutils_command = path_to_tvcutils + " validate_bed"
tvcutils_command += ' --reference "%s"' % options.reference
tvcutils_command += ' --target-regions-bed "%s"' % options.ptrim_bed
tvcutils_command += ' --effective-bed "%s"' % os.path.join(options.outdir,'effective_regions.bed')
RunCommand(tvcutils_command,'Write effective bed')
parameters = {}
if options.paramfile:
try:
json_file = open(options.paramfile, 'r')
parameters = json.load(json_file)
json_file.close()
if 'pluginconfig' in parameters:
parameters = parameters['pluginconfig']
except:
printtime('ERROR: No parameter file found at: ' + options.paramfile)
sys.exit(1)
# TVC
printtime('Calling small INDELs and SNPs using tvc ...')
meta_tvc_args = parameters.get('meta',{}).get('tvcargs','tvc')
if meta_tvc_args == 'tvc':
tvc_command = path_to_tvc
else:
tvc_command = meta_tvc_args
tvc_command += ' --output-dir %s' % options.outdir
tvc_command += ' --reference %s' % options.reference
if options.numthreads:
tvc_command += ' --num-threads %s' % options.numthreads
if options.ptrim_bed:
tvc_command += ' --target-file %s' % options.ptrim_bed
tvc_command += ' --trim-ampliseq-primers on'
elif options.bedfile:
tvc_command += ' --target-file %s' % options.bedfile
if options.postprocessed_bam:
#postprocessed_bam_tmp = options.postprocessed_bam + '.tmp.bam'
#tvc_command += ' --postprocessed-bam %s' % postprocessed_bam_tmp
tvc_command += ' --postprocessed-bam %s' % options.postprocessed_bam
if options.paramfile:
tvc_command += ' --parameters-file %s' % options.paramfile
if options.errormotifsfile:
tvc_command += ' --error-motifs %s' % options.errormotifsfile
if options.normbamfile:
tvc_command += ' --input-bam %s,%s' % ((options.bamfile, options.normbamfile))
tvc_command += ' --sample-name "%s"' % (options.testsamplename)
else:
tvc_command += ' --input-bam %s' % (options.bamfile)
tvc_command += ' --output-vcf small_variants.vcf'
if options.hotspot_vcf:
tvc_command += ' --input-vcf %s' % options.hotspot_vcf
if parameters and parameters['torrent_variant_caller'].get('process_input_positions_only','0') == '1' and parameters['freebayes'].get('use_input_allele_only','0') == '1':
tvc_command += ' --process-input-positions-only on'
tvc_command += ' --use-input-allele-only on'
RunCommand(tvc_command, 'Call Hotspots Only')
else:
RunCommand(tvc_command,'Call small indels and SNPs')
long_indel_command = path_to_tvcassembly
long_indel_command += ' --reference %s' % options.reference
if options.normbamfile:
long_indel_command += ' --input-bam %s,%s' % ((options.bamfile, options.normbamfile))
long_indel_command += ' --sample-name "%s"' % (options.testsamplename)
else:
long_indel_command += ' --input-bam %s' % (options.bamfile)
if options.bedfile:
long_indel_command += ' --target-file %s' % options.bedfile
long_indel_command += ' --output-vcf %s/indel_assembly.vcf' % options.outdir
if options.paramfile:
long_indel_command += ' --parameters-file %s' % options.paramfile
RunCommand(long_indel_command, "Assemble long indels")
if options.postprocessed_bam:
#bamsort_command = 'samtools sort -m 2G -l1 -@6 %s %s' % (postprocessed_bam_tmp, options.postprocessed_bam[:-4])
#RunCommand(bamsort_command,'Sort postprocessed bam')
bamindex_command = 'samtools index %s' % options.postprocessed_bam
RunCommand(bamindex_command,'Index postprocessed bam')
#RunCommand('rm -f ' + postprocessed_bam_tmp, 'Remove unsorted postprocessed bam')
if options.generate_gvcf == "on":
unify_command = 'samtools depth'
if options.bedfile != None:
unify_command += ' -b ' + options.bedfile
if options.postprocessed_bam:
unify_command += ' ' + options.postprocessed_bam
else:
unify_command += ' ' + ' '.join(options.bamfile.split(','))
unify_command += ' | '
else:
unify_command = ''
unify_command += path_to_tvcutils + ' unify_vcf'
unify_command += ' --novel-tvc-vcf %s/small_variants.vcf' % options.outdir
if os.path.exists("%s/indel_assembly.vcf" % options.outdir):
unify_command += ' --novel-assembly-vcf %s/indel_assembly.vcf' % options.outdir
if options.hotspot_vcf:
unify_command += ' --hotspot-annotation-vcf "%s"' % options.hotspot_vcf
unify_command += ' --output-vcf %s/TSVC_variants.vcf.gz' % options.outdir
unify_command += ' --reference-fasta %s' % options.reference
if os.path.exists(options.outdir + '/tvc_metrics.json'):
unify_command += ' --tvc-metrics %s/tvc_metrics.json' % options.outdir
if options.bedfile:
unify_command += ' --target-file "%s"' % options.bedfile
if options.generate_gvcf == "on":
unify_command += ' --input-depth stdin'
if parameters and 'gen_min_coverage' in parameters.get('freebayes', {}):
unify_command +=' --min-depth ' + str(parameters['freebayes']['gen_min_coverage'])
RunCommand(unify_command, 'Unify variants and annotations from all sources (tvc,IndelAssembly,hotpots)')
# Generate uncompressed vcf file
RunCommand('gzip -dcf "%s/TSVC_variants.vcf.gz" > "%s/TSVC_variants.vcf"' % (options.outdir,options.outdir), 'Generate uncompressed vcf')
if options.generate_gvcf == "on" and os.path.exists(options.outdir+'/TSVC_variants.genome.vcf.gz'):
RunCommand('gzip -dcf "%s/TSVC_variants.genome.vcf.gz" > "%s/TSVC_variants.genome.vcf"' % (options.outdir,options.outdir), 'Generate uncompressed gvcf')
if __name__ == '__main__':
main()
|
