Heads up:
We’re moving the GATK website, docs and forum to a new platform. Read the full story and breakdown of key changes on this blog.
Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.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.

Require Genotypes in VCF file in order to output Fst statistics

macmac bodøMember

Hello,
I am using a VCF file obtained from "GenotypeGVCFs" to calculate Fst statistics with VCF tools.
I get this error message:
Error: Require Genotypes in VCF file in order to output Fst statistics.
I suppose something went wrong when creating the VCF file with GATK:
I first used "HaplotypeCaller" for -ERC GVCF cohort analysis workflow; and then I combined all the resulting gVCF files into one with "GenotypeGVCFs".
So, an idea about what this error message means?
Thank you.

Best Answer

Answers

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    @mac, please post some example records from the HaplotypeCaller results and the GenotypeGVCFs results.

  • macmac bodøMember

    @shlee Hi :)
    For the haplotype caller, I used the command:
    java -jar GenomeAnalysisTK.jar -R MasurcaAssembly.fasta -T HaplotypeCaller -I CF_160_12_realigned.bam --emitRefConfidence GVCF -o CF_160_12.raw.snps.indels.g.vcf

    For the GenotypeGVCFs I used the command:
    java -jar GenomeAnalysisTK.jar -T GenotypeGVCFs -R MasurcaAssembly.fasta --variant CF_160_12.raw.snps.indels.g.vcf --variant...(70 more)... -o output.vcf

    What do you mean by "example records of the results"?
    I cannot upload the resulting files here apparently.

    Thank you.

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭
    edited July 2017

    @mac,

    You can post a few records directly to the forum using copy-paste.

  • macmac bodøMember
    edited July 2017

    For the results of HaplotypeCaller:

    ##fileformat=VCFv4.2
    ##ALT=<ID=NON_REF,Description="Represents any possible alternative allele at this location">
    ##FILTER=<ID=LowQual,Description="Low quality">
    ##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
    ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are fil
    tered)">
    ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
    ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
    ##FORMAT=<ID=MIN_DP,Number=1,Type=Integer,Description="Minimum DP observed within the GVCF block">
    ##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternat
    e alleles are phased in relation to one another">
    ##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a give
    n sample (but not across samples) connects records within a phasing group">
    ##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in th
    e VCF specification">
    ##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact T
    est to detect strand bias.">
    ##GATKCommandLine.HaplotypeCaller=<ID=HaplotypeCaller,Version=3.7-0-gcfedb67,Date="Sat Jul 15 22:03:53 CEST 2017",Epoc
    h=1500149033890,CommandLineOptions="analysis_type=HaplotypeCaller input_file=[CF_176_13_realigned.bam] showFullBamList
    =false read_buffer_size=null read_filter=[] disable_read_filter=[] intervals=null excludeIntervals=null interval_set_r
    ule=UNION interval_merging=ALL interval_padding=0 reference_sequence=MasurcaAssembly.fasta nonDeterministicRandomSeed=
    false disableDithering=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=
    null downsample_to_coverage=500 baq=OFF baqGapOpenPenalty=40.0 refactor_NDN_cigar_string=false fix_misencoded_quality_
    scores=false allow_potentially_misencoded_quality_scores=false useOriginalQualities=false defaultBaseQualities=-1 perf
    ormanceLog=null BQSR=null quantize_quals=0 static_quantized_quals=null round_down_quantized=false disable_indel_quals=
    false emit_original_quals=false preserve_qscores_less_than=6 globalQScorePrior=-1.0 secondsBetweenProgressUpdates=10 v
    alidation_strictness=SILENT remove_program_records=false keep_program_records=false sample_rename_mapping_file=null un
    safe=null disable_auto_index_creation_and_locking_when_reading_rods=false no_cmdline_in_header=false sites_only=false 
    never_trim_vcf_format_field=false bcf=false bam_compression=null simplifyBAM=false disable_bam_indexing=false generate
    _md5=false num_threads=1 num_cpu_threads_per_data_thread=1 num_io_threads=0 monitorThreadEfficiency=false num_bam_file
    _handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_w
    ith_unindexed_bam=false generateShadowBCF=false variant_index_type=DYNAMIC_SEEK variant_index_parameter=-1 reference_w
    indow_stop=0 phone_home= gatk_key=null tag=NA logging_level=INFO log_to_file=null help=false version=false likelihoodC
    alculationEngine=PairHMM heterogeneousKmerSizeResolution=COMBO_MIN dbsnp=(RodBinding name= source=UNBOUND) dontTrimAct
    iveRegions=false maxDiscARExtension=25 maxGGAARExtension=300 paddingAroundIndels=150 paddingAroundSNPs=20 comp=[] anno
    tation=[StrandBiasBySample] excludeAnnotation=[ChromosomeCounts, FisherStrand, StrandOddsRatio, QualByDepth] group=[St
    andardAnnotation, StandardHCAnnotation] debug=false useFilteredReadsForAnnotations=false emitRefConfidence=GVCF bamOut
    put=null bamWriterType=CALLED_HAPLOTYPES emitDroppedReads=false disableOptimizations=false annotateNDA=false useNewAFC
    alculator=false heterozygosity=0.001 indel_heterozygosity=1.25E-4 heterozygosity_stdev=0.01 standard_min_confidence_threshold_for_calling=-0.0 standard_min_confidence_threshold_for_emitting=30.0 max_alternate_alleles=6 max_genotype_count=1024 max_num_PL_values=100 input_prior=[] sample_ploidy=2 genotyping_mode=DISCOVERY alleles=(RodBinding name= source=UNBOUND) contamination_fraction_to_filter=0.0 contamination_fraction_per_sample_file=null p_nonref_model=null exactcallslog=null output_mode=EMIT_VARIANTS_ONLY allSitePLs=true gcpHMM=10 pair_hmm_implementation=VECTOR_LOGLESS_CACHING pair_hmm_sub_implementation=ENABLE_ALL always_load_vector_logless_PairHMM_lib=false phredScaledGlobalReadMismappingRate=45 noFpga=false sample_name=null kmerSize=[10, 25] dontIncreaseKmerSizesForCycles=false allowNonUniqueKmersInRef=false numPruningSamples=1 recoverDanglingHeads=false doNotRecoverDanglingBranches=false minDanglingBranchLength=4 consensus=false maxNumHaplotypesInPopulation=128 errorCorrectKmers=false minPruning=2 debugGraphTransformations=false allowCyclesInKmerGraphToGeneratePaths=false graphOutput=null kmerLengthForReadErrorCorrection=25 minObservationsForKmerToBeSolid=20 GVCFGQBands=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 70, 80, 90, 99] indelSizeToEliminateInRefModel=10 min_base_quality_score=10 includeUmappedReads=false useAllelesTrigger=false doNotRunPhysicalPhasing=false keepRG=null justDetermineActiveRegions=false dontGenotype=false dontUseSoftClippedBases=false captureAssemblyFailureBAM=false errorCorrectReads=false pcr_indel_model=CONSERVATIVE maxReadsInRegionPerSample=10000 minReadsPerAlignmentStart=10 mergeVariantsViaLD=false activityProfileOut=null activeRegionOut=null activeRegionIn=null activeRegionExtension=null forceActive=false activeRegionMaxSize=null bandPassSigma=null maxReadsInMemoryPerSample=30000 maxTotalReadsInMemory=10000000 maxProbPropagationDistance=50 activeProbabilityThreshold=0.002 min_mapping_quality_score=20 filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">
    ##GVCFBlock0-1=minGQ=0(inclusive),maxGQ=1(exclusive)
    
    Post edited by shlee on
  • macmac bodøMember
    edited July 2017

    For result of GenotypeGVCFs:

    ##fileformat=VCFv4.2
    ##ALT=<ID=NON_REF,Description="Represents any possible alternative allele at this location">
    ##FILTER=<ID=LowQual,Description="Low quality">
    ##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
    ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are fil
    tered)">
    ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
    ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
    ##FORMAT=<ID=MIN_DP,Number=1,Type=Integer,Description="Minimum DP observed within the GVCF block">
    ##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternat
    e alleles are phased in relation to one another">
    ##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a give
    n sample (but not across samples) connects records within a phasing group">
    ##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in th
    e VCF specification">
    ##FORMAT=<ID=RGQ,Number=1,Type=Integer,Description="Unconditional reference genotype confidence, encoded as a phred qu
    ality -10*log10 p(genotype call is wrong)">
    ##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact T
    est to detect strand bias.">
    ##GATKCommandLine.GenotypeGVCFs=<ID=GenotypeGVCFs,Version=3.7-0-gcfedb67,Date="Wed Jul 19 08:34:16 CEST 2017",Epoch=15
    00446056658,CommandLineOptions="analysis_type=GenotypeGVCFs input_file=[] showFullBamList=false read_buffer_size=null 
    read_filter=[] disable_read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=AL
    L interval_padding=0 reference_sequence=MasurcaAssembly.fasta nonDeterministicRandomSeed=false disableDithering=false 
    maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=1
    000 baq=OFF baqGapOpenPenalty=40.0 refactor_NDN_cigar_string=false fix_misencoded_quality_scores=false allow_potential
    ly_misencoded_quality_scores=false useOriginalQualities=false defaultBaseQualities=-1 performanceLog=null BQSR=null qu
    antize_quals=0 static_quantized_quals=null round_down_quantized=false disable_indel_quals=false emit_original_quals=fa
    lse preserve_qscores_less_than=6 globalQScorePrior=-1.0 secondsBetweenProgressUpdates=10 validation_strictness=SILENT 
    remove_program_records=false keep_program_records=false sample_rename_mapping_file=null unsafe=null disable_auto_index
    _creation_and_locking_when_reading_rods=false no_cmdline_in_header=false sites_only=false never_trim_vcf_format_field=
    false bcf=false bam_compression=null simplifyBAM=false disable_bam_indexing=false generate_md5=false num_threads=1 num
    _cpu_threads_per_data_thread=1 num_io_threads=0 monitorThreadEfficiency=false num_bam_file_handles=null read_group_bla
    ck_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false gene
    rateShadowBCF=false variant_index_type=DYNAMIC_SEEK variant_index_parameter=-1 reference_window_stop=0 phone_home= gat
    k_key=null tag=NA logging_level=INFO log_to_file=null help=false version=false variant=[(RodBindingCollection [(RodBin
    ding name=variant source=CF_160_12.raw.snps.indels.g.vcf)]), (RodBindingCollection [(RodBinding name=variant2 source=C
    F_160_13.raw.snps.indels.g.vcf)]), (RodBindingCollection [(RodBinding name=variant3 source=CF_160_14.raw.snps.indels.g
    .vcf)]), (RodBindingCollection [(RodBinding name=variant4 source=CF_176_11.raw.snps.indels.g.vcf)]), (RodBindingCollec
    tion [(RodBinding name=variant5 source=CF_176_13.raw.snps.indels.g.vcf)]), (RodBindingCollection [(RodBinding name=varNP13.raw.snps.indels.g.vcf)]), (RodBindingCollection [(RodBinding name=variant12 source=CF_NP14.raw.snps.indels.g.vcf)
    ]), (RodBindingCollection [(RodBinding name=variant13 source=CF_US12.raw.snps.indels.g.vcf)]), (RodBindingCollection [
    (RodBinding name=variant14 source=CF_US1.raw.snps.indels.g.vcf)]), (RodBindingCollection [(RodBinding name=variant15 s
    ource=CF_US2.raw.snps.indels.g.vcf)]), (RodBindingCollection [(RodBinding name=variant16 source=CF_Is7.raw.snps.indels
    .g.vcf)]), (RodBindingCollection [(RodBinding name=variant17 source=CF_Is8.raw.snps.indels.g.vcf)]), (RodBindingCollec
    tion [(RodBinding name=variant18 source=CF_Lure17.raw.snps.indels.g.vcf)]), (RodBindingCollection [(RodBinding name=va
    riant19 source=CF_Lure18.raw.snps.indels.g.vcf)]), (RodBindingCollection [(RodBinding name=variant20 source=CF_Lure19.
    raw.snps.indels.g.vcf)]), (RodBindingCollection [(RodBinding name=variant21 source=CF_Skj33.raw.snps.indels.g.vcf)]), 
    (RodBindingCollection [(RodBinding name=variant22 source=CF_Skj34.raw.snps.indels.g.vcf)]), (RodBindingCollection [(RodBinding name=variant23 source=CF_Skj35.raw.snps.indels.g.vcf)])] out=/home/mac/Cfin_Captures_For_LA/AllTogether.vcf includeNonVariantSites=false uniquifySamples=false annotateNDA=false useNewAFCalculator=false heterozygosity=0.001 indel_heterozygosity=1.25E-4 heterozygosity_stdev=0.01 standard_min_confidence_threshold_for_calling=10.0 standard_min_confidence_threshold_for_emitting=30.0 max_alternate_alleles=6 max_genotype_count=1024 max_num_PL_values=100 input_prior=[] sample_ploidy=2 annotation=[] group=[StandardAnnotation] dbsnp=(RodBinding name= source=UNBOUND) filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">
    ##GATKCommandLine.HaplotypeCaller=<ID=HaplotypeCaller,Version=3.7-0-gcfedb67,Date="Sat Jul 15 15:00:43 CEST 2017",Epoch=1500123643823,CommandLineOptions="analysis_type=HaplotypeCaller input_file=[CF_160_13_realigned.bam] showFullBamList=false read_buffer_size=null read_filter=[] disable_read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=MasurcaAssembly.fasta nonDeterministicRandomSeed=false disableDithering=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=
    ##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as 
    listed">
    ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
    ##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
    ##INFO=<ID=BaseQRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt Vs. Ref base quali
    ties">
    ##INFO=<ID=ClippingRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref number 
    of hard clipped bases">
    ##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
    ##INFO=<ID=DS,Number=0,Type=Flag,Description="Were any of the samples downsampled?">
    ##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
    ##INFO=<ID=ExcessHet,Number=1,Type=Float,Description="Phred-scaled p-value for exact test of excess heterozygosity">
    ##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
    ##INFO=<ID=HaplotypeScore,Number=1,Type=Float,Description="Consistency of the site with at most two segregating haplot
    ypes">
    ##INFO=<ID=InbreedingCoeff,Number=1,Type=Float,Description="Inbreeding coefficient as estimated from the genotype like
    lihoods per-sample when compared against the Hardy-Weinberg expectation">
    ##INFO=<ID=MLEAC,Number=A,Type=Integer,Description="Maximum likelihood expectation (MLE) for the allele counts (not ne
    cessarily the same as the AC), for each ALT allele, in the same order as listed">
    ##INFO=<ID=MLEAF,Number=A,Type=Float,Description="Maximum likelihood expectation (MLE) for the allele frequency (not n
    ecessarily the same as the AF), for each ALT allele, in the same order as listed">
    ##INFO=<ID=MQ,Number=1,Type=Float,Description="RMS Mapping Quality">
    ##INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">
    ##INFO=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">
    ##INFO=<ID=RAW_MQ,Number=1,Type=Float,Description="Raw data for RMS Mapping Quality">
    ##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias">
    ##INFO=<ID=SOR,Number=1,Type=Float,Description="Symmetric Odds Ratio of 2x2 contingency table to detect strand bias">
    ##contig=<ID=ctg7180000254689,length=2174>
    
    Post edited by shlee on
  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    @mac,

    It seems you've only posted the VCF headers. Can you post variant records?

  • macmac bodøMember

    @shlee ,
    Ok, sorry.
    So, the following comes from HaplotypeCaller:

    contig=<ID=ctg7180000254689,length=2174>

    contig=<ID=ctg7180000254703,length=369>

    contig=<ID=ctg7180000254709,length=468>

    contig=<ID=ctg7180000254717,length=311>

    contig=<ID=ctg7180000254724,length=4625>

    contig=<ID=ctg7180000254726,length=445>

    And the following comes from GenotypeGVCFs:

    contig=<ID=ctg7180000254689,length=2174>

    contig=<ID=ctg7180000254703,length=369>

    contig=<ID=ctg7180000254709,length=468>

    contig=<ID=ctg7180000254717,length=311>

    contig=<ID=ctg7180000254724,length=4625>

    contig=<ID=ctg7180000254726,length=445>

    contig=<ID=ctg7180000254748,length=621>

    contig=<ID=ctg7180000254761,length=323>

    contig=<ID=ctg7180000254778,length=492>

    Is that what you want?
    Thank you for your help.

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭
    edited July 2017

    @mac,

    The contig lines are still a part of the header. Do your VCFs have anything additional after the #CHROM POS ... line? If not, then your error occurs because you have no variant records.

  • macmac bodøMember

    Hi @shlee ,
    As I used the nohup option, I got these information when running HaplotypeCaller. But I do not understand why I get the warning message...
    INFO 09:15:50,245 SAMDataSource$SAMReaders - Initializing SAMRecords in serial
    INFO 09:15:51,463 SAMDataSource$SAMReaders - Done initializing BAM readers: total time 1.22
    INFO 09:15:52,337 HCMappingQualityFilter - Filtering out reads with MAPQ < 20
    INFO 09:29:51,414 GenomeAnalysisEngine - Preparing for traversal over 1 BAM files
    INFO 09:46:04,452 GenomeAnalysisEngine - Done preparing for traversal
    INFO 09:46:04,453 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING]
    INFO 09:46:04,453 ProgressMeter - | processed | time | per 1M | | total | remaining
    INFO 09:46:04,453 ProgressMeter - Location | active regions | elapsed | active regions | completed | runtime | runtime
    INFO 09:46:04,454 HaplotypeCaller - Standard Emitting and Calling confidence set to 0.0 for reference-model confidence output
    INFO 09:46:04,454 HaplotypeCaller - All sites annotated with PLs forced to true for reference-model confidence output
    WARN 09:46:04,492 InbreedingCoeff - Annotation will not be calculated. InbreedingCoeff requires at least 10 unrelated samples.
    INFO 09:46:34,458 ProgressMeter - Starting 0.0 30.0 s 49.6 w 100.0% 30.0 s 0.0 s
    INFO 09:46:50,744 HaplotypeCaller - Using global mismapping rate of 45 => -4.5 in log10 likelihood units
    INFO 09:47:04,462 ProgressMeter - ctg7180000258688:1165 2066681.0 60.0 s 29.0 s 4.0% 25.2 m 24.2 m
    INFO 09:47:34,465 ProgressMeter - ctg7180000264831:815 6756145.0 90.0 s 13.0 s 13.0% 11.6 m 10.1 m
    INFO 09:48:04,469 ProgressMeter - ctg7180000271535:817 1.1023367E7 120.0 s 10.0 s 21.2% 9.5 m 7.5 m
    INFO 09:48:34,474 ProgressMeter - ctg7180000277236:596 1.5365097E7 2.5 m 9.0 s 29.5% 8.5 m 6.0 m
    INFO 09:49:04,481 ProgressMeter - ctg7180000283291:401 1.8815109E7 3.0 m 9.0 s 36.1% 8.3 m 5.3 m
    INFO 09:49:34,486 ProgressMeter - ctg7180000296475:722 2.199497E7 3.5 m 9.0 s 42.2% 8.3 m 4.8 m
    INFO 09:50:04,491 ProgressMeter - ctg7180000313687:571 2.5903419E7 4.0 m 9.0 s 49.7% 8.0 m 4.0 m
    INFO 09:50:34,494 ProgressMeter - ctg7180000327208:380 2.8634046E7 4.5 m 9.0 s 54.9% 8.2 m 3.7 m
    INFO 09:51:04,496 ProgressMeter - ctg7180000387334:210 3.11569E7 5.0 m 9.0 s 59.8% 8.4 m 3.4 m
    INFO 09:51:34,499 ProgressMeter - ctg7180000431408:700 3.5529998E7 5.5 m 9.0 s 68.2% 8.1 m 2.6 m
    INFO 09:52:04,501 ProgressMeter - ctg7180000437416:1300 3.9157402E7 6.0 m 9.0 s 75.1% 8.0 m 119.0 s
    INFO 09:52:34,504 ProgressMeter - ctg7180000443271:584 4.3488883E7 6.5 m 8.0 s 83.4% 7.8 m 77.0 s
    INFO 09:53:04,507 ProgressMeter - ctg7180000449259:353 4.713757E7 7.0 m 8.0 s 90.4% 7.7 m 44.0 s
    INFO 09:53:34,510 ProgressMeter - ctg7180000455080:922 5.1271351E7 7.5 m 8.0 s 98.4% 7.6 m 7.0 s
    Using AVX accelerated implementation of PairHMM
    INFO 09:53:39,503 VectorLoglessPairHMM - libVectorLoglessPairHMM unpacked successfully from GATK jar file
    INFO 09:53:39,503 VectorLoglessPairHMM - Using vectorized implementation of PairHMM
    INFO 09:53:39,504 VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0
    INFO 09:53:39,504 PairHMM - Total compute time in PairHMM computeLikelihoods() : 0.0
    INFO 09:53:39,504 HaplotypeCaller - Ran local assembly on 214607 active regions
    INFO 09:53:39,785 ProgressMeter - done 5.211649E7 7.6 m 8.0 s 100.0% 7.6 m 0.0 s

    INFO 09:53:39,785 ProgressMeter - Total runtime 455.33 secs, 7.59 min, 0.13 hours

    Done. There were 1 WARN messages, the first 1 are repeated below.

    WARN 09:46:04,492 InbreedingCoeff - Annotation will not be calculated. InbreedingCoeff requires at least 10 unrelated samples.

    (END)

  • RaziRazi Member
    edited August 13
    Hi,
    I have the same problem to calculate Fst exactly same as Mac,
    Error: Require Genotypes in VCF file in order to output Fst statistics.
    I have 3+5 samples (for 2 population).

    I used -ERC GVCF for creating vcf files:

    java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R ref.fasta -I input.bam --genotyping_mode DISCOVERY -stand_emit_conf 10 -stand_call_conf 30 -ERC GVCF -variant_index_type LINEAR -variant_index_parameter 128000 -o output.vcf

    and then I used GenotypeGVCFs to merge vcf files:

    java -jar GenomeAnalysisTK.jar -T GenotypeGVCFs -R ref.fasta --variant (8 files) -o list.vcf

    and for Fst calculation:
    vcftools --vcf list.vcf --weir-fst-pop p1.txt --weir-fst-pop p2.txt --out p1_vs_p2


    But I didn't understand how should I solve that? I can't increase my samples.

    This is header of GenotypeGVCFs result file:

    ##fileformat=VCFv4.1
    ##ALT=<ID=NON_REF,Description="Represents any possible alternative allele at this location">
    ##FILTER=<ID=LowQual,Description="Low quality">
    ##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
    ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
    ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
    ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
    ##FORMAT=<ID=MIN_DP,Number=1,Type=Integer,Description="Minimum DP observed within the GVCF block">
    ##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phase$
    ##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not a$
    ##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification$
    .........
    .......................
    ##contig=<ID=scaffold_0,length=20988840>
    ##contig=<ID=scaffold_1,length=17964029>
    ##contig=<ID=scaffold_2,length=16154463>
    ##contig=<ID=scaffold_3,length=13658650>
    ##contig=<ID=scaffold_4,length=13151975>
    ##contig=<ID=scaffold_5,length=12650179>
    ##contig=<ID=scaffold_6,length=12071984>
    ##contig=<ID=scaffold_7,length=11537834>
    ##contig=<ID=scaffold_8,length=11485610>
    ##contig=<ID=scaffold_9,length=11075858>
    ##contig=<ID=scaffold_10,length=10894907>
    ##contig=<ID=scaffold_11,length=10634252>
    .................
    .....................
    ##contig=<ID=scaffold_6451,length=882>
    ##reference=file:///c3se/NOBACKUP/users/raziyeh/sylvia/Database_AllPath/SylviaAllPath.scf.fasta
    #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT ind1
    scaffold_0 12533 . TAAAAAAAA T 172.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.406;ClippingRankSum=0.406;DP=52$
    scaffold_0 12584 . T A 241.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.00;ClippingRankSum=0.133;DP=82;FS=0.000$
    scaffold_0 12613 . C T 270.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=1.94;ClippingRankSum=-8.830e-01;DP=95;FS=$
    scaffold_0 13052 . AATC A 172.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=1.92;ClippingRankSum=0.322;DP=96;FS=0.000$
    scaffold_0 13224 . C T 137.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.198;ClippingRankSum=0.198;DP=92;FS=0.00$
    scaffold_0 13654 . TA T 233.80 . AC=2;AF=1.00;AN=2;DP=89;FS=0.000;GQ_MEAN=21.00;MLEAC=2;MLEAF=1.00;MQ=42.0$
    scaffold_0 13706 . A AT 50.56 . AC=2;AF=1.00;AN=2;BaseQRankSum=0.248;ClippingRankSum=0.365;DP=98;FS=0.000$
    scaffold_0 13812 . G A 69.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.406;ClippingRankSum=0.406;DP=101;FS=0.0$
    scaffold_0 13853 . T C 133.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.00;ClippingRankSum=-1.380e+00;DP=81;FS=$
    scaffold_0 13971 . C T 118.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.278;ClippingRankSum=1.25;DP=87;FS=18.19$
    scaffold_0 13972 . A G 118.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.278;ClippingRankSum=0.597;DP=84;FS=18.1$
    scaffold_0 13989 . A G 114.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=-1.560e-01;ClippingRankSum=0.289;DP=92;FS$
    scaffold_0 14014 . G A 114.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=-5.360e-01;ClippingRankSum=0.825;DP=97;FS
    ...........
    ********************************************************************************
    And this is vcf file example:
    ##fileformat=VCFv4.1
    ##ALT=<ID=NON_REF,Description="Represents any possible alternative allele at this location">
    ##FILTER=<ID=LowQual,Description="Low quality">
    ##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
    ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
    ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
    ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
    ##FORMAT=<ID=MIN_DP,Number=1,Type=Integer,Description="Minimum DP observed within the GVCF block">
    ##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phase$
    ##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not a$
    ##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification$
    ##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect stran$
    ##GATKCommandLine=<ID=HaplotypeCaller,Version=3.3-0-g37228af,Date="Fri Jul 26 21:33:17 CEST 2019",Epoch=1564169597602,CommandLineOptions=$
    ##GVCFBlock=minGQ=0(inclusive),maxGQ=1(exclusive)
    ##GVCFBlock=minGQ=1(inclusive),maxGQ=2(exclusive)
    ##GVCFBlock=minGQ=10(inclusive),maxGQ=11(exclusive)
    ##GVCFBlock=minGQ=11(inclusive),maxGQ=12(exclusive)
    ##GVCFBlock=minGQ=12(inclusive),maxGQ=13(exclusive)
    ..............
    ......................
    contig=<ID=scaffold_2,length=16154463>
    ##contig=<ID=scaffold_3,length=13658650>
    ##contig=<ID=scaffold_4,length=13151975>
    ##contig=<ID=scaffold_5,length=12650179>
    ##contig=<ID=scaffold_6,length=12071984>
    ##contig=<ID=scaffold_7,length=11537834>
    ##contig=<ID=scaffold_8,length=11485610>
    ##contig=<ID=scaffold_9,length=11075858>
    ##contig=<ID=scaffold_10,length=10894907>
    ##contig=<ID=scaffold_11,length=10634252>
    ##contig=<ID=scaffold_12,length=10486095>
    ..................
    .........................
    ##contig=<ID=scaffold_6451,length=882>
    ##reference=file:///c3se/NOBACKUP/users/raziyeh/sylvia/Database_AllPath/SylviaAllPath.scf.fasta
    #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT ind1
    scaffold_0 1 . A <NON_REF> . . END=12441 GT:DP:GQ:MIN_DP:PL 0/0:0:0:0:0,0,0
    scaffold_0 12442 . T <NON_REF> . . END=12445 GT:DP:GQ:MIN_DP:PL 0/0:1:3:1:0,3,37
    scaffold_0 12446 . A <NON_REF> . . END=12451 GT:DP:GQ:MIN_DP:PL 0/0:2:6:2:0,6,83
    scaffold_0 12452 . G <NON_REF> . . END=12454 GT:DP:GQ:MIN_DP:PL 0/0:2:5:2:0,6,62
    scaffold_0 12455 . A <NON_REF> . . END=12455 GT:DP:GQ:MIN_DP:PL 0/0:2:6:2:0,6,88
    scaffold_0 12456 . C <NON_REF> . . END=12457 GT:DP:GQ:MIN_DP:PL 0/0:2:5:2:0,6,62
    scaffold_0 12458 . A <NON_REF> . . END=12464 GT:DP:GQ:MIN_DP:PL 0/0:2:6:2:0,6,73
    scaffold_0 12465 . T <NON_REF> . . END=12483 GT:DP:GQ:MIN_DP:PL 0/0:3:9:3:0,9,119
    scaffold_0 12484 . T <NON_REF> . . END=12486 GT:DP:GQ:MIN_DP:PL 0/0:4:12:4:0,12,174
    scaffold_0 12487 . G <NON_REF> . . END=12491 GT:DP:GQ:MIN_DP:PL 0/0:4:3:4:0,3,45
    scaffold_0 12492 . G <NON_REF> . . END=12532 GT:DP:GQ:MIN_DP:PL 0/0:5:0:4:0,0,0
    scaffold_0 12533 . TAAAAAAAA T,<NON_REF> 163.73 . BaseQRankSum=-0.322;ClippingRankSum=-1.517;DP=8;MLEAC=1,0$
    scaffold_0 12542 . A <NON_REF> . . END=12542 GT:DP:GQ:MIN_DP:PL 0/0:8:0:8:0,0,0
    scaffold_0 12543 . A <NON_REF> . . END=12544 GT:DP:GQ:MIN_DP:PL 0/0:8:24:8:0,24,360
    scaffold_0 12545 . A <NON_REF> . . END=12552 GT:DP:GQ:MIN_DP:PL 0/0:9:27:9:0,27,394
    scaffold_0 12553 . A <NON_REF> . . END=12558 GT:DP:GQ:MIN_DP:PL 0/0:10:30:10:0,30,419
    scaffold_0 12559 . C <NON_REF> . . END=12573 GT:DP:GQ:MIN_DP:PL 0/0:11:33:11:0,33,462
    scaffold_0 12574 . G <NON_REF> . . END=12575 GT:DP:GQ:MIN_DP:PL 0/0:12:36:12:0,36,489
    scaffold_0 12576 . T <NON_REF> . . END=12581 GT:DP:GQ:MIN_DP:PL 0/0:13:39:13:0,39,573
    scaffold_0 12582 . A <NON_REF> . . END=12583 GT:DP:GQ:MIN_DP:PL 0/0:14:42:14:0,42,624
    scaffold_0 12584 . T A,<NON_REF> 241.77 . BaseQRankSum=0.000;ClippingRankSum=0.133;DP=14;MLEAC=1,0;MLEAF=0.$
    scaffold_0 12585 . T <NON_REF> . . END=12585 GT:DP:GQ:MIN_DP:PL 0/0:15:45:15:0,45,666
    scaffold_0 12586 . A <NON_REF> . . END=12589 GT:DP:GQ:MIN_DP:PL 0/0:16:48:16:0,48,676
    scaffold_0 12590 . T <NON_REF> . . END=12597 GT:DP:GQ:MIN_DP:PL 0/0:17:45:16:0,45,67
    ................
    ...........................


    Would you please help me? I need the Fst result soon. Thanks

    Best,
    Razi
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @Razi

    You are using GATK 3 and we do not support that anymore. Please upgrade to the latest version of GATK v4.1.3.0

  • RaziRazi Member
    Hi,
    I post a comment before about this error here. You suggest me to use GATK 4. I used GATK/4.1.2.0:

    1. for creating gvcf:
    gatk HaplotypeCaller -R ref.scf.fasta -I input.bam -ERC GVCF -O out.g.vcf

    2. for combine:
    gatk CombineGVCFs -R ref.scf.fasta --variant 1.g.vcf --variant 2.g.vcf --variant 3.g.vcf --variant 4.g.vcf --variant 5.g.vcf --variant 6.g.vcf --variant 7.g.vcf --variant 8.g.vcf -O outt.vcf

    3. for genotype:
    gatk GenotypeGVCFs -R ref.scf.fasta --variant outt -O outt.vcflist

    4. for Fst:
    vcftools --vcf outt.vcflist --weir-fst-pop 1.txt --weir-fst-pop 2.txt --out 1_vs_2


    There was no error till step 3 and successfully done. But, I got error in step 4:
    Error: Require Genotypes in VCF file in order to output Fst statistics.

    Would you please help me? I really get tired because of Fst and always get error :(

    I really don't know what should I do again?


    Best,
    Razi
  • RaziRazi Member

    Hi @bhanuGandham

    In our server we already have the GATK/4.1.2.0 but you suggested to use GATK v4.1.3.0.

    Are there major changes from v4.1.2 to v4.1.3? Should I upgrade to 4.1.3?

    Could you please help me to solve this problem?

    Regards,
    Razi

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin
  • RaziRazi Member

    Hi @bhanuGandham

    I used the last version (GATK/4.1.3.0-Java-1.8) and the last version of vcf tools (VCFtools/0.1.16), but I still have the same error!!!! :(

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @Razi

    As mentioned in your above post, looks like the error you see is in the 4th step where you are using vcftools, am I right? Unfortunately we do not support vcftools. I did however see a post on www.biostars.org with a similar error as yours. Take a look: https://www.biostars.org/p/342914/
    If that doesn't help maybe post your question in www.seqanswers.com www.biostars.org www.stackoverflow.com or https://bioinformatics.stackexchange.com/

Sign In or Register to comment.