We've moved!
This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!

Can I use the fasta alternate reference maker to make a reference genome containing IUPAC codes?

amywebbamywebb Member
edited August 2012 in Ask the GATK team

I have been using the fasta alternate reference maker to switch the reference and variant alleles in the human reference genome. I would like to use it to change the reference allele to the IUPAC ambiguity code instead, but I get an error when I include anything but A, G, C, or T as the variant allele in my vcf file. Is there a way to work around this with the fasta alternate reference maker?

Thanks,
Amy Webb

Answers

  • ebanksebanks Broad InstituteMember, Broadie, Dev ✭✭✭✭

    Hi Amy,
    Sorry, no - we don't support IUPAC codes for now. As always: we'd be happy to incorporate a patch if this is something you (or someone else) wanted to implement!

  • manimani Member

    Hi Eric,

    Is there any plan for supporting IUPAC codes in ‘FastaAlternateReferenceMaker’ in the near future? I am generating consensus sequences from vcf files for population genetics studies and I have just gone through the documentation and forum discussions on FastaAlternateReferenceMaker. It appears that you don’t support multi-allelic loci. The documentation is amply clear on this: “if there are multiple variants that start at a site, it chooses one of them randomly”. It would be nice if FastaAlternateReferenceMaker could take multiple alleles and outputs IUPAC codes for such loci.

    Cheers,
    Mani

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Mani,

    Unfortunately that's not something we can devote resources to in the foreseeable future. However, if someone from the community were to implement this, we would be happy to look at a patch.

  • manimani Member

    Hi Geraldine,

    Thanks for the instant reply. Yes, job for me, then!

    Cheers,
    Mani

  • bmoneybmoney Member

    Mani did you ever find an answer for this? Or did you ever write some code? Thanks.

    -Brian

  • volaviivolavii Member

    i`am dealing with the same problem. But isnt it possible to track all heterozygous snps in the vcf-file and exchange the "wrong" nucleotide in the consensus sequence by the respective iupac base? thats my plan... but i dont know if it is a good one^^

  • yeamanyeaman Member

    Will there be any modifications to the alternate reference maker so that instead of a random base, you can specify that it inserts the most common base when there is a polymorphic site? I am working with a highly polymorphic species without a reference genome (lodgepole pine) that we have aligned to a relative (loblolly pine). I tried the alternate reference maker approach using an existing VCF I had which included a large number of individuals, but it resulted in a severely degraded genome because every time there was a low-frequency SNP, I assume there was a 50/50 chance that the minor allele was incorporated into the genome. In many cases, the alternate allele became the reference, and the calls we got from aligning to this adjusted reference were significantly worse. I suppose one alternative would be for me to use a VCF based on only one individual, which would greatly reduce the frequency at which this happens (I probably should have taken this approach at the outset).

  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @yeaman‌

    Hi,

    You can produce a vcf of sites that show variation above a certain frequency threshold, and use that to make the alternate reference. There is no GATK tool that does this exactly, but you can make a script to do it. You may also want to consider using VariantsToTable, since the output is easier to parse in a custom script than a vcf. https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_variantutils_VariantsToTable.php
    Once you have a list of sites, you can convert back to vcf format with VariantsToVCF. https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_variantutils_VariantsToVCF.php

    -Sheila

  • yeamanyeaman Member

    Thanks Sheila. I think that would mostly work, but if I understand correctly, it would still have the problem of substituting the minor allele for the reference in ~50% of the cases where the minor allele is present and passes the frequency threshold. More importantly, it would also fail to adjust the reference for any SNPs where the study species is polymorphic with MAF < cutoff, and both alleles are different than the reference. For cases where there is quite a bit of divergence between the study species and the closest reference, this could be a problem for quite a few SNPs, I think?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @yeaman, it's up to you to tune your thresholds / add logic to your script to be smart about which allele gets selected for each site depending on what you observe in your population and what you're trying to achieve. That's not something GATK can do for you, and we have no plans to add such functionality in the immediate future.

  • yeamanyeaman Member

    Ok, thanks for the help!

  • everestial007everestial007 GreensboroMember ✭✭

    It may not be a priority but here is a list of things that would make "alternatefastamaker" a good and diverse tool.

    • ability to select a sample GT from the vcf file to make the alternate fasta.
    • ability to make a chain file when indels are substituted. The lift over tools from ensemble is quite complicated for non bio-informatic people.
    • ability to substitute rare vs. common allele.
    • ability to substitute at a particular interval.

    I am only making suggestions. Hope it will happen in the future.

    Thanks,

  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @everestial007
    Hi,

    Thanks!
    1) You can use SelectVariants to select the genotypes that you are interested in and use -removeUnusedAlternates to remove the alleles you are not interested in.
    2) What do you mean by chain file for indels?
    3) See suggestion for number 1. You can also use SelectVariants to select on AF.
    4) You should be able to use -L.

    -Sheila

  • everestial007everestial007 GreensboroMember ✭✭

    Hi @Sheila

    2) I think alternatefastamaker doesn't take indels into consideration?
    Chain files are basically a log of the change in base position when InDels are patched into the reference genome (these are helpful when converting the alignments back to reference equivalent). For my one of the scaffold the chain file looks like:

    chain 1000 scaffold_1118 3008 + 0 3008 scaffold_1118 2982 + 0 2982 scaffold_1118
    66 0 1 A . T 66
    130 0 2 A . TG 196
    166 1 0 A C . 362
    157 9 0 A TTGAGAGAT . 520
    77 2 0 C TT . 606
    87 0 1 A . C 695
    188 12 0 A GGTTGAACCGCG . 883
    312 1 0 A C . 1207
    124 1 0 A T . 1332
    1542 4 0 C TTCT . 2875
    129

    For more details see the link: https://genome.ucsc.edu/goldenpath/help/chain.html
    There is only one tool until today (liftOver) http://genomewiki.ucsc.edu/index.php/LiftOver_Howto
    and is very complicated to non-bioinformatics people since it requires using BLAT and converting the chromosomes into small bins.
    I think if InDels could be incorporated in alternatefastamaker along with chain function, GATK would be more funtastic.

    3) I am not sure but I think Select variants will only do partial amount of job when combined with Alternatefastamaker. Say, I want to make a population level genome, for which I have called variants from genome reseq data from 100 individuals. And, I want to select only variants that are present in at least 50% of the sample at each site to make population level genome. I can select the variants using -min50 but when making the alternate reference it will pick the variants randomly, rite? - Not sure if I am putting it very clearly.

    Thanks,

    • Bishwa K.
  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @everestial007
    Hi Bishwa,

    I just tested out FastaAlternateReferenceMaker on a deletion, and I think it works properly/in the way we want it to.
    For example, let's say my reference is CTCTT at positions 1-5, but the sample in the VCF is C (TCTT is deleted). My output fasta from FastaAlternateReferenceMaker is C when I run it on positions 1-5. So, the deleted sites are not included in the output Fasta file. I think that is what the tool should do.

    FastaAlternateReferenceMaker does choose alternate alleles randomly, but if you select the alternate alleles that are represented most frequently, you should not have to worry, right?

    -Sheila

  • everestial007everestial007 GreensboroMember ✭✭

    Hi @Sheila

    I just checked and it says fastaalternatereference maker works only for SNPs and simple indels. I am not sure how much basepair substitution is allowed.

    I realized that when I made alternate genome using fastaAlternate there was no change in the size of the final genome for two different population. Its my guess, that there is equal amount of insertion vs. deletion allowed. - So, its working to some extent. But, when I used the other tools for making population level genome, my genome size changed by almost 3 Mb.
    I hope big indels could be incorporated in later version. Also, creation of a chain file should be a priority if possible.

    I will give a try by pooling most common variants.

    And thanks for the update !

  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @everestial007
    Hi Bishwa,

    If you could post example regions of where the tools differ, that would be useful.

    Thanks,
    Sheila

Sign In or Register to comment.