vcf-expression-annotator updates

README.md

@@ -15,7 +15,7 @@ A tool that will add the data from bam-readcount files to the VCF sample column.
 
 **vcf-expression-annotator**
 
-A tool that will add the data from several expression tools’ output files to the VCF INFO column. Supported tools are StringTie, Kallisto, and Cufflinks. There also is a `custom` option to annotate with data from any tab-delimited file.
+A tool that will add the data from several expression tools’ output files to the VCF FORMAT column, on a per-sample basis (use `-s` to select the sample for multi-sample VCFs). Supported tools are StringTie, Kallisto, and Cufflinks. There also is a `custom` option to annotate with data from any tab-delimited file.
 
 **vcf-info-annotator**
 

docs/index.rst

@@ -15,7 +15,8 @@ annotate VCF files with data from other tools.
 
 **vcf-expression-annotator**
     A tool that will add the data from several expression tools' output files
-    to the VCF INFO column. Supported tools are StringTie, Kallisto,
+    to the VCF FORMAT column, on a per-sample basis (use ``-s`` to select the
+    sample for multi-sample VCFs). Supported tools are StringTie, Kallisto,
     and Cufflinks. There also is a ``custom`` option to annotate with data
     from any tab-delimited file.
 

docs/vcf_expression_annotator.rst

@@ -36,6 +36,13 @@ By default the output VCF will be written to a ``.tx.vcf`` or ``.gx.vcf`` file n
 your input VCF file. You can set a different output file name using the
 ``--output-vcf`` parameter.
 
+By default, running the VCF Expression Annotator against a VCF that already has a
+``GX``/``TX`` FORMAT header raises an error. Use ``--force`` to annotate such a VCF
+anyway - for example, to fill in a sample that wasn't previously annotated. Even with
+``--force``, an error is raised if doing so would overwrite an existing non-blank
+``GX``/``TX`` value for the target sample; use ``--overwrite`` (which requires
+``--force``) to allow replacing those existing values.
+
 Usage
 -----
 

tests/test_vcf_expression_annotator.py

@@ -101,6 +101,83 @@ def test_error_already_GX_annotated(self):
             vcf_expression_annotator.main(command)
         self.assertTrue('is already gene expression annotated. GX format header already exists.' in str(context.exception))
 
