Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We appreciate your help!

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.

Meaning of warning (<10 samples for inbreeding calculation) while joint genotyping (GenotypeGVCFs)

sp580sp580 GermanyMember

Hello!
I have combined my gvcfs produced by HaplotypeCaller into one vcf file using CombineGVCFs.

When I use said vcf file for joint genotypeing with GenotypeGVCFs, I get often the following warnings:

WARN InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples

This is makes me wonder if I have done the previous steps correctly, so I will give you some background to what I am doing.

1/ We are working with 6 different lines of mouse
2/ For each line, 10 animals were sequenced.
3/ Variants were called per sample with HaplotypeCaller in GVCF mode (after preparing BAMs as per GATK's best practices)
4/ gvcf files were combined with CombineGVCFs regardless of group membership (60 gvcfs)
5/ Joint genotyping with GenotypeGVCFs

I am now wondering if the fact that on step 4 I did not produce one combined vcf per group (i.e. 6 vcfs, each resulting from 10 gvcfs) could have something to do with this warning.

Answers

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @sp580,

    4/ gvcf files were combined with CombineGVCFs regardless of group membership (60 gvcfs)

    This WARN occurs for cohorts with < 10 samples. Are you grouping the ten samples per mouse line as one sample (10/line = 1 sample)? Or is each animal considered a separate sample (6x10=60 samples) where you have 60 unique sample names? You can look at the column labels in your combined GVCF or final callset VCF to check the number of samples. In either case, given you have at least ten samples, you should not get the WARN, so this is curious.

    On a side note, I think the InbreedingCoeff annotation is not of interest to you given typical mouse strains are inbred by design.

  • sp580sp580 GermanyMember
    edited September 2018

    Hi @shlee,
    thanks for the quick response.

    Yes, I have unique samples, one per animal (6x10=60 samples). All were combined into the same vcf file with CombineGVCFs using the following script:

    GATK=path/to/gatk-4.0.6.0
    GVCF=path/to/gvcfs
    OUT=path/to/output
    REF=path/to/reference_genome_ensembl
    
    FLS=$(ls -1 $GVCF | grep ".gz$" | awk -F- '{print "--variant "$1"-"$2}') # insert flag "--variant" to sample name
    
    cd $GVCF
    $GATK/gatk CombineGVCFs \
            -R $REF/Mus_musculus.GRCm38.dna.primary_assembly.fa \
            $(echo $FLS) \# to pass all 60 files as "--variant sample-name.g.vcf.gz"
            -O $OUT/cohort.g.vcf
    

    Followed by GenotypeGVCFs:

    GATK=path/to/gatk-4.0.6.0
    GVCF=path/to/combinedGVCF
    REF=path/to/reference_genome_ensembl
    OUT=path/to/output
    $GATK/gatk GenotypeGVCFs \
            -R $REF/Mus_musculus.GRCm38.dna.primary_assembly.fa \
            -V $GVCF/cohort.g.vcf \
            -O $OUT/joint_genotypes.vcf
    

    In either case, given you have at least ten samples, you should not get the WARN, so this is curious.

    Does it make sense to figure out where the warnings are coming from?; if so, what would you suggest?

    You can look at the column labels in your combined GVCF or final callset VCF to check the number of samples.

    This is how the the column names looks like, all 60 samples are there:

    #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT H07738-L1 H07739-L1 H07740-L1 H07741-L1 H07742-L1 H07743-L1 H07744-L1 H07745-L1 H07746-L1 H07747-L1 H07748-L1 H07749-L1 H07750-L1 H07751-L1 H07752-L1 H07753-L1 H07754-L1 H07755-L1 H07756-L1 H07757-L1 H07758-L1 H07759-L1 H07760-L1 H07761-L1 H07762-L1 H07763-L1 H07764-L1 H07765-L1 H07766-L1 H07767-L1 H07768-L1 H07769-L1 H07770-L1 H07771-L1 H07772-L1 H07773-L1 H07774-L1 H07775-L1 H07776-L1 H07777-L1 H07778-L1 H07779-L1 H07780-L1 H07781-L1 H07782-L1 H07783-L1 H07784-L1 H07785-L1 H07786-L1 H07787-L1 H07788-L1 H07789-L1 H07790-L1 H07791-L1 H07792-L1 H07793-L1 H07794-L1 H07795-L1 H07796-L1 H07797-L1

    Is it OK to have combined all 60 gvcf files into one vcf instead of doing it group-wise (i.e. 6 combined vcfs)?

    EDIT: I added the command to execute GenotypeGVCFs

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭
    edited September 2018

    @sp580,

    Can you check that the annotation is calculated and present in your callset? You can grep 'InbreedingCoeff=' $OUT/joint_genotypes.vcf to pull out any such records. Here is an example record:

    chr18   12197179    rs8095123   G   C   5070.50 PASS    AC=5;AC_Orig=21;AF=0.417;AF_Orig=0.438;AN=12;AN_Orig=48;BaseQRankSum=0.00;ClippingRankSum=0.00;DB;DP=174;DP_Orig=726;ExcessHet=8.9104;FS=38.985;InbreedingCoeff=-0.2699;MQ=54.71;MQRankSum=-3.445e+00;QD=9.94;ReadPosRankSum=-3.600e-02;SOR=1.924;VQSLOD=-5.680e+00;culprit=MQRankSum   GT:AD:DP:GQ:PL  0/1:20,6:26:99:142,0,645    0/0:34,0:34:99:0,102,1016   0/0:27,0:27:81:0,81,903 0/1:33,10:43:99:232,0,1013  0/1:23,9:32:99:194,0,696    1/1:0,12:12:36:411,36,0
    

    If the annotation is present, then the WARN should not have been emitted by GenotypeGVCFs and you can safely ignore it. I can convey to the developers that the WARN is being emitted unexpectedly and therefore causing unnecessary alarm.

    If the annotation is absent, then this is a bug in GenotypeGVCFs and we'll have to get to the root of it. For starters, you are using GATK4.0.6.0 and we'd ask that you use the latest v4.0.8.1 to see if the bug recapitulates. If it does, we would then ask if you could please submit a bug report, following directions in https://software.broadinstitute.org/gatk/guide/article?id=1894.

    P.S. Even if you choose not to submit a bug report, please let us know if your callset contains values for the InbreedingCoeff. Thank you.

    Post edited by shlee on
  • sp580sp580 GermanyMember
    edited September 2018

    Morning @shlee
    The job is still running, but here are two matches to 'InbreedingCoeff='

    grep "InbreedingCoeff=" joint_genotypes_after_CombineGVCFs.vcf | head -1
    1       3000097 .       CTTTTTTTTTTTT   C       236.78  .       AC=3;AF=0.044;AN=68;BaseQRankSum=-1.440e+00;ClippingRankSum=0.00;DP=235;ExcessHet=0.1158;FS=0.000;InbreedingCoeff=-0.0067;MLEAC=4;MLEAF=0.059;MQ=26.93;MQRankSum=0.00;QD=23.68;ReadPosRankSum=-1.440e+00;SOR=0.223  GT:AD:DP:GQ:PGT:PID:PL  0/0:2,0:2:6:.:.:0,6,79  ./.:4,0:4:.:.:.:0,0,0   0/0:4,0:4:3:.:.:0,3,45  0/0:4,0:4:0:.:.:0,0,80      0/0:6,0:6:6:.:.:0,6,90  0/0:6,0:6:3:.:.:0,3,45  0/0:1,0:1:3:.:.:0,3,41  0/0:3,0:3:0:.:.:0,0,31  ./.:2,0:2:.:.:.:0,0,0   ./.:2,0:2:.:.:.:0,0,0   0/0:2,0:2:3:.:.:0,3,45      0/0:8,0:8:9:.:.:0,9,135 ./.:2,0:2:.:.:.:0,0,0   ./.:3,0:3:.:.:.:0,0,0   0/0:5,0:5:3:.:.:0,3,45  0/0:2,0:2:3:.:.:0,3,45  0/0:3,0:3:0:.:.:0,0,30  ./.:6,0:6:.:.:.:0,0,0   0/0:5,0:5:6:.:.:0,6,90      0/0:4,0:4:12:.:.:0,12,159       ./.:1,0:1:.:.:.:0,0,0   0/0:3,0:3:3:.:.:0,3,45  0/0:8,0:8:6:.:.:0,6,90  ./.:2,0:2:.:.:.:0,0,0   0/0:4,0:4:9:.:.:0,9,135 0/0:8,0:8:3:.:.:0,3,45      0/0:11,0:11:15:.:.:0,15,225     ./.:3,0:3:.:.:.:0,0,0   ./.:4,0:4:.:.:.:0,0,0   0/0:4,0:4:3:.:.:0,3,45  ./.:1,0:1:.:.:.:0,0,0   ./.:1,0:1:.:.:.:0,0,0   ./.:1,0:1:.:.:.:0,0,0   ./.:2,0:2:.:.:.:0,0,0       0/0:5,0:5:6:.:.:0,6,90  ./.:5,0:5:.:.:.:0,0,0   ./.:2,0:2:.:.:.:0,0,0   ./.:0,0:0:.:.:.:0,0,0   ./.:3,0:3:.:.:.:0,0,0   ./.:5,0:5:.:.:.:0,0,0   ./.:1,0:1:.:.:.:0,0,0       0/0:7,0:7:6:.:.:0,6,90  0/0:3,0:3:0:.:.:0,0,39  0/0:6,0:6:12:.:.:0,12,180       1/1:0,4:4:15:1|1:3000097_CTTTTTTTTTTTT_C:184,15,0       0/1:3,3:6:95:0|1:3000097_CTTTTTTTTTTTT_C:117,0,95   0/0:2,0:2:3:.:.:0,3,45  0/0:3,0:3:3:.:.:0,3,45  0/0:3,0:3:0:.:.:0,0,39  0/0:1,0:1:3:.:.:0,3,36  0/0:5,0:5:6:.:.:0,6,90  ./.:8,0:8:.:.:.:0,0,0   ./.:2,0:2:.:.:.:0,0,0   0/0:8,0:8:9:.:.:0,9,135     0/0:5,0:5:9:.:.:0,9,135 0/0:2,0:2:6:.:.:0,6,76  ./.:3,0:3:.:.:.:0,0,0   ./.:0,0:0:.:.:.:0,0,0   ./.:5,0:5:.:.:.:0,0,0   ./.:2,0:2:.:.:.:0,0,0
    
    
     grep "InbreedingCoeff=" joint_genotypes_after_CombineGVCFs.vcf | sed -n 25p
    1       3001579 .       A       T       28375.56        .       AC=58;AF=0.483;AN=120;BaseQRankSum=0.129;ClippingRankSum=0.00;DP=1588;ExcessHet=0.0000;FS=0.914;InbreedingCoeff=0.7901;MLEAC=58;MLEAF=0.483;MQ=59.63;MQRankSum=0.00;QD=32.50;ReadPosRankSum=0.109;SOR=0.662 GT:AD:DP:GQ:PL  1/1:0,35:35:99:1350,105,0       0/0:10,0:10:30:0,30,354 1/1:0,20:20:60:763,60,0 0/1:18,10:28:99:310,0,577   0/1:13,8:21:99:262,0,400        0/0:33,0:33:63:0,63,1120        0/0:27,0:27:73:0,73,1038        0/1:11,13:24:99:379,0,354       0/1:15,14:29:99:412,0,467       1/1:0,24:24:72:937,72,0     0/0:23,0:23:63:0,63,859 0/0:34,0:34:99:0,99,1169        0/0:28,0:28:81:0,81,1215        0/0:27,0:27:81:0,81,999 0/0:32,0:32:90:0,90,1350        0/0:16,0:16:42:0,42,630 0/0:18,0:18:51:0,51,765     0/0:37,0:37:99:0,99,1485        0/0:25,0:25:72:0,72,1080        0/0:25,0:25:75:0,75,951 1/1:0,26:26:78:948,78,0 1/1:0,14:14:42:525,42,0 1/1:0,27:27:81:1054,81,0   1/1:0,20:20:60:667,60,0  1/1:0,24:24:72:843,72,0 1/1:0,36:36:99:1232,108,0       1/1:0,31:31:93:1180,93,0        1/1:2,26:28:4:903,4,0   1/1:0,23:23:69:803,69,0 1/1:1,27:28:45:957,45,0 0/0:19,0:19:54:0,54,810     0/0:20,0:20:60:0,60,704 0/0:20,0:20:57:0,57,855 0/0:31,0:31:90:0,90,1350        0/0:36,0:36:99:0,99,1395        0/0:36,0:36:99:0,102,1530       0/0:24,0:24:57:0,57,855     0/0:20,0:20:60:0,60,730 0/0:29,0:29:75:0,75,1125        0/0:18,0:18:51:0,51,765 0/0:9,0:9:27:0,27,280   0/0:15,0:15:39:0,39,585 0/0:28,0:28:81:0,81,1215        1/1:0,43:43:99:1604,129,0   0/0:40,0:40:99:0,104,1436       0/1:16,15:31:99:463,0,525       0/1:7,7:14:99:171,0,196 1/1:0,13:13:39:481,39,0 0/0:21,0:21:60:0,60,777 1/1:0,22:22:66:776,66,0 1/1:0,58:58:99:2139,174,0   1/1:0,38:38:99:1433,114,0       1/1:0,27:27:81:1034,81,0        1/1:0,35:35:99:1363,105,0       1/1:0,22:22:66:837,66,0 1/1:0,28:28:84:1073,84,0        1/1:0,34:34:99:1240,102,0  1/1:0,16:16:48:647,48,0  1/1:0,27:27:81:1039,81,0        1/1:0,27:27:81:902,81,0
    

    If the annotation is present, then the WARN should not have been emitted by GenotypeGVCFs and you can safely ignore it.

    The problem is that I cannot find a way to identify the records for which a warning was emited (unless you meant that the annotation is present in any record). Based on the log file, there is not much information I can use to retrieve those entries:

    17:03:27.643 INFO  ProgressMeter -        Current Locus  Elapsed Minutes    Variants Processed  Variants/Minute
    17:03:27.760 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    17:03:28.469 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    17:03:28.471 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    17:03:28.474 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    17:03:28.477 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    17:03:37.957 INFO  ProgressMeter -            1:3057218              0.2                 56000         325802.4
    17:03:55.019 INFO  ProgressMeter -            1:3076254              0.5                 75000         164377.6
    17:04:05.044 INFO  ProgressMeter -            1:3145003              0.6                143000         229405.6
    17:04:15.125 INFO  ProgressMeter -            1:3232316              0.8                230000         290642.6
    17:04:16.165 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    17:04:25.153 INFO  ProgressMeter -            1:3321442              1.0                319000         332811.7
    17:04:27.999 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    17:04:35.218 INFO  ProgressMeter -            1:3414153              1.1                411000         364927.9
    17:04:45.220 INFO  ProgressMeter -            1:3509331              1.3                504000         389806.3
    17:04:55.225 INFO  ProgressMeter -            1:3603441              1.5                598000         409673.2
    17:05:05.293 INFO  ProgressMeter -            1:3703620              1.6                698000         428878.6
    17:05:15.400 INFO  ProgressMeter -            1:3801761              1.8                796000         443219.5
    17:05:25.403 INFO  ProgressMeter -            1:3899901              2.0                894000         455502.7
    17:05:35.469 INFO  ProgressMeter -            1:3998028              2.1                992000         465633.0
    17:05:45.484 INFO  ProgressMeter -            1:4100382              2.3               1089000         474024.4
    17:05:55.521 INFO  ProgressMeter -            1:4195487              2.5               1184000         480396.0
    17:06:05.558 INFO  ProgressMeter -            1:4287673              2.6               1276000         484817.8
    17:06:15.562 INFO  ProgressMeter -            1:4384784              2.8               1373000         490596.6
    17:06:25.675 INFO  ProgressMeter -            1:4481909              3.0               1470000         495416.6
    17:06:30.658 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    

    In any case, the tool reports warnings for some records and not for others, even though all of them have the same experimental design.

    Thanks!

    EDIT: running GenotypeGVCFs with GATK4.0.8.1

    These are the first few lines of the log file, simlar as before:

    10:54:50.698 INFO  ProgressMeter -        Current Locus  Elapsed Minutes    Variants Processed  Variants/Minute
    10:54:50.866 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    10:54:51.817 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    10:54:51.819 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    10:54:51.822 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    10:54:51.827 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    10:55:07.449 INFO  ProgressMeter -            1:3007059              0.3                  6000          21492.5
    10:55:17.470 INFO  ProgressMeter -            1:3064234              0.4                 63000         141192.3
    10:55:27.485 INFO  ProgressMeter -            1:3128310              0.6                127000         207138.4
    10:55:37.495 INFO  ProgressMeter -            1:3195267              0.8                193000         247451.8
    10:55:43.791 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    10:55:47.599 INFO  ProgressMeter -            1:3268385              0.9                266000         280487.2
    10:55:57.701 INFO  ProgressMeter -            1:3339461              1.1                337000         301777.5
    10:55:59.010 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    10:56:07.708 INFO  ProgressMeter -            1:3412150              1.3                409000         318659.9
    10:56:17.720 INFO  ProgressMeter -            1:3478237              1.5                475000         327507.2
    10:56:27.722 INFO  ProgressMeter -            1:3562392              1.6                557000         344450.9
    10:56:37.796 INFO  ProgressMeter -            1:3643526              1.8                638000         357429.6
    10:56:47.838 INFO  ProgressMeter -            1:3729664              2.0                724000         370838.3
    10:56:57.856 INFO  ProgressMeter -            1:3822796              2.1                817000         385504.6
    10:57:07.957 INFO  ProgressMeter -            1:3911919              2.3                906000         396039.6
    10:57:18.005 INFO  ProgressMeter -            1:4002030              2.5                996000         405683.4
    10:57:28.059 INFO  ProgressMeter -            1:4092372              2.6               1081000         412173.3
    10:57:38.106 INFO  ProgressMeter -            1:4175458              2.8               1164000         417184.4
    10:57:48.155 INFO  ProgressMeter -            1:4256637              3.0               1245000         420949.4
    10:57:58.196 INFO  ProgressMeter -            1:4337736              3.1               1326000         424324.5
    10:58:08.306 INFO  ProgressMeter -            1:4436849              3.3               1425000         432674.8
    10:58:18.366 INFO  ProgressMeter -            1:4524541              3.5               1512000         436853.2
    10:58:19.051 WARN  InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    
    Post edited by sp580 on
  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @sp580,

    WARN InbreedingCoeff - Annotation will not be calculated, must provide at least 10 samples
    

    It appears that I have misunderstood how this particular WARN is used. Apologies. I've double-checked a stdout for a larger cohort where 99% of sites include the inbreeding coefficient annotation. This stdout still emits a number of these WARNs. If we delve into the code (I am learning JAVA on my own time), we see the contexts for which this WARN is emitted. Basically, the tool does not count any sample where the genotype is null (./.) towards the sample count. So for your sites, those lacking the annotation will have fewer than ten genotypes.

  • sp580sp580 GermanyMember

    Thanks @shlee
    that makes sense, I checked my records and saw too that when >50 samples (out of 60) are ./. the annotation is not present. So I guess those sites will be remove later on during filtering, and I should not be worried about them.

    Good thing it was not a bug :smile:

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Yep, glad it's not a bug. 🙂

  • masaki396masaki396 Member

    Hi @shlee ,

    I have encountered the same warning messages. Is there anyway we can suppress the warning messages?

    Thanks.

    Masaki

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @masaki396

    You can use --verbosity option and set it to Error to only see the log for errors.

  • OggieMOggieM Member

    After carefully reading this whole thread, I have two questions - how does one ever have 10+ genotypes for Inbreeding Coefficient if they are all deleted by CombineGVCF?
    More interestingly, how did the OP even have a few locations where this was calculated? As far as I can tell, in my files all the genotypes are wiped clean.

  • Tiffany_at_BroadTiffany_at_Broad Cambridge, MAMember, Administrator, Broadie, Moderator admin

    Hi @OggieM - I will respond to you shortly! Sorry, I'm digging through a bit of a backlog at the moment. If anyone else can chime in, please do!

  • tsato3tsato3 Cambridge, MAMember, Broadie, Dev

    Hello @OggieM,

    What do you mean by "deleted?" If you are referring to the "./." genotype, it means "no call," meaning that for that sample there was not enough data to determine its genotype (e.g. coverage was zero).

    What does OP stand for? If you're seeing that most samples do not have genotypes at most sites, it likely means that something went awry in the previous steps in the pipeline.

    Hope this helps,

Sign In or Register to comment.