Question: SnpEff output summary for hundreds of vcfs
0
gravatar for moorem
3.5 years ago by
moorem220
United Kingdom
moorem220 wrote:

Hi all! 

I was wondering whether anyone is aware of any existing tools to summarise snpEffs vcf output.

I have a few hundred vcf files from snpEff and want to know how often there is, for example a high impact*, mutation in the same gene. So the output would indicate in how many of my vcf files was there a high impact mutation. 

*As per snpEffs definitions. Relevant to microbes: stop_gained, stop_lost, frameshift

Thanks in advance for any help with this! 

ADD COMMENTlink modified 3.5 years ago by Pierre Lindenbaum119k • written 3.5 years ago by moorem220
4
gravatar for Pierre Lindenbaum
3.5 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum119k wrote:

I use my tool "GroupByGene" https://github.com/lindenb/jvarkit/wiki/GroupByGene to group this kind of VCF information to a gene. All samples must be grouped in the same filtered VCF (filtered= keep "high impact" variants ).

ADD COMMENTlink written 3.5 years ago by Pierre Lindenbaum119k

This looks really promising and an impressive suite of tools overall! 

The issue I'm having with the output right now is I'm getting crazy numbers. I'll check I haven't done anything wrong with snpEff but otherwise have you seen this before:

java -jar ../dist-1.139/groupbygene.jar file1_snpEff.vcf file2_snpEff.vcf > out.txt
produces:

#chrom  min.POS max.POS gene.name       gene.type       samples.affected        count.variations        LES1    LES2
NC_002516       2832    13283   PA0006  snpeff-gene-name        2       134     86      93

ADD REPLYlink written 3.5 years ago by moorem220

what's crazy ? everytime you have a variation (synonymous, non-synonymous, etc...) in this gene, the number is incremented.

ADD REPLYlink written 3.5 years ago by Pierre Lindenbaum119k

yes, I too want to know what "crazy" is - right now if someone asked I couldn't even put an order of magnitude on what "crazy" is

ADD REPLYlink written 3.5 years ago by Istvan Albert ♦♦ 80k

I often don't use taht tool. If you suspect an error, please report an issue: https://github.com/lindenb/jvarkit/issues

ADD REPLYlink written 3.5 years ago by Pierre Lindenbaum119k

I suspect it's my issue, otherwise I will do, thanks again for your help! 

ADD REPLYlink written 3.5 years ago by moorem220

The example gene is only 537bp so 86 or 93 would be a crazy amount between this organism and the reference. It must be an error with my data as biologically that is not what would be expected --certainly not for every single gene. The sum of the mutations reported-within genes is also greater than the number of mutations in the entire genome so this seems to be a fair assumption. 

ADD REPLYlink written 3.5 years ago by moorem220

not agree: min.POS =2832 max.POS=13283 so length > 10kb 

ADD REPLYlink written 3.5 years ago by Pierre Lindenbaum119k

thanks, I'll check this! 

ADD REPLYlink written 3.5 years ago by moorem220

Thanks for the reply! 

I've interpreted count.variations, as you say, as counting all the variations in each gene. The total mutations present in all genes with at least one variation ends up as 297,882, but there are only 30,599 mutations overall. 

Have I misunderstood? :/ I think there has to be an issue with the reference that I've used with snpEff. 

ADD REPLYlink written 3.5 years ago by moorem220

could you give an example command line ? I'm trying to look if my genes of interest are present in all of my VCF files. I have few hundreds of annotated VCF files . 

ADD REPLYlink written 3.5 years ago by jan110

Hi , I have an error when I run the program

 

2015-11-17 18:05:30 "Your input file has a malformed header: unexpected tag count 6 in line <ID=GERM,Number=1,Type=Integer,Description="Counts of donor occurs this mutation, total recorded donor number",Source=/opt/local/genomeinfo/qannotate/icgc_germline_qsnp_PUBLIC.vcf,FileDate=null>"
htsjdk.tribble.TribbleException$InvalidHeader: Your input file has a malformed header: unexpected tag count 6 in line <ID=GERM,Number=1,Type=Integer,Description="Counts of donor occurs this mutation, total recorded donor number",Source=/opt/local/genomeinfo/qannotate/icgc_germline_qsnp_PUBLIC.vcf,FileDate=null

 

The annotation was done by in-house qSNP and qAnnonate. Is there a way I could fix this?

 

ADD REPLYlink written 3.4 years ago by jan110
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1668 users visited in the last hour