Bug Bulletin: we have identified a bug that affects indexing when producing gzipped VCFs. This will be fixed in the upcoming 3.2 release; in the meantime you need to reindex gzipped VCFs using Tabix.

Pedigree Analysis

delangeldelangel Posts: 71GSA Member mod
edited October 2013 in Methods and Workflows

Workflow

To call variants with the GATK using pedigree information, you should base your workflow on the Best Practices recommendations -- the principles detailed there all apply to pedigree analysis.

But there is one crucial addition: you should make sure to pass a pedigree file (PED file) to all GATK walkers that you use in your workflow. Some will deliver better results if they see the pedigree data.

At the moment there are two of the standard annotations affected by pedigree:

  • Allele Frequency (computed on founders only)
  • Inbreeding coefficient (computed on founders only)

Note that you will need at least 10 founders to compute the inbreeding coefficient.

Trio Analysis

In the specific case of trios, an additional GATK walker, PhaseByTransmission, should be used to obtain trio-aware genotypes as well as phase by descent.

Important note

The annotations mentioned above have been adapted for PED files starting with GATK v.1.6. If you already have VCF files generated by an older version of the GATK or have not passed a PED file while running the UnifiedGenotyper or VariantAnnotator, you should do the following:

  • Run the latest version of the VariantAnnotator to re-annotate your variants.
  • Re-annotate all the standard annotations by passing the argument -G StandardAnnotation to VariantAnnotator. Make sure you pass your PED file to the VariantAnnotator as well!
  • If you are using Variant Quality Score Recalibration (VQSR) with the InbreedingCoefficient as an annotation in your model, you should re-run VQSR once the InbreedingCoefficient is updated.

PED files

The PED files used as input for these tools are based on PLINK pedigree files. The general description can be found here.

For these tools, the PED files must contain only the first 6 columns from the PLINK format PED file, and no alleles, like a FAM file in PLINK.

Post edited by Geraldine_VdAuwera on

