Добавление в файл GVCF экзомов 1000G, чтобы gatk VariantRecalibrator работал с небольшим образцом

У меня есть данные секвенирования небольшого ампликона размером 500 п.н. из нескольких образцов. Лучшие принципы GATK предполагают запуск VariantRecalibrator с файлами GVCF, которые я генерирую. Я пытаюсь заставить это работать, но получаю сообщение об ошибке «Найдены аннотации с нулевой дисперсией». Чтение руководства gatk и другие должности (например , 1, ЭГ2), я считаю , что это происходит потому , что у меня нет достаточного количества образцов в моем файле VCF. Я не знаю, как «дополнить» мой файл vcf данными от 1000G или exac! Любая помощь очень ценится.

Вот мой конвейер, чтобы добраться до этого момента, и ошибка, которую я вижу внизу.

      #====================================================
# 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.

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

0 ответов

Другие вопросы по тегам