I have RNA-seq datasets, and I applied GATK variant calling procedure for RNA-seq (https://software.broadinstitute.org/gatk/documentation/article.php?id=3891). I have not found anything wrong with it yet. However, when I took the .vcf file output from GATK and used snpEff with GRCh 37.75 to assess the impact of variants, it seems either the genomic coordinates were not handled correctly or the variant impact assessment mechanism was off. For instance, one variant output from GATK was documented like the following:

    chr19	7975586	.	C	T	275.78	.	AC=2;AF=1.00;AN=2;DP=9;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=30.64;SOR=4.615	GT:AD:DP:GQ:PL	1/1:0,9:9:26:304,26,0

This denotes the variant as T/T at chr19:7975586, which corresponds to codon GAT (D or Asp) in an exon of MAP2K7. The reference allele is C/C, which forms the codon GAC coding for Asp as well. In other words, this variant does not change the underlying amino acid, it's a synonymous variant, and the impact of this variant should be very low if any?

However, snpEff generated the following results:

    chr19	7975586	.	C	T	275.78	.	AC=2;AF=1.00;AN=2;DP=9;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=30.64;SOR=4.615;ANN=T|structural_interaction_variant|HIGH|MAP2K7|ENSG00000076984|interaction|2DYL:A_166-A_207:ENST00000397983|protein_coding|7/12|c.621C>T||||||,T|synonymous_variant|LOW|MAP2K7|ENSG00000076984|transcript|ENST00000545011|protein_coding|6/11|c.699C>T|p.Asp233Asp|764/3497|699/1386|233/461||,T|synonymous_variant|LOW|MAP2K7|ENSG00000076984|transcript|ENST00000397981|protein_coding|6/11|c.573C>T|p.Asp191Asp|675/3430|573/1281|191/426||,T|synonymous_variant|LOW|MAP2K7|ENSG00000076984|transcript|ENST00000397983|protein_coding|7/12|c.621C>T|p.Asp207Asp|681/3396|621/1308|207/435||,T|synonymous_variant|LOW|MAP2K7|ENSG00000076984|transcript|ENST00000397979|protein_coding|6/11|c.573C>T|p.Asp191Asp|627/3361|573/1260|191/419||,T|upstream_gene_variant|MODIFIER|MAP2K7|ENSG00000076984|transcript|ENST00000498118|retained_intron||n.-398C>T|||||398|,T|downstream_gene_variant|MODIFIER|CTD-3193O13.13|ENSG00000268149|transcript|ENST00000595655|antisense||n.*1947G>A|||||1947|,T|downstream_gene_variant|MODIFIER|MAP2K7|ENSG00000076984|transcript|ENST00000475022|retained_intron||n.*3635C>T|||||3635|,T|non_coding_transcript_exon_variant|MODIFIER|MAP2K7|ENSG00000076984|transcript|ENST00000468058|retained_intron|5/10|n.768C>T||||||,T|non_coding_transcript_exon_variant|MODIFIER|MAP2K7|ENSG00000076984|transcript|ENST00000494348|retained_intron|2/3|n.496C>T||||||,T|non_coding_transcript_exon_variant|MODIFIER|MAP2K7|ENSG00000076984|transcript|ENST00000465324|retained_intron|2/2|n.365C>T||||||	GT:AD:DP:GQ:PL	1/1:0,9:9:26:304,26,0

snpEff annotated this variant as structural_interaction_variant with HIGH impact. How and why? As said above, shouldn't this simply be a meaningless variant? Is GRCh 37.75 the right match for hg19?

Another example is like the following:

GATK output --

    chr16	3070708	.	TG	T	173.77	.	AC=2;AF=1.00;AN=2;DP=10;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=24.82;SOR=2.584	GT:AD:DP:GQ:PL	1/1:0,7:7:21:211,21,0

snpEff output --
   

     chr16	3070708	.	TG	T	173.77	.	AC=2;AF=1.00;AN=2;DP=10;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=24.82;SOR=2.584;ANN=T|frameshift_variant|HIGH|TNFRSF12A|ENSG00000006327|transcript|ENST00000575124|protein_coding|1/3|c.318delG|p.Ser107fs|341/1133|318/681|106/226||INFO_REALIGN_3_PRIME,T|upstream_gene_variant|MODIFIER|CLDN6|ENSG00000184697|transcript|ENST00000328796|protein_coding||c.-4687delC|||||2521|,T|upstream_gene_variant|MODIFIER|THOC6|ENSG00000131652|transcript|ENST00000326266|protein_coding||c.-3619delG|||||3323|,T|upstream_gene_variant|MODIFIER| ... (omitted)

The question is: this is a variant (indel) located either in an intronic region or in the 5' UTR of a TNFRSF12A transcript, how can this be interpreted as a frameshift variant?

Thanks for any suggestions!


In the second example, the ensEMBL transcript ENST00000575124 has a protein coding region that extends into the intronic regions of other transcripts of this gene, and so explains the frameshift. You can see the extended protein coding region of transcript ENST00000575124 in the [UCSC genome browser][1] .

Not sure about the first example. A [structural_interaction_variant][2] suggests this amino acid is important for the protein structure, or for an interaction with another protein and so is flagged by SNPeff as an important amino acid site. However, in your case there is a synonymous change at this amino acid and so I would think your variant probably has no impact.

  [1]: http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr16%3A3070208-3071207&hgsid=601505481_CnDb4BMbTN9pxLc19NWtMW0DGQHp
  [2]: http://www.sequenceontology.org/browser/current_svn/term/SO:0002093