This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to call denovo and rare variants from a trios cohort?

Hi all, I have two objectives for my analysis:

  1. To find out denovo and rare variants from autism trios analysis Trios: (Father, Mother, Proband)
  2. To detect frequently occuring mutation across autism patient.

I have used GATK best practice to identify the germline variants across 15 trios (Total 45 samples [15X3]) and I have generated gVCF file for all the 45 samples.

Can some one guide me how to call de novo variants from this trios cohort?

gatk gvcf trio denovo variants
• 4,881 views • link
updated by Jeremy Leipzig Jeremy Leipzig • written by DareDevil DareDevil

2 answers

run GATK VariantAnnotator with '-A PossibleDeNovo --pedigree the.ped '

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

You'll need a VCF with the trio(s), not 45 gVCFs, so GenotypeGVCFs those together (or at least into 15 trio vcfs)

Annotate with snpEff

java -Xmx4G -jar snpEff.jar -i vcf -o vcf GRCh37.69 trio.vcf > trio_snpeff.vcf

and then to get variant frequencies, annotate with gnomAD using VCFanno vcfanno: issues to annotate VCF with gnomAD r3.0

Gemini is a nice tool to query VCFs

gemini load -v trio_snpeff.vcf -p trio.ped trio.db
gemini de_novo --columns "chrom,start,end,etc..." trio.db

Bear in mind in this approach you'll detect a ton of de novo mutations, but only 50 or so in a genome or just 1 or 2 in an exome are real. So you will likely have to use a lot of filters to get rid of artefacts. In most papers they still Sanger the candidates.

• 0 views • link
written by Jeremy Leipzig Jeremy Leipzig

@Jeremy Leipzig

Are there much difference between gVCF and VCF.

In GATK forum they appear to be similar for me.

Are you sure the analysis can not be carried out with gVCF ?

Thanks

• 0 views • link
written by DareDevil DareDevil

a gVCF can only handle one individual, so you couldn't use them call trios with tools that expect a trio VCF (correction: yes there are multi-sample gVCFs that briefly arise from CombineGVCFs but they are never found in the wild)

• 0 views • link
written by Jeremy Leipzig Jeremy Leipzig

Are there any tools to convert gvcf to vcf file. I found another question here

They use following command

gzip -dc sample.genome.vcf.gz | extract_variants | bgzip -c > sample.vcf.gz

is it fine to do like this?

@Jeremy Leipzig

@ Pierre Lindenbaum

• 0 views • link
written by DareDevil DareDevil

you said

I have used GATK best practice to identify the germline variants across 15 trios

and then

Are there any tools to convert gvcf to vcf file

?????!!!!

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Hi @Pierre Lindenbaum, I used following command to generate the vcf files. Is it correct to call the vcf files this way

Variant calling using Haplotype caller algorithm

gatk HaplotypeCaller --input sample1_sorted_int.deduped.bqsr.bam --reference $reference --emit-ref-confidence GVCF --dbsnp $DBSNP --output sample1_genome.vcf.gz -L $interval

Calling only variants from a Genome VCF file

gatk GenotypeGVCFs -R $reference --variant sample1_genome.vcf.gz -O sample1_genotypeGvcf.vcf -L $interval
• 0 views • link
written by DareDevil DareDevil

Is it correct to call the vcf files this way

no you first need to call CombineGVCFs for all your samples and then GenotypeGVCFs

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

@Pierre Lindenbaum

Can you clarify my doubt ? there are gvcf files from mother, father and proband. (15*3=45 files)

How I need to combine the gvcf file:

  1. Combine all 45 gvcf files

  2. Combine 15 proband's gvcf files

  3. Combine gvcf of single trios, similarly for all 15 trios separately ?
• 0 views • link
written by DareDevil DareDevil

@Pierre Lindenbaum, I appreciate your help

Hi, I tried combinegvcf followed by GenotypeGVCFs

gatk CombineGVCFs \
    -R $reference \
    -L $interval \
    --variant 001C1.vcf.gz \
    --variant 001c2.vcf.gz \
    --variant 001P1.vcf.gz \
        .....
        .....
        -O cohort.g.vcf.gz

