Locally annotating SNP IDs and Gene names of called variants

Question

[entropy]

I have GWAS results after variant calling. The VCF file only had CHR (1:22) and POS (12345678 etc) information but the ID column has all ".", namely no rsIDs in it. After GWAS analysis I have a list of variants with CHR and POS information along with their statistical results from Plink2 as a tab delimited text file.

I want to annotate these variants in command line using a tool like VEP if possible. I want to get rsIDs and corresponding gene symbols (GRCh38) and other important annotations. VEP documentation fro command line is a bit not clear to me for my multiple annotation tasks.

Could you please provide me an example command along with example input file names for running VEP for this goal? If VEP is not sufficient for the task alone, other tools like snpEff is also ok.

My GWAS result text file has a header and data like below:

CHR    POS       ID     P         BETA 
4      123456     .     0.005     23.4
... and so on

[entropy]

Hi Emily I see you suggest VEP for similar questions. It would be great if you could also give an example commadn for my query? Thanks.

Answer 1

[barslmn]

VEP has a great web interface You can try using that.

[entropy]

Yes, but I want to make VEP as part of my pipeline and thus would like to use the command line version to automate the process. The documentation was a bit confusing for me and could not see an example command that shows how make annotations for both rsID and Gene symbols at he same time given that there is only chr id and coordinate information.

[barslmn]

Have you installed the VEP? What was the command you tried?

[entropy]

Yes I downloaded the VEP cache and you can see the details of it below.

I have also craffetd a command that mostly working but it does not return rsIDs. Instead it returned COSMIC IDs. I am currently looking for a way to get rsIDs instead of COSMIC. Below is the command that mostly worked and one example output of VEP:

vep -i vep_ex_input.txt --cache --dir_cache /home/.vep --offline --fork 5 --sift b --canonical  --symbol  --check_existing --assembly GRCh38 --tab --polyphen b --fields "#Uploaded_variation,SYMBOL,Existing_variation,Consequence,IMPACT,SIFT,PolyPhen,CANONICAL,CLIN_SIG,PHENO,Location,Allele,Gene,Feature,Feature_type,GENE_PHENO,BIOTYPE,VARIANT_CLASS,ENSP,NEAREST,STRAND,SYMBOL_SOURCE,UPLOADED_ALLELE,REF_ALLELE,USED_REF,GIVEN_REF,SOMATIC,MAX_AF,AF,EUR_AF,AMR_AF,PUBMED,CELL_TYPE" -o outputs/annotated_output.txt

An an example output of VEP:

-   ZYZ2        COSV10011132    3_prime_UTR_variant,NMD_transcript_variant  MODIFIER    -   -   -   -   1   8:123456789     A   ENSG00000001233 ENST00000123457 Transcript  -   -   -   -   -   -1  HGNC    -   -   -   -   1   -   -   -   -   -   -

And here is the VEP information:

$ vep --show_cache_info --dir_cache /home/.vep
Smartmatch is experimental at /app/ensembl-vep/modules/Bio/EnsEMBL/VEP/AnnotationSource/File.pm line 472.
gnomADg v3.1.2
gnomADe r2.1.1
ClinVar 202301
polyphen        2.2.3
sift    6.2.1
1000genomes     phase3
assembly        GRCh38.p14
COSMIC  97
regbuild        1.0
genebuild       2014-07
dbSNP   154
HGMD-PUBLIC     20204
gencode GENCODE 44

[barslmn]

8:123456789

In this case there is no rsid for that variant. https://gnomad.broadinstitute.org/region/8-123456769-123456809?dataset=gnomad_r4

Could you try a variant we know to have an rsid like this: 17 43044391 G A

https://www.ncbi.nlm.nih.gov/snp/rs12516

[entropy]

I have run for all the genes, which are 15660 genes in my data and I could not see a single rsID. The "Existing_variation" columns has all the variant IDs. It has either "-" or "COSV12345" like content.

[barslmn]

Only other thing comes to my mind is to the with input format.

Could you format your input file into one of these formats: https://www.ensembl.org/info/docs/tools/vep/vep_formats.html

[entropy]

I think I already did so. Here is the head example of my input data format:

1 12345678  12345678    C/T      +
4  32345678 32345678    G/A      +
2 52345678  52345678    T/C      +

[barslmn]

Could you try a variant we know to have an rsid like this: 17 43044391 G A

https://www.ncbi.nlm.nih.gov/snp/rs12516

Could you still try running this to see if rsid is returned?

It would be like this in VEP format:

17 43044391 43044391 G/A +