Comments

  • sboylesboyle Posts: 21Member

    I've been looking all over for how to add a PED file to my VariantAnnotator run. I don't see an explanation on the VariantAnnotator page or here. I've tried to use the -list function to look at possible annotations, but I don't see ped file as an option. How should I pass a ped file to VariantAnnotator to re-annotate variants from an old version GATK run that wasn't originally run with a ped file?

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,293Administrator, GSA Member admin

    PED files are passed through an engine argument, they're not tool-specific: see http://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_sting_gatk_CommandLineGATK.html

    Geraldine Van der Auwera, PhD

  • sboylesboyle Posts: 21Member

    Thank you for the very prompt response. The program is now running away happily (hopefully :).

  • alirezakjalirezakj Posts: 48Member
    edited April 2013

    Hi Geraldine, can you tell me whats wrong with my PED file? Here it is attached! Ped file (A family with, father(NA00001), mother(NA00002), sonNA00003), daugther(NA00004)) the son has ALS disease for example.

    • (number sign)Family ID
    • (number sign)Individual ID
    • (number sign)Paternal ID
    • (number sign)Maternal ID
    • (number sign)Sex (1=male; 2=female; other=unknown)
    • (number sign)Phenotype (-9=missing; 0=missing; 1=unaffected; 2=affected)

    • FAM001 NA00001 0 0 1 1

    • FAM001 NA00002 0 0 2 1
    • FAM001 NA00003 NA00001 NA00002 1 2
    • FAM001 NA00004 NA00001 NA00002 2 1
    Post edited by alirezakj on
  • alirezakjalirezakj Posts: 48Member

    Can I put the cousin also as: - FAM002 NA00005 0 0 2 1 Can I have only sisters and brothers in a ped file without having mother and father? and last question: is better to put the ped file at the time of variant calling with unifiedgenotyper?

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,293Administrator, GSA Member admin

    Hi @alirezakj,

    Actually I think PED files can only contain trios, so if you want to phase siblings you have to put them in as different families (though obviously with the same parents). So FAM1 would have Mom, Dad and Kid1, FAM2 would have Mom, Dad and Kid2, and so on. Not sure what to do about cousins though.

    Some GATK tools use the PED files and some don't. The simplest is to pass your PED file to every tool; those that can use it will do so, and those that can't will just ignore it.

    Geraldine Van der Auwera, PhD

  • alirezakjalirezakj Posts: 48Member
    edited April 2013

    Thank you so much Geraldine, very helpful. Three more questions:

    1. for a disease with complex genetics if a member of a family dose not show the phenotype, what number should be used in the phenotype field of the Ped file (0=missing or 1=unaffected)?
    2. should I always set -pedValidationType as "SILENT" when I pass the Ped file to walkers, why most people set it as SILENT?
    3. Why in the phenotype field missing can be both set as 0 and -9, why not just 0?

    Thank you so much you are being so helpful.

    Post edited by Geraldine_VdAuwera on
  • alirezakjalirezakj Posts: 48Member

    From your explanation I understand, for a family of four (father "F", mother "M", son "S", daughter "D" affected) the Ped file should look like the following:

    • FAM1 F 0 0 1 0
    • FAM1 M 0 0 2 0
    • FAM1 S F M 1 0
    • FAM2 F 0 0 1 0
    • FAM2 M 0 0 2 0
    • FAM2 D F M 2 2

    Am I right? Thanks

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,293Administrator, GSA Member admin
    edited April 2013

    To your previous three questions:

    1. Set 1 if you know they are unaffected, 0 if you don't know.

    2. By default validation is set to STRICT. Some people choose to use SILENT for various reasons, for example if they are using a BAM file containing a large cohort of individuals, but they are only analyzing one family trio. If they used STRICT, the program would complain that all the other samples are lacking pedigree information. If you don't care about those other samples then this is a bother. These are the validation options:

      STRICT
      Require if a pedigree file is provided at all samples in the VCF or BAM files have a corresponding entry in the pedigree file(s).
      SILENT
      Do not enforce any overlap between the VCF/BAM samples and the pedigree data

    3. This is typically a programmer's decision to allow different values to mean the same thing; in this case I don't know why they chose this. It doesn't really matter, just pick one and be sure to always use the same one in all your work to avoid confusion.

    Post edited by Geraldine_VdAuwera on

    Geraldine Van der Auwera, PhD

  • MaxMax Posts: 20Member
    edited May 2013

    Hi, I just created a ped file in order to use it for the variant calling. However I'm receiving the following error:

    INFO 12:45:21,838 PedReader - Reading PED file /Data/samples.ped with missing fields: []

    INFO 12:45:21,949 PedReader - Phenotype is other? false

    I already saw some comments about this error which might occur due to a missformed ped file. However I could not find the problem with my ped file. It looks like this, separated by tabs.

    FAM001 S00002 0 0 1 1

    FAM001 S00003 0 0 2 1

    FAM001 S00001 S00002 S00003 1 2

    FAM002 ....

    Is there any option to get a more detailed error message or can someone tell me whats wrong about my file ? I have to say that I created it by hand, since there are only 8 trios in it.

    Post edited by Max on
  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,293Administrator, GSA Member admin

    Hmm, any chance that you have spaces instead of tabs somewhere in there?

    Geraldine Van der Auwera, PhD

  • MaxMax Posts: 20Member

    I could not find any spaces, only tabs. Just checked it again.

    I also created a minimal PED file, containing only the first trio. However, this resulted in the same warning message. I also tried to experiment with the values, changed the separator to spaces, without success.

  • CarneiroCarneiro Posts: 271Administrator, GSA Member admin

    Hey MAX that doesn't look like an error, it's just an info line, what is the output of the program?

  • MaxMax Posts: 20Member

    Hi Carneiro,

    well the output looks normal, but I wasnt sure if everything worked fine since I saw this warning message. Therefore I thougt there might be a problem with my ped file.

  • CarneiroCarneiro Posts: 271Administrator, GSA Member admin
  • FarahFarah Posts: 8Member

    Hi,

    Its my first time on the GATK pipeline and I'm happy to say I've made it to having a vcf files. My issue is I now need to filter variants based on inheritance. I have trios but they are not what is usually meant my trios. I have two affected siblings and a parent. So my trio is like this; child1, child2, parent. I need to find homozygous or heterozygous variants that both children inherited from the parent. Can I do this by the GATK -T PhasebyTransmission or no? If no, do I use Annovar or any suggestions as to what I'd use?

    Thanks very much!

  • CarneiroCarneiro Posts: 271Administrator, GSA Member admin

    There is no tool designed to do specifically that, but you can use the toolkit (CombineVariants and SelectVariants are a good start). Or, you can write your own rod walker to do so very easily.

  • KlausNZKlausNZ Posts: 28Member

    @Geraldine_VdAuwera said: Hi alirezakj, Actually I think PED files can only contain trios, so if you want to phase siblings you have to put them in as different families (though obviously with the same parents). So FAM1 would have Mom, Dad and Kid1, FAM2 would have Mom, Dad and Kid2, and so on.

    Geraldine, is this a confirmed workaround? I'm getting "##### ERROR MESSAGE: Inconsistent values detected for [Father] for field Family_ID value1 [FAM1] value2 [FAM2]" (ID's in [brackets] replaced)" when trying this approach. Performing two successive PhaseByTransmission runs gets around this (commenting out the other trio each time (just commented with '#' at start of line, so problem unlikely caused by ped file format error (though I wish as it's be the easiest fix ;-)). Seems to work

    However, what happens if the phase for a parent genotype in the second run doesn't agree with the phase determined in the first run? Does the second run simply over-write the first? ('academic' question only, have yet to observe this case)

    Thanks

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,293Administrator, GSA Member admin

    Hi Klaus,

    I've never needed to do it myself, but I'm told that's how it's supposed to work. It may be that you need to pass the trios in separately -- that would fit with what you're seeing.

    No idea about your hypothetical, but let me know if you run into the case, I'd be interested to find out what happens.

    Geraldine Van der Auwera, PhD

  • byb121byb121 Posts: 12Member

    @Geraldine_VdAuwera said:

    Hi Geraldine,

    I think my questions fit here other than starting a new topic. I have several questions related to the pedigree files usage with no clear answer can be found within the forum. Hope you could help me on these.

    I noticed that with --pedigree, one can pass several pedigree files a time. I guess a few things are not clear here: 1, It's clear that each pedigree file can only have one phenotype, but how about each pedigree file for sample has different phenotype? Will GATK consider the pedigrees are under a same phenotype or it knows that in different file there is a different phenotype?

    2, The note of this doc mentioned what to do to adapt ped info for the variants called with UnifiedGenotyper or VariantAnnotator, how about variants called with HaplotypeCaller? Can I do the same?

    3, For trios in a bigger family, if I want to do the phasing on each trio later, should I put each trio in a separate ped file while calling variants with HaplotypeCaller?

    4, I use Queue and scala script to execute HaplotypeCaller, could you shed a light on how to take ped files as an input parameter for the scala?

    Many thanks,

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,293Administrator, GSA Member admin

    Hi @byb121,

    1. I'm sorry but I don't understand your first question, can you please clarify what you mean? Perhaps if you describe an example of use case (with hypothetical details) it will help.

    2. The HaplotypeCaller will treat pedigrees the same way as UnifiedGenotyper, so the advice for UG is also applicable for HC.

    3. As I recall you can use a single ped file for calling variants; only the phasing tools require that you break up the family into separate trio peds.

    4. You can add a pedigree input the same way as you would add any input file argument to your script. What are you currently using as a base to develop your scala script? FYI, next week we are holding a workshop that includes details on how to write scala scripts for Queue. The presentations will be online by Monday 21st; be sure to have a look at them as they will include a lot of helpful documentation on this topic.

    Geraldine Van der Auwera, PhD

  • byb121byb121 Posts: 12Member
    edited October 2013

    Hi @Geraldine_VdAuwera,

    Thanks a lot for your reply. Sorry I didn't make myself clear on the first question. Here's an example for it: There are 30 exomes, 10 of them are for studying disease A ( you will certainly have affected and unaffected ones), another 10 are for disease B, 5 are for disease C, the rest 5 are singulars for 5 other different diseases. The situation is actually quite common in our lab now. As far as I know, you can't use 1 pedigree file for two different diseases (Not sure I am correct on it). So when calling variants I will have to supply 3 ped files (the 5 singulars can be without ped info) for disease A, B and C. My question is: Will GATK know that the 3 files are for 3 diseases but not one, so correct annotations (eg: Allele Frequency and Inbreeding coefficient) will be assigned?

    For better variant calling and recalibration results, we need at least 30 samples as suggested in best practise. So there's also a situation where we add 20 or more additional exomes when we really want to call variants from 10 other samples. If some of the additional exomes are related, should the pedigree file be used while calling variants? I guess it's 'Yes', but a confirmation from you expert is a better comfort.

    Cheers

    Post edited by byb121 on
  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,293Administrator, GSA Member admin

    Hi again and sorry for the delayed response. If you mean that you want to calculate frequency/inbreeding per disease group, this is a post-processing step that is performed after variant calling.

    Geraldine Van der Auwera, PhD

  • KIMJIHYEKIMJIHYE koreaPosts: 2Member

    Hi, I'm trying to run PhaseByTransmission, but I keep getting NullPointerException in the stack trace, and a Code exception error message. I've copy-pasted the standard error output below. The vcf file was produced using GATK UnifiedGenotyper from recalibrated bam files (with -ped option), and the variants were recalibrated as well. The "analysis ready" vcf file was then annotated using VariantAnnotator (-G StandardAnnotation and -ped parameters as recommended above). I've pasted the ped file contents below as well.

    Can you advise me on anything that I can try? Thank you.


    ped file

    F1     case      dad   mom  1       2
    F1     dad        0       0       1       1
    F1     mom      0       0       2       1
    

    stack trace

    INFO 02:57:49,282 HelpFormatter - -------------------------------------------------------------------------------- INFO 02:57:49,284 HelpFormatter - The Genome Analysis Toolkit (GATK) v2.8-1-g932cd3a, Compiled 2013/12/06 16:47:15 INFO 02:57:49,284 HelpFormatter - Copyright (c) 2010 The Broad Institute INFO 02:57:49,284 HelpFormatter - For support and documentation go to http://www.broadinstitute.org/gatk INFO 02:57:49,287 HelpFormatter - Program Args: -T PhaseByTransmission -R /home/adminrig/src/GATK.2.0/resource.bundle/2.8/b37/human_g1k_v37.fasta -V trio.annot2.vcf -ped input.ped -o trio.annot2.pbt.vcf INFO 02:57:49,287 HelpFormatter - Date/Time: 2014/04/04 02:57:49 INFO 02:57:49,288 HelpFormatter - -------------------------------------------------------------------------------- INFO 02:57:49,288 HelpFormatter - -------------------------------------------------------------------------------- INFO 02:57:49,297 ArgumentTypeDescriptor - Dynamically determined type of trio.annot2.vcf to be VCF INFO 02:57:49,785 GenomeAnalysisEngine - Strictness is SILENT INFO 02:57:49,939 GenomeAnalysisEngine - Downsampling Settings: Method: BY_SAMPLE, Target Coverage: 1000 INFO 02:57:49,959 RMDTrackBuilder - Loading Tribble index from disk for file trio.annot2.vcf INFO 02:57:50,036 PedReader - Reading PED file input.ped with missing fields: [] INFO 02:57:50,113 PedReader - Phenotype is other? false INFO 02:57:50,162 GenomeAnalysisEngine - Preparing for traversal INFO 02:57:50,176 GenomeAnalysisEngine - Done preparing for traversal INFO 02:57:50,177 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING] INFO 02:57:50,178 ProgressMeter - Location processed.sites runtime per.1M.sites completed total.runtime remaining INFO 02:58:20,180 ProgressMeter - 7:95001555 2.63e+04 30.0 s 19.0 m 42.8% 70.0 s 40.0 s

    ERROR ------------------------------------------------------------------------------------------
    ERROR stack trace

    java.lang.NullPointerException at org.broadinstitute.sting.gatk.walkers.phasing.PhaseByTransmission$TrioPhase.getPhasedGenotype(PhaseByTransmission.java:431) at org.broadinstitute.sting.gatk.walkers.phasing.PhaseByTransmission$TrioPhase.getPhasedGenotypes(PhaseByTransmission.java:390) at org.broadinstitute.sting.gatk.walkers.phasing.PhaseByTransmission.phaseTrioGenotypes(PhaseByTransmission.java:753) at org.broadinstitute.sting.gatk.walkers.phasing.PhaseByTransmission.map(PhaseByTransmission.java:853) at org.broadinstitute.sting.gatk.walkers.phasing.PhaseByTransmission.map(PhaseByTransmission.java:126) at org.broadinstitute.sting.gatk.traversals.TraverseLociNano$TraverseLociMap.apply(TraverseLociNano.java:267) at org.broadinstitute.sting.gatk.traversals.TraverseLociNano$TraverseLociMap.apply(TraverseLociNano.java:255) at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.executeSingleThreaded(NanoScheduler.java:274) at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.execute(NanoScheduler.java:245) at org.broadinstitute.sting.gatk.traversals.TraverseLociNano.traverse(TraverseLociNano.java:144) at org.broadinstitute.sting.gatk.traversals.TraverseLociNano.traverse(TraverseLociNano.java:92) at org.broadinstitute.sting.gatk.traversals.TraverseLociNano.traverse(TraverseLociNano.java:48) at org.broadinstitute.sting.gatk.executive.LinearMicroScheduler.execute(LinearMicroScheduler.java:99) at org.broadinstitute.sting.gatk.GenomeAnalysisEngine.execute(GenomeAnalysisEngine.java:313) at org.broadinstitute.sting.gatk.CommandLineExecutable.execute(CommandLineExecutable.java:113) at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:245) at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:152) at org.broadinstitute.sting.gatk.CommandLineGATK.main(CommandLineGATK.java:91)

    ERROR ------------------------------------------------------------------------------------------
    ERROR A GATK RUNTIME ERROR has occurred (version 2.8-1-g932cd3a):
    ERROR
    ERROR This might be a bug. Please check the documentation guide to see if this is a known problem.
    ERROR If not, please post the error message, with stack trace, to the GATK forum.
    ERROR Visit our website and forum for extensive documentation and answers to
    ERROR commonly asked questions http://www.broadinstitute.org/gatk
    ERROR
    ERROR MESSAGE: Code exception (see stack trace for error itself)
    ERROR ------------------------------------------------------------------------------------------
  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,293Administrator, GSA Member admin

    Have you checked if there is a blank line in your file? I think I've seen this happen before, and you do have the line "Reading PED file input.ped with missing fields: [] " in the log output.

    Geraldine Van der Auwera, PhD

  • KIMJIHYEKIMJIHYE koreaPosts: 2Member

    I don't see any unnecessary spaces or blank lines. I don't believe that my ped file format is incorrect. Doesn't "missing fields: []" mean that I don't have any missing fields within the ped file? Here is the output of "cat -A input.ped":

    F1^Icase^Idad^Imom^I1^I2$
    F1^Idad^I0^I0^I1^I1$
    F1^Imom^I0^I0^I2^I1$

    Do you think there could be other reasons for the error? Could it be a bug? Thank you again for your help Geraldine. I really appreciate it.

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,293Administrator, GSA Member admin

    Hmm, that looks fine. Can you please validate your VCF to make sure there is nothing wrong with it?

    Geraldine Van der Auwera, PhD

Sign In or Register to comment.