3177: store input pdh value even when input dataclass is File. Update two tests to...
[arvados.git] / crunch_scripts / GATK2-VariantFiltration
1 #!/usr/bin/env python
2
3 import arvados
4 import os
5 import re
6
7 arvados.job_setup.one_task_per_input_file(if_sequence=0, and_end_task=True)
8
9 this_job = arvados.current_job()
10 this_task = arvados.current_task()
11 gatk_path = arvados.util.tarball_extract(
12     tarball = this_job['script_parameters']['gatk_binary_tarball'],
13     path = 'gatk')
14 bundle_path = arvados.util.collection_extract(
15     collection = this_job['script_parameters']['gatk_bundle'],
16     path = 'gatk-bundle',
17     files = ['human_g1k_v37.dict', 'human_g1k_v37.fasta', 'human_g1k_v37.fasta.fai'])
18 this_task_input = this_task['parameters']['input']
19
20 input_file = list(arvados.CollectionReader(this_task_input).all_files())[0]
21
22 # choose vcf temporary file names
23 vcf_in = os.path.join(arvados.current_task().tmpdir,
24                       os.path.basename(input_file.name()))
25 vcf_out = re.sub('(.*)\\.vcf', '\\1-filtered.vcf', vcf_in)
26
27 # fetch the unfiltered data
28 vcf_in_file = open(vcf_in, 'w')
29 for buf in input_file.readall():
30     vcf_in_file.write(buf)
31 vcf_in_file.close()
32
33 stdoutdata, stderrdata = arvados.util.run_command(
34     ['java', '-Xmx1g',
35      '-jar', os.path.join(gatk_path,'GenomeAnalysisTK.jar'),
36      '-T', 'VariantFiltration', '--variant', vcf_in,
37      '--out', vcf_out,
38      '--filterExpression', 'QD < 2.0',
39      '--filterName', 'GATK_QD',
40      '--filterExpression', 'MQ < 40.0',
41      '--filterName', 'GATK_MQ',
42      '--filterExpression', 'FS > 60.0',
43      '--filterName', 'GATK_FS',
44      '--filterExpression', 'MQRankSum < -12.5',
45      '--filterName', 'GATK_MQRankSum',
46      '--filterExpression', 'ReadPosRankSum < -8.0',
47      '--filterName', 'GATK_ReadPosRankSum',
48      '-R', os.path.join(bundle_path, 'human_g1k_v37.fasta')],
49     cwd=arvados.current_task().tmpdir)
50
51 # store the filtered data
52 with open(vcf_out, 'rb') as f:
53     out = arvados.CollectionWriter()
54     while True:
55         buf = f.read()
56         if len(buf) == 0:
57             break
58         out.write(buf)
59 out.set_current_file_name(os.path.basename(vcf_out))
60
61 this_task.set_output(out.finish())