Allow vcf-expression-annotator to work on multiple samples or work on a VCF that has previously been annotated

[wanqiangdehuoguo]

For example,

$ head genetpm.tsv

  GeneID          N190533 T190533
ENSG00000000003  0.743   12.1  
ENSG00000000005  0.0232   0.115
ENSG00000000419 46.4     43.4  
ENSG00000000457  5.22     6.26 
ENSG00000000460  9.80     4.45 
ENSG00000000938  2.22    31.5  

190533.vep.vcf is a mulit sample vcf for somatic mutation:

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  N190533 T190533
chr1    1041944 .       C       G       .       PASS AS_FilterStatus=SITE;AS_SB_TABLE=209,110|41,22;DP=397;ECNT=1;GERMQ=93;MBQ=20,20;MFRL=229,220;MMQ=60,60;MPOS=37;NALOD=2.01;NLOD=30.01;POPAF=6.00;ROQ=93;TLOD=136.09;CSQ=G|splice_polypyrimidine_tract_variant&intron_variant|LOW|AGRN|ENSG00000188157|Transcript|ENST00000379370.7|protein_coding||6/35|ENST00000379370.7:c.1178-12C>G|||||||||1||HGNC|HGNC:329|1|||MAGR......VVVGRHPLHLLEDAVTKPELRPCPTP      GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/0:142,0:9.708e-03:142:48,0:47,0:100,0:90,52,0,0       0/1:177,63:0.26 2:240:61,21:67,23:134,47:119,58,41,22

The error is:

$ vcf-expression-annotator --sample-name N190533 --id-column GeneID --expression-column N190533 --output-vcf 190533.vep1.vcf 190533.vep.vcf genetpm.tsv custom gene
WARNING:root:69 of 1300 genes did not have an expression entry for their gene id.
$ vcf-expression-annotator --sample-name T190533 --id-column GeneID --expression-column T190533 --output-vcf 190533.vep2.vcf 190533.vep1.vcf genetpm.tsv custom gene
Traceback (most recent call last):
  File "/home/sym/.conda/envs/pVACtools/bin/vcf-expression-annotator", line 8, in <module>
    sys.exit(main())
  File "/home/sym/.conda/envs/pVACtools/lib/python3.10/site-packages/vatools/vcf_expression_annotator.py", line 191, in main
    (vcf_reader, is_multi_sample) = create_vcf_reader(args)
  File "/home/sym/.conda/envs/pVACtools/lib/python3.10/site-packages/vatools/vcf_expression_annotator.py", line 90, in create_vcf_reader
    raise Exception("ERROR: VCF {} is already gene expression annotated. GX format header already exists.".format(args.input_vcf))
Exception: ERROR: VCF 190533.vep1.vcf is already gene expression annotated. GX format header already exists.

And same error is

$ vcf-expression-annotator --sample-name N190533,T190533 --id-column GeneID --expression-column N190533,T190533 --output-vcf 190533.vep1.vcf 190533.vep.vcf genetpm.tsv custom gene
Traceback (most recent call last):
  File "/home/sym/.conda/envs/pVACtools/bin/vcf-expression-annotator", line 8, in <module>
    sys.exit(main())
  File "/home/sym/.conda/envs/pVACtools/lib/python3.10/site-packages/vatools/vcf_expression_annotator.py", line 191, in main
    (vcf_reader, is_multi_sample) = create_vcf_reader(args)
  File "/home/sym/.conda/envs/pVACtools/lib/python3.10/site-packages/vatools/vcf_expression_annotator.py", line 84, in create_vcf_reader
    raise Exception("ERROR: VCF {} does not contain a sample column for sample {}.".format(args.input_vcf, args.sample_name))
Exception: ERROR: VCF 190533.vep.vcf does not contain a sample column for sample N190533,T190533.

[wanqiangdehuoguo]

I found a compensation method, but it's pretty complicated. I wish there was an easier way.
The method is

vcf-expression-annotator --sample-name N190533 --id-column GeneID --expression-column N190533 --output-vcf 190533.vep1.vcf 190533.vep.vcf genetpm.tsv custom gene
vcf-expression-annotator --sample-name T190533 --id-column GeneID --expression-column T190533 --output-vcf 190533.vep2.vcf 190533.vep.vcf genetpm.tsv custom gene

bcftools view -s N190533 190533.vep2.vcf -o N190533.vcf
bcftools view -s T190533 190533.vep1.vcf -o T190533.vcf
bgzip N190533.vcf
bgzip T190533.vcf
tabix -p vcf N190533.vcf.gz
tabix -p vcf T190533.vcf.gz
bcftools merge N190533.vcf.gz T190533.vcf.gz -o 190533.vep.gx.vcf

[susannasiebert]

It looks like you VCF is maybe a matched tumor-normal VCF. On our end, we don't usually annotated the normal sample as the normal expression is not used in our downstream tools (e.g. pVACtools). However, I understand that there might be other use cases that require this information for both samples. Have you tried removing the GX header from the 190533.vep1.vcf VCF before running it a second time for T190533? Not sure if that would work but let me know. If that doesn't work, I'm happy to consider an update to the tool to support this use case.

[wanqiangdehuoguo]

Got it—my VCF is a matched tumor-normal file. For personal reasons, I’d like to append gene-expression (GX) annotations to it for both samples. I’ve already worked around the issue as described above, but it would be more convenient if vcf-expression-annotator could natively support multi-sample annotation. The script solution works：

#!/bin/bash  

# Define input parameters
inputvcf=$1  # Input VCF file
geneexpressionsv=$2  # Gene expression VCF file
samplename1=$3  # Sample 1 name
samplename2=$4  # Sample 2 name
id_column=$5  # ID column in the VCF file
outputvcf=$6  # Output VCF file

# Activate pVACtools environment
source /opt/local/anaconda3/etc/profile.d/conda.sh
conda activate pVACtools

# Run vcf-expression-annotator for sample 1 annotation
vcf-expression-annotator --sample-name $samplename1 --id-column $id_column --expression-column $samplename1 --output-vcf temp_vep1.vcf $inputvcf $geneexpressionsv custom gene

# Run vcf-expression-annotator for sample 2 annotation
vcf-expression-annotator --sample-name $samplename2 --id-column $id_column --expression-column $samplename2 --output-vcf temp_vep2.vcf $inputvcf $geneexpressionsv custom gene

# Use bcftools to extract data for sample 1 and sample 2
bcftools view -s $samplename1 temp_vep1.vcf -o ${samplename1}.vcf
bcftools view -s $samplename2 temp_vep2.vcf -o ${samplename2}.vcf

# Compress and index VCF files
bgzip ${samplename1}.vcf
bgzip ${samplename2}.vcf
tabix -p vcf ${samplename1}.vcf.gz
tabix -p vcf ${samplename2}.vcf.gz

# Merge the VCF files of both samples
bcftools merge ${samplename1}.vcf.gz ${samplename2}.vcf.gz -o $outputvcf

# Remove temporary files
rm temp_vep1.vcf temp_vep2.vcf
rm ${samplename1}.vcf.gz ${samplename2}.vcf.gz
rm ${samplename1}.vcf.gz.tbi ${samplename2}.vcf.gz.tbi

# Deactivate conda environment
conda deactivate

[susannasiebert]

Makes sense. I will keep this issue open to track development efforts for this feature.