crunch_scripts/GATK2-VariantFiltration

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