From f1092d8a9f09ff6a347a3933b40ed8e8beeb952d Mon Sep 17 00:00:00 2001 From: Chris Miller Date: Mon, 6 Jul 2026 17:26:34 -0500 Subject: [PATCH] vcf-expression-annotator updates --- README.md | 2 +- docs/index.rst | 3 +- docs/vcf_expression_annotator.rst | 7 +++ tests/test_vcf_expression_annotator.py | 77 ++++++++++++++++++++++++++ vatools/vcf_expression_annotator.py | 67 ++++++++++++++++------ 5 files changed, 137 insertions(+), 19 deletions(-) diff --git a/README.md b/README.md index 7d85854..85459f6 100644 --- a/README.md +++ b/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** diff --git a/docs/index.rst b/docs/index.rst index 73669b7..7754b35 100644 --- a/docs/index.rst +++ b/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. diff --git a/docs/vcf_expression_annotator.rst b/docs/vcf_expression_annotator.rst index 1d28846..a0ceda6 100644 --- a/docs/vcf_expression_annotator.rst +++ b/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 ----- diff --git a/tests/test_vcf_expression_annotator.py b/tests/test_vcf_expression_annotator.py index 91a8b1b..2b2e7ea 100644 --- a/tests/test_vcf_expression_annotator.py +++ b/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 = [ diff --git a/vatools/vcf_expression_annotator.py b/vatools/vcf_expression_annotator.py index 1818feb..63bba8d 100644 --- a/vatools/vcf_expression_annotator.py +++ b/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)