Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

FastaAlternateReferenceMaker gives the input sequence back without applying changes from vcf

LisetteMKLisetteMK The NetherlandsMember

Dear all,

I have sequencing reads (mapped with tophat/bowtie2 to reference genome) and from those created vcf files using samtools.
Now I have tried to get the fasta nucleotide sequence of a certain gene, for which I ran the following command with GATK v3.5-0

java -jar GATK -T FastaAlternateReferenceMaker -R Documents/sequencing/Pfalciparum.genome.fasta -o Documents/sequencing/CSP/CSP_NF166_GATK.fa -L Pf3D7_03_v3:221323-222516 -V Documents/sequencing/NF166.vcf

The output file is identical to the reference sequence at this locus, but from visual inspection of the aligned reads (bam file in IGV) I know where changes are expected. Any hints on what I might be doing wrong will be much appreciated.

Cheers,
Lisette

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @LisetteMK
    Hi Lissette,

    Can you post an example record from the VCF and the corresponding sequence from the output of FastaAlternateReferenceMaker? For example, if there is a variant in the VCF at position 10, please post the VCF record at position 10 and positions 1-20 from the output fasta.

    Thanks,
    Sheila

  • LisetteMKLisetteMK The NetherlandsMember

    Hi Sheila - yes sure.

    vcf at one position of base change:

    CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NF135sorted.bam

    (...)
    Pf3D7_03_v3 221615 . C T,<*> 0 .DP=69;I16=0,0,43,3,0,0,2465,147625,0,0,1254,43296,0,0,1017,24263;QS=0,1,0;VDB=0.42664;SGB=-0.693147;MQSB=0.0762321;MQ0F=0 PL 255,138,0,255,138,255

    fasta:

    1 Pf3D7_03_v3:221323

    (...)TTTACAGCA C TGTTGGCATT

    ...so the C (pos 10) is still in the output although it should really be a T...

    Thank you so much for looking at this. Let me know if you need anything else.

    L

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @LisetteMK
    Hi L,

    Hmm. I am wondering if the issue is that your VCF does not have the genotype (GT) field. That is necessary for FastaAlternateReferenceMaker. Can you try running ValidateVariants on your VCF? Also, does this happen at all variants sites and not just the variant sites with * as an alternate allele?

    Thanks,
    Sheila

  • LisetteMKLisetteMK The NetherlandsMember

    Hi again,

    I ran the validation tool on the vcf file. It said "no failures" for NF166.
    And yes, all variant sites remain the same like in the original sequence, I have checked for two different vcf files now. Apparently, there is always the * among the ALT, except for INDELs.
    Sorry for this lengthy discussion and thanks again for your input.

    Best
    Lisette

  • pdexheimerpdexheimer Member ✭✭✭✭

    I think that the * allele is your problem - in the lingo of VCFs, that's a 'symbolic allele'. FastaAlternateReferenceMaker only considers SNPs and simple insertions or deletions, every other variant gets ignored. And the presence of a symbolic allele will be enough to classify that variant as symbolic rather than a SNP...

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @LisetteMK
    Hi Lisette,

    Can you submit a bug report? Instructions are here.

    Thanks,
    Sheila

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @pdexheimer
    Hi Phillip,

    The tool should not be skipping over the * allele. I have seen cases where the * allele will be present in the output of FastaAlternateReferenceMaker. I'm just not sure what is going on in Lisette's case.

    -Sheila

  • LisetteMKLisetteMK The NetherlandsMember

    Yes, I'll submit one. Thank you!

    L

  • gatrgatr Member

    @Sheila I have the same problem and have been looking for a solution. Interestingly, I did not have this problem when I used the raw SNP VCF to generate alternate reference. When I used the filtered variants the genomes came out all identical. Any update on this?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @gatr @LisetteMK @pdexheimer
    Hi everyone,

    It seems the correct behavior is that FastaAlternateReferenceMaker should indeed be skipping over the * allele. However, we have seen cases where it includes the * allele in the output fasta. We are aware of this bug and are working to fix it, however, it is not high priority right now.

    @gatr
    Yes, if you are using a filtered VCF, the failing sites will not be included in the final output. If you would like to include failing sites, you need to use the unfiltered VCF.

    -Sheila

  • gatrgatr Member

    @Sheila I am struggling to find a solution for a persistent issue for me. I have whole genome data for 14 individuals, which i put through HaplotypeCaller to generate g.vcf files. Then through GenotypeGVCF and then SelectVariants. I can see that the jointgenotype vcf shows variants among the individuals. These same variants can also be viewed in the filtered vcfs split by sample (IGV screenshot attached). But when I run each of the individual filtered .vcf files through FastaAlternateReferenceMaker, each of the output fastas are identical! No variants come through to the sequences. I have seen the previous response on this forum and there are no * symbolic alleles in the vcf file. These filtered .vcfs are the ones I want to use to generate a genome fasta from (not the unfiltered VCF), because whats the point of doing all these filtering if I cannot use it in my intended analyses? I need genome Fasta's with the HQ variants so that I can run phylogenetic analyses on these. What are my options? Unfortunately I have a second different project with 20 genomes, and these issues are relevant there as well. So I really need to figure how to move forward.

    If there is no current way around this, what do you think about hard filtering the g.vcf files for individual samples, and running it through this pipeline? That would skip over the GenotypeGVCF step and go straight to SelectVariants and filter sites based on MQ, QD etc. Please advice. thanks

    PS: I have done the samtools mpileup >> consensus fasta route and I do not like it - the results look just weird. Also the consensus sequences not suitable for part of the analyses I want to do - as the phasing of such sequence data becomes another intensive and time consuming step. Not to mention, having been 'raised' on GATK best practices pipelines, I am unable to 'trust' some of the alternates out there - not sure if that is justified or not.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
    edited August 2016

    @gatr
    Hi,

    Can you please post the exact command you ran? Are you using the latest version of GATK?
    Can you please post some example VCF records and the corresponding section of the output FASTA?

    Thanks,
    Sheila

  • SamDDeckerSamDDecker Ghent, BelgiumMember

    Hi @Sheila,

    I have the exact same problem as described above by @gatr.
    I called variants on 48 samples using HaplotyeCaller/GenotypeGVCFs, hard-filtered, threw out all SNPs that did not PASS the filtering, subsetted the multisample VCF file to get one VCF for each sample to use in FastaAlternateReferenceMaker.

    What seems to happen is that, even though a SNP is genotyped as homozygous REF, the output FASTA contains the ALT allele.

    Attached is an IGV screenshot showing an example site (contig__404:10002), where the output sequence for sample "1" should be A, while instead the output gives me a G (the ALT allele in that position).

    Here is the VCF record of this example site :

    contig_404 10002 . A G 556158.93 PASS AC=0;AF=0.00;AN=2;BaseQRankSum=1.14;ClippingRankSum=0.00;DP=116;ExcessHet=0.0000;FS=0.000;InbreedingCoeff=0.9964;MQ=59.98;MQRankSum=0.00;QD=28.87;ReadPosRankSum=1.47;SOR=1.167 GT:AD:DP:GQ:PL 0/0:116,0:116:99:0,120,1800

    And the full command :

    java -jar $gatk \
    -T FastaAlternateReferenceMaker \
    -R clc_w64b65_cleanmerged_dbg2olc_15_75_23_added.fasta \
    -o Sample1_rbcL.fasta \
    -L S1_rbcL_interval.list \
    -V PASSonly_SNPs_Run2.vcf

    Do you have any suggestions?

    Thanks a lot in advance!
    Sam

  • SamDDeckerSamDDecker Ghent, BelgiumMember

    Oh and I have been using GATK3.8-0 :smile:

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
    edited October 2017

    @SamDDecker
    Hi Sam,

    Can you try using --use_IUPAC_sample {sample_name} in your FastaAlternateReferenceMaker command? You can just test a small section and see if
    that helps.

    Thanks,
    Sheila

  • SamDDeckerSamDDecker Ghent, BelgiumMember

    Hi @Sheila,

    In the example I gave earlier, adding the -IUPAC argument indeed did the trick! It now prints the reference allele on that site! I will extract some more sequences, see if it works for all instances and post an update!

    Thanks for the help!
    Sam

Sign In or Register to comment.