+    def test_error_overwrite_requires_force(self):
+        with self.assertRaises(SystemExit) as context:
+            vcf_expression_annotator.main([
+                os.path.join(self.test_data_dir, 'input.vcf'),
+                os.path.join(self.test_data_dir, 'genes.fpkm_tracking'),
+                'cufflinks',
+                'gene',
+                '--overwrite',
+            ])
+        self.assertEqual(context.exception.code, 2)
+
+    def test_error_force_without_overwrite_blocks_non_blank_value(self):
+        with self.assertRaises(Exception) as context:
+            command = [
+                os.path.join(self.test_data_dir, 'multiple_samples.gx.vcf'),
+                os.path.join(self.test_data_dir, 'genes.fpkm_tracking'),
+                'cufflinks',
+                'gene',
+                '-s', 'H_NJ-HCC1395-HCC1395',
+                '--force',
+            ]
+            vcf_expression_annotator.main(command)
+        self.assertTrue('already has a non-blank GX value' in str(context.exception))
+        self.assertTrue('Use --overwrite to replace existing values.' in str(context.exception))
+
+    def test_force_annotates_sample_with_blank_value(self):
+        temp_path = tempfile.TemporaryDirectory()
+        output_vcf = os.path.join(temp_path.name, 'output.vcf')
+        command = [
+            os.path.join(self.test_data_dir, 'multiple_samples.gx.vcf'),
+            os.path.join(self.test_data_dir, 'genes.fpkm_tracking'),
+            'cufflinks',
+            'gene',
+            '-s', 'H_NJ-HCC1395-HCC1396',
+            '--force',
+            '-o', output_vcf,
+        ]
+        vcf_expression_annotator.main(command)
+        with open(output_vcf) as f:
+            data_line = [l for l in f if not l.startswith('#')][0]
+        fields = data_line.rstrip('\n').split('\t')
+        gx_index = fields[8].split(':').index('GX')
+        self.assertEqual(fields[9].split(':')[gx_index], 'ENSG00000225255|0.0')
+        self.assertEqual(fields[10].split(':')[gx_index], 'ENSG00000225255|0.0')
+        temp_path.cleanup()
+
+    def test_force_overwrite_replaces_existing_value(self):
+        # multiple_samples.gx.vcf already has GX=ENSG00000225255|0.0 for this sample.
+        # Use an expression file with a different FPKM for that gene so the output can
+        # only match if the value was actually recomputed and overwritten, not just left as-is.
+        temp_path = tempfile.TemporaryDirectory()
+        updated_expression_file = os.path.join(temp_path.name, 'genes.fpkm_tracking')
+        with open(os.path.join(self.test_data_dir, 'genes.fpkm_tracking')) as infile, open(updated_expression_file, 'w') as outfile:
+            for line in infile:
+                fields = line.rstrip('\n').split('\t')
+                if fields[0] == 'ENSG00000225255':
+                    fields[9] = '42'
+                outfile.write('\t'.join(fields) + '\n')
+        output_vcf = os.path.join(temp_path.name, 'output.vcf')
+        command = [
+            os.path.join(self.test_data_dir, 'multiple_samples.gx.vcf'),
+            updated_expression_file,
+            'cufflinks',
+            'gene',
+            '-s', 'H_NJ-HCC1395-HCC1395',
+            '--force',
+            '--overwrite',
+            '-o', output_vcf,
+        ]
+        vcf_expression_annotator.main(command)
+        with open(output_vcf) as f:
+            data_line = [l for l in f if not l.startswith('#')][0]
+        fields = data_line.rstrip('\n').split('\t')
+        gx_index = fields[8].split(':').index('GX')
+        self.assertEqual(fields[9].split(':')[gx_index], 'ENSG00000225255|42.0')
+        temp_path.cleanup()
+
     def test_error_id_column_nonexistent_in_file(self):
         with self.assertRaises(Exception) as context:
             command = [

vatools/vcf_expression_annotator.py

@@ -52,6 +52,13 @@ def to_array(dictionary):
         array.append("{}|{}".format(key, value))
     return sorted(array)
 
+def is_blank_format_value(value):
+    if value is None:
+        return True
+    if isinstance(value, list):
+        return len(value) == 0 or all(v in (None, '.', '') for v in value)
+    return value in ('.', '')
+
 def parse_expression_file(args, vcf_reader, vcf_writer):
     if args.format == 'stringtie' and args.mode == 'transcript':
         df_all = read_gtf(args.expression_file, usecols=['reference_id', 'transcript_id', 'TPM', 'feature'])
@@ -85,28 +92,31 @@ def create_vcf_reader(args):
     if 'CSQ' not in vcf_reader.header.info_ids():
         vcf_reader.close()
         raise Exception("ERROR: VCF {} is not VEP-annotated. Please annotate the VCF with VEP before running this tool.".format(args.input_vcf))
-    if args.mode == 'gene' and 'GX' in vcf_reader.header.format_ids():
+    if args.mode == 'gene' and 'GX' in vcf_reader.header.format_ids() and not args.force:
         vcf_reader.close()
-        raise Exception("ERROR: VCF {} is already gene expression annotated. GX format header already exists.".format(args.input_vcf))
-    elif args.mode == 'transcript' and 'TX' in vcf_reader.header.format_ids():
+        raise Exception("ERROR: VCF {} is already gene expression annotated. GX format header already exists. Use --force to annotate anyway.".format(args.input_vcf))
+    elif args.mode == 'transcript' and 'TX' in vcf_reader.header.format_ids() and not args.force:
         vcf_reader.close()
-        raise Exception("ERROR: VCF {} is already transcript expression annotated. TX format header already exists.".format(args.input_vcf))
-    return vcf_reader, is_multi_sample
+        raise Exception("ERROR: VCF {} is already transcript expression annotated. TX format header already exists. Use --force to annotate anyway.".format(args.input_vcf))
+    sample_name = args.sample_name if is_multi_sample else vcf_reader.header.samples.names[0]
+    return vcf_reader, sample_name
 
 def create_vcf_writer(args, vcf_reader):
     (head, sep, tail) = args.input_vcf.rpartition('.vcf')
     new_header = vcf_reader.header.copy()
     if args.mode == 'gene':
-        new_header.add_format_line(OrderedDict([('ID', 'GX'), ('Number', '.'), ('Type', 'String'), ('Description', 'Gene Expressions')]))
+        if 'GX' not in new_header.format_ids():
+            new_header.add_format_line(OrderedDict([('ID', 'GX'), ('Number', '.'), ('Type', 'String'), ('Description', 'Gene Expressions')]))
         output_file = ('').join([head, '.gx.vcf', tail])
     elif args.mode == 'transcript':
-        new_header.add_format_line(OrderedDict([('ID', 'TX'), ('Number', '.'), ('Type', 'String'), ('Description', 'Transcript Expressions')]))
+        if 'TX' not in new_header.format_ids():
+            new_header.add_format_line(OrderedDict([('ID', 'TX'), ('Number', '.'), ('Type', 'String'), ('Description', 'Transcript Expressions')]))
         output_file = ('').join([head, '.tx.vcf', tail])
     if args.output_vcf:
         output_file = args.output_vcf
     return vcfpy.Writer.from_path(output_file, new_header)
 
-def add_expressions(entry, is_multi_sample, sample_name, df, items, tag, id_column, expression_column, ignore_ensembl_id_version, missing_expressions_count):
+def add_expressions(entry, sample_name, df, items, tag, id_column, expression_column, ignore_ensembl_id_version, missing_expressions_count, overwrite):
     subset = None
     if ignore_ensembl_id_version:
         items_without_version = [re.sub(r'\.[0-9]+$', '', item) for item in items]
@@ -117,18 +127,25 @@ def add_expressions(entry, is_multi_sample, sample_name, df, items, tag, id_colu
         subset = df[df[id_column].isin(items)]
     expressions = subset[[id_column, expression_column]].groupby(id_column).sum().to_dict()[expression_column]
     missing_expressions_count += (len(items) - len(expressions))
-    if is_multi_sample:
-        entry.FORMAT += [tag]
-        entry.call_for_sample[sample_name].data[tag] = to_array(expressions)
-    else:
-        entry.add_format(tag, to_array(expressions))
+    call = entry.call_for_sample[sample_name]
+    existing_value = call.data.get(tag)
+    if not is_blank_format_value(existing_value) and not overwrite:
+        raise Exception(
+            "ERROR: Sample {} already has a non-blank {} value ({}) at {}:{}. Use --overwrite to replace existing values.".format(
+                sample_name, tag, existing_value, entry.CHROM, entry.POS
+            )
+        )
+    if tag not in entry.FORMAT:
+        entry.FORMAT.append(tag)
+    call.data[tag] = to_array(expressions)
     return (entry, missing_expressions_count)
 
 def define_parser():
     parser = argparse.ArgumentParser(
         "vcf-expression-annotator",
         description = "A tool that will add the data from several expression tools' output files" +
-                      "to the VCF INFO column. Supported tools are StringTie, Kallisto, " +
+                      "to the VCF FORMAT column, on a per-sample basis (use -s to select the " +
+                      "sample for multi-sample VCFs). Supported tools are StringTie, Kallisto, " +
                       "and Cufflinks. There also is a ``custom`` option to annotate with data " +
                       "from any tab-delimited file."
     )
@@ -175,6 +192,20 @@ def define_parser():
         help='Assumes that the final period and number denotes the Ensembl ID version and ignores it (i.e. for "ENST00001234.3" - ignores the ".3").',
         action="store_true"
     )
+    parser.add_argument(
+        "--force",
+        action="store_true",
+        help="Allow annotating a VCF that already has a GX/TX FORMAT header (e.g. to fill in "
+            +"a sample that wasn't previously annotated). By default, running against an "
+            +"already-annotated VCF raises an error."
+    )
+    parser.add_argument(
+        "--overwrite",
+        action="store_true",
+        help="Requires --force. Allow overwriting existing non-blank GX/TX values for the "
+            +"target sample. Without this flag, --force will still raise an error if the "
+            +"sample already has a non-blank value for a variant."
+    )
 
     return parser
 
@@ -187,8 +218,10 @@ def main(args_input = sys.argv[1:]):
             raise Exception("--id-column is not set. This is required when using the `custom` format.")
         if args.expression_column is None:
             raise Exception("--expression-column is not set. This is required when using the `custom` format.")
+    if args.overwrite and not args.force:
+        parser.error("--overwrite requires --force")
 
-    (vcf_reader, is_multi_sample) = create_vcf_reader(args)
+    (vcf_reader, sample_name) = create_vcf_reader(args)
     format_pattern = re.compile('Format: (.*)')
     csq_format = format_pattern.search(vcf_reader.header.get_info_field_info('CSQ').description).group(1).split('|')
 
@@ -215,12 +248,12 @@ def main(args_input = sys.argv[1:]):
         if args.mode == 'gene':
             genes = list(genes)
             if len(genes) > 0:
-                (entry, missing_expressions_count) = add_expressions(entry, is_multi_sample, args.sample_name, df, genes, 'GX', id_column, expression_column, args.ignore_ensembl_id_version, missing_expressions_count)
+                (entry, missing_expressions_count) = add_expressions(entry, sample_name, df, genes, 'GX', id_column, expression_column, args.ignore_ensembl_id_version, missing_expressions_count, args.overwrite)
                 entry_count += len(genes)
         elif args.mode == 'transcript':
             transcript_ids = list(transcript_ids)
             if len(transcript_ids) > 0:
-                (entry, missing_expressions_count) = add_expressions(entry, is_multi_sample, args.sample_name, df, transcript_ids, 'TX', id_column, expression_column, args.ignore_ensembl_id_version, missing_expressions_count)
+                (entry, missing_expressions_count) = add_expressions(entry, sample_name, df, transcript_ids, 'TX', id_column, expression_column, args.ignore_ensembl_id_version, missing_expressions_count, args.overwrite)
                 entry_count += len(transcript_ids)
         vcf_writer.write_record(entry)