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.

supporting dataset for CalculateGenotypePosteriors

Dear team,

I am relatively new to the GATK environment, so please forgive me if I missed something obvious. I realize that similar question came up before, but I did not find an answer that solved my problem.

I am trying to run CalculateGenotypePosteriors with a supporting dataset. In the tool documentation you use 1000G.phase3.integrated.sites_only.no_MATCHED_REV.hg38.vcf.gz, which I downloaded from the GATK bundle

console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/


When I run
gatk CalculateGenotypePosteriors -R Homo_sapiens_assembly38.fasta -V in.vcf.gz -O out.vcf.gz -supporting 1000G.phase3.integrated.sites_only.no_MATCHED_REV.hg38.vcf.gz

the result is a user error

A USER ERROR has occurred: Input files reference and features have incompatible contigs: Found contigs with the same name but different lengths:
contig reference = chr15 / 101991189
contig features = chr15 / 90338345.

All my alignment, variant calling and genotyping has been done with the same Homo_sapiens_assembly38.fasta file (obtained from the GATK bundle). I am using GATK 4.0.6.0.

Running
gatk ValidateVariants -R Homo_sapiens_assembly38.fasta -V in.vcf.gz --dbsnp GATK-bundle/dbsnp_138.hg38.vcf.gz

completed without an error. So my question is: Is there a problem with this supporting input file (1000G.phase3.integrated.sites_only.no_MATCHED_REV.hg38.vcf.gz)? Is there another file I could use? Are there other tests I could use to check for the integrity of my input vcf files?

Best wishes,

Georg

