0

我从几个样本中获得了一个 500 bp 的小扩增子的测序数据。GATK 最佳原则建议在我生成的 GVCF 文件上运行VariantRecalibrator。我试图让这个工作，但我收到一个关于“找到零方差的注释”的错误。阅读 gatk 手册和其他帖子（ eg1 ， eg2 ），我相信这是因为我的 vcf 文件中没有足够的样本。我不知道如何用 1000G 或 exac 的数据“填充”我的 vcf 文件！非常感谢任何帮助。

这是我到达这一点的管道，以及我在底部看到的错误。

#====================================================
# Alignment with BWA mem
#====================================================
bwa mem -M -t 1 -R $(echo "@RG\tID:$ID\tSM:$ID"_"$SM\tLB:$ID"_"$SM\tPL:ILLUMINA") "$REF" "$R1" "$R2" > "$OUT"/results/alignment/${SN}_bwa.sam

samtools view -Sb "$OUT"/results/alignment/${SN}_bwa.sam > "$OUT"/results/alignment/${SN}_bwa.bam

#====================================================
# Sort and mark duplicated with picard
#====================================================

picard SortSam I="$OUT"/results/alignment/${SN}_bwa.bam O="$OUT"/results/alignment/${SN}_bwa_sorted.bam SORT_ORDER=coordinate USE_JDK_DEFLATER=true USE_JDK_INFLATER=true

picard MarkDuplicates I="$OUT"/results/alignment/${SN}_bwa_sorted.bam O="$OUT"/results/alignment/${SN}_bwa_sorted_marked.bam M="$OUT"/results/alignment/marked_dup_metrics.txt USE_JDK_DEFLATER=true USE_JDK_INFLATER=true


#====================================================
# Run base-recalibrator and applyBSQR to recalibrate base qualities
#====================================================

# Base calibration
gatk --java-options "-Xmx4g" BaseRecalibrator -I "$OUT"/results/alignment/${SN}_bwa_sorted_marked.bam -R "$REF" --known-sites "$knownINDELS" -O "$OUT"/results/alignment/recal_data.table

gatk --java-options "-Xmx4g" AnalyzeCovariates -bqsr "$OUT"/results/alignment/recal_data.table -plots "$OUT"/results/alignment/AnalyzeCovariates.pdf

gatk --java-options "-Xmx4g" ApplyBQSR -R "$REF" -I "$OUT"/results/alignment/${SN}_bwa_sorted_marked.bam --bqsr-recal-file "$OUT"/results/alignment/recal_data.table -O "$OUT"/results/alignment/${SN}_bwa_sorted_marked_recalibrated.bam

#====================================================
# Index the bam file
#====================================================

samtools index "$OUT"/results/alignment/${SN}_bwa_sorted_marked_recalibrated.bam


#====================================================
# Run HaplotypeCaller
#====================================================

gatk --java-options "-Xmx4g" HaplotypeCaller \
      --intervals "$INTERVALS" \
      -R "$REF" \
      -I "$OUT"/results/alignment/${SN}_bwa_sorted_marked_recalibrated.bam \
      -O "$OUT"/results/variants/${SN}_g.vcf.gz \
      -ERC GVCF




#====================================================
# Genotype GVCFs
#====================================================

gatk --java-options "-Xmx4g" GenotypeGVCFs  --allow-old-rms-mapping-quality-annotation-data  -R "/Volumes/Seagate Expansion Drive/refs/hg38/gatk download/Homo_sapiens_assembly38.fasta" -V "$OUT"/results/variants/cohort.g.vcf.gz -O "$OUT"/results/variants/cohort.vcf.gz



#====================================================
# Normalise and index the vcf
#====================================================

# Normalise indels
bcftools norm "$OUT"/results/variants/cohort.vcf.gz "/Volumes/Seagate Expansion Drive/refs/hg38/gatk download/Homo_sapiens_assembly38.fasta" -m -both -o "$OUT"/results/variants/cohort.norm.vcf.gz -O z

# Index the vcf
tabix -p vcf "$OUT"/results/variants/cohort.norm.vcf.gz


#====================================================
# Variant recalibration
#====================================================

gatk VariantRecalibrator \
   -R "/Volumes/Seagate Expansion Drive/refs/hg38/gatk download/Homo_sapiens_assembly38.fasta" \
   -V "$OUT"/results/variants/cohort.norm.vcf.gz \
   -AS \
   --resource:hapmap,known=false,training=true,truth=true,prior=15.0 "/Volumes/Seagate Expansion Drive/refs/hg38/gatk download/hapmap_3.3.hg38.vcf.gz" \
   --resource:omni,known=false,training=true,truth=false,prior=12.0 "/Volumes/Seagate Expansion Drive/refs/hg38/gatk download/1000G_omni2.5.hg38.vcf.gz" \
   --resource:1000G,known=false,training=true,truth=false,prior=10.0 "/Volumes/Seagate Expansion Drive/refs/hg38/gatk download/1000G_phase1.snps.high_confidence.hg38.vcf.gz" \
   --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 "/Volumes/Seagate Expansion Drive/refs/hg38/gatk download/Homo_sapiens_assembly38.dbsnp138.vcf" \
   -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
   -mode SNP \
   -O "$OUT"/results/variants/output.AS.recal \
   --tranches-file "$OUT"/results/variants/output.AS.tranches \
   --rscript-file "$OUT"/results/variants/output.plots.AS.R

这是我看到的错误：

***********************************************************************

A USER ERROR has occurred: Bad input: Found annotations with zero variance. They must be excluded before proceeding.

***********************************************************************
4

0 回答 0