Next

gatk --java-options "-Xmx16g" GenotypeGVCFs  --reference $reference --variant cohort.g.vcf.gz --output output.vcf.gz --dbsnp $DBSNP --intervals $interval

It gave me an output like this,

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  001C1   001P1   001P2   001c2   003C    003P1   003P2   004C    004P1   004P2   007C    007P1   007P2   009C    009P1   009P2   010C    010P1   010P2   013C    013P1   013P2   022C    022P1   022P2   023C    023P1   023P2   030C    030P1   030P2   034C    034P1   034P2   036C    036P1   036P2   042C    042P1   042P2   044C    044P1   044P2
chr1    13273   .   G   C   1660.27 .   AC=11;AF=0.196;AN=56;BaseQRankSum=0.464;DP=295;ExcessHet=0.0920;FS=0.000;InbreedingCoeff=0.3988;MLEAC=17;MLEAF=0.304;MQ=42.17;MQRankSum=-1.139e+00;QD=13.61;ReadPosRankSum=0.848;SOR=0.680  GT:AD:DP:GQ:PL  0/0:9,0:9:0:0,0,170 0/1:8,7:15:99:135,0,153 1/1:0,23:23:69:601,69,0 0/1:12,6:18:94:94,0,270 0/0:14,0:14:42:0,42,376 1/1:0,11:11:33:277,33,0 0/0:20,0:20:60:0,60,537 0/0:4,0:4:12:0,12,115   0/1:6,5:11:97:97,0,115  0/0:2,0:2:6:0,6,51  0/0:8,0:8:24:0,24,199   ./.:0,0:0:.:0,0,0   0/0:15,0:15:45:0,45,363 0/0:11,0:11:33:0,33,292 0/0:14,0:14:42:0,42,349 ./.:0,0:0:.:0,0,0   0/0:9,0:9:27:0,27,227   0/0:9,0:9:27:0,27,249   0/0:6,0:6:18:0,18,150   0/0:6,0:6:18:0,18,155   0/0:7,0:7:21:0,21,174   0/1:5,9:14:82:178,0,82  ./.:0,0:0:.:0,0,0   ./.:0,0:0:.:0,0,0   0/0:7,0:7:21:0,21,203   ./.:0,0:0:.:0,0,0   ./.:0,0:0:.:0,0,0   ./.:0,0:0:.:0,0,0   ./.:0,0:0:.:0,0,0   0/0:3,0:3:9:0,9,81  ./.:0,0:0:.:0,0,0   0/0:2,0:2:6:0,6,66  0/0:14,0:14:42:0,42,382 ./.:0,0:0:.:0,0,0   ./.:0,0:0:.:0,0,0   0/0:8,0:8:24:0,24,190   ./.:0,0:0:.:0,0,0   ./.:0,0:0:.:0,0,0   0/0:5,0:5:15:0,15,125   ./.:0,0:0:.:0,0,0   1/1:0,7:7:21:180,21,0   0/1:15,8:23:99:129,0,357    ./.:0,0:0:.:0,0,0

How do I call denovo variants from each trio set here?

• 0 views • link
written by DareDevil DareDevil

see my answer. I told you to use -A PossibleDeNovo

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

May I know whether the pedigree file which I created is correct?

FamID   IndID   PatID   MatID   Sex     Phenotype
1       in1     Pat1    Mat1    1       2
1       Pat1    0       0       1       0
1       Mat1    0       0       2       0
2       in2     Pat2    Mat2    1       2
2       Pat2    0       0       1       0
2       Mat2    0       0       2       0

where Pat is fatherID;

Mat is MotherID;

in is childID

• 0 views • link
written by DareDevil DareDevil

Ya, this format is correct, but you need to add # in the first line as:

#FamID   IndID   PatID   MatID   Sex     Phenotype
1       in1     Pat1    Mat1    1       2
1       Pat1    0       0       1       0
1       Mat1    0       0       2       0
2       in2     Pat2    Mat2    1       2
2       Pat2    0       0       1       0
2       Mat2    0       0       2       0
• 0 views • link
written by svp svp

Log in to answer this question.

More posts like this