Answers

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @gwotto,

    Looks like you've run into the same issue as that discussed in this thread. Can you tell us from where and when you downloaded this file? If it's been a while, can you download it again?

    If it is no longer available from the original source, does there appear to be another file that can take its place? If not, please check out our cloud GRCh38 files at https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0. You should be able to click and download any file in this bucket.

  • gwottogwotto Member
    Dear @shlee,

    thanks a lot for your answer. It appears to be the same problem as the in the thread from 2016. I downloaded the file a week ago from the repository you indicated. I also tried the file 1000G_phase1.snps.high_confidence.hg38.vcf.gz, but that one does not contain genotypes, so it gave an error too... I suppose CalculateGenotypePosteriors is widely used, so I wonder what other users are using?

    Best wishes,

    Georg
  • SChaluvadiSChaluvadi Member, Broadie, Moderator admin

    @gwotto I am not positive if this is what you are looking for but I have found from the 1000G website, this list of per chromosome phase 3 vcf files. There are a few README files that might be able to help you determine which set of files you could use, if applicable.

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @gwotto,

    I can confirm that 1000G.phase3.integrated.sites_only.no_MATCHED_REV.hg38.vcf.gz from the cloud bucket contains the chr15 contig with the expected length 90338345. In fact, from the header line, we see the following assembly information:

     ##contig=<ID=chr15,assembly=GCF_000001405.26,length=90338345>
    

    The GCF_000001405.26 has my memory tickling and I think these are from NCBI remap. At least, there is a remap_api.pl perl module I've used in the past that uses this naming scheme to indicate references. I believe this particular reference is the GRCh38 reference. But this doesn't help you. I would suggest checking out the NCBI online remap service at https://www.ncbi.nlm.nih.gov/genome/tools/remap for clues as to what may be going on.

    Perhaps @SChaluvadi can help follow up as well with folks who may know the origin of this file on this side.

  • john156john156 Member
    Hey everyone
    I also tried downloading this file from the google cloud bucket, and it seems the issue persists, i.e. the file abruptly cuts of at
    chr15 90276957 rs529800547 G C 100 PASS AC=1;AF=0.000199681;AFR_AF=0.0000;AMR_AF=0.0000;AN=5008;ASP;DP=19063;EAS_AF=0.0000;EUR_AF=0.

    In that previous thread it was mentioned that it might be fixed after GATK4 release, so I was wondering this was taken care of somewhere else, or it's still on the roadmap?
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    HI @john156

    Please post the link of the file in the google bucket you are trying to download and I will look into it for you.

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @john156

    We have not heard from you in 2 business days and so we are not closing this issue. Please write to us if you have any other questions.

  • gwottogwotto Member
    Dear @bhanuGandham, I can not speak for @john156, but I am pretty sure that he refers to the same file that I pointed at at the beginning of this thread (and a related error). If you can not see the whole thread, the file is from
    console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/
    and is called
    1000G.phase3.integrated.sites_only.no_MATCHED_REV.hg38.vcf.gz
    Best wishes,
    Georg
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Oops sorry I didn't scroll up. My bad!
    Thank you @gwotto, I will look into this issue.

  • AdelaideRAdelaideR Unconfirmed, Member, Broadie, Moderator admin
  • gwottogwotto Member
    Dear @AdelaideR, see my response there. I tried so many different files now, converting them back and fort, but they are either from the wrong assembly, the wrong vcf version or there are other problems, as I pointed out. So, is it reasonable to assume that there is no file available that works as support file for CalculateGenotypePosterior with hg38 and GATK 4? This brings me back to my original questions: what are other people using? What is used in the Broad pipelines?
    Best wishes,
    Georg
  • gwottogwotto Member
    edited March 19
    Hi,

    I think I was finally able to solve that problem, after coming across the answer of @hindrik here:

    gatkforums.broadinstitute.org/gatk/discussion/8694/1000g-phase3-integrated-sites-only-no-matched-rev-hg38-vcf-corrupted#latest

    "It turns out that chr15 has been skipped, chr16 replaced chr15 and that chr17 has been used twice."

    Indeed, when I open 1000G.phase3.integrated.sites_only.no_MATCHED_REV.hg38_1.vcf, I see this (had to remove all angle brackets, to be able to show it here):

    ##contig= ID=chr15,assembly=GCF_000001405.26,length=90338345
    ##contig= ID=chr16,assembly=GCF_000001405.26,length=83257441
    ##contig= ID=chr17,assembly=GCF_000001405.26,length=83257441

    So I replace the lengths of chr15 and chr16 with the correct ones (which I have from the vcf file of my aligned data), resulting in these lines:

    ##contig= ID=chr15,assembly=GCF_000001405.26,length=101991189
    ##contig= ID=chr16,assembly=GCF_000001405.26,length=90338345
    ##contig= ID=chr17,assembly=GCF_000001405.26,length=83257441

    As I said, I do not have too much experience in with GATK, so I do not know if this editing has any unwanted side effects, but at least the error I mentioned above does not come up anymore.

    As I mentioned earlier, linking to a broken file can cause a lot of confusion, especially for beginners, so you should think of fixing, or at least removing it....

    Best wishes,

    Georg
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @gwotto

    Thank you for helping us resolve this, I will look into it.

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin
    edited April 8

    Hi @gwotto

    I have a question, are you able to see chr16 variants in your output file? I ask to verify if the error is only with the contig names in the header or if the chr16 variants are missing as well in the -supporting 1000G.phase3.integrated.sites_only.no_MATCHED_REV.hg38.vcf.gz

  • gwottogwotto Member
    Dear @bhanuGandham thanks a lot for pointing this out. I checked, and indeed the file 1000G.phase3.integrated.sites_only.no_MATCHED_REV.hg38.vcf.gz
    contains only variants up to chr15. This has escaped my attention, because I was happy that the program was running after my modifications of the file. The output file however contains all the variants. I suppose the calculation however is affected by chromosomes > chr15 missing.
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    @gwotto

    Thank you for the update. You are right. 1000G.phase3.integrated.sites_only.no_MATCHED_REV.hg38.vcf.gz file needs to be recreated and is missing >chr15 variants. I am looking into fixing this. In the mean time, as mentioned by @SChaluvadi
    use the 1000G website, this list of per chromosome phase 3 vcf files.

Sign In or Register to comment.