|
发表于 2023-6-22 16:21:35
|
查看: 5471 |
回复: 1
1.变异检测
- #生成gvcf
- time gatk HaplotypeCaller --emit-ref-confidence GVCF -R /share/home/xiehs/data/GATK/hg38/Homo_sapiens_assembly38.fasta -I merge.sorted.markdup.BQSR.bam -O merge.HC.g.vcf.gz
- #合并gvcf
- time gatk GenotypeGVCFs -R /share/home/xiehs/data/GATK/hg38/Homo_sapiens_assembly38.fasta -V merge.HC.g.vcf.gz -O merge.HC.vcf.gz
复制代码
2.结果过滤
2.1 VQSR
准备的已知变异集作为训练集,可以是 Hapmap、OMNI,1000G,dbsnp,瓶中基因组计划等这些国际性项目的数据,然后利用训练集对每一个位点进行过滤。利用 VariantRecalibrator工具进行机器学习,ApplyVQSR 工具进行处理。 VQSR 过滤 SNP 和 InDel 分别进行,首先处理 SNP,得到结果后,在进行 InDel 处理。
- #处理SNP
- # --max-gaussians默认值为8,报错提示需要降低
- gatk VariantRecalibrator --max-gaussians 6 -R /share/home/xiehs/data/GATK/hg38/Homo_sapiens_assembly38.fasta -V merge.HC.vcf.gz --resource:hapmap,known=false,training=true,truth=true,prior=15.0 /share/home/xiehs/data/GATK/hg38/hapmap_3.3.hg38.vcf.gz --resource:omni,known=false,training=true,truth=false,prior=12.0 /share/home/xiehs/data/GATK/hg38/1000G_omni2.5.hg38.vcf.gz --resource:1000G,known=false,training=true,truth=false,prior=10.0 /share/home/xiehs/data/GATK/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 /share/home/xiehs/data/GATK/hg38/dbsnp_146.hg38.vcf.gz -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -mode SNP -O merge.HC.snps.recal --tranches-file output.tranches --rscript-file output.plots.R
-
- gatk ApplyVQSR -R /share/home/xiehs/data/GATK/hg38/Homo_sapiens_assembly38.fasta -V merge.HC.vcf.gz -O merge.HC.snps.VQSR.vcf.gz --recal-file merge.HC.snps.recal --tranches-file merge.HC.snps.tranches -mode SNP
复制代码 1、HapMap,它来自国际人类单倍体型图计划,数据集包含了大量家系数据,并且有非常严格的质控和严密的实验验证,因此它的准确性是目前公认最高的。
2、Omni,这个数据源自 Illumina 的 Omni 基因型芯片,它的验证结果常常作为基因型的金标准。
3、1000G 千人基因组计划(1000 genomes project)质控后的变异数据,质控后,它包含的绝大部分都是真实的变异,但由于没办法做全面的实验验证,并不能排除含有少部分假阳的结果。
4、dbSNP。dbSNP 收集的数据,实际都是研究者们发表了相关文章提交上来的变异,这些变异很多是没做过严格验证的
处理 InDel
- #处理InDel
- gatk VariantRecalibrator -R /share/home/xiehs/data/GATK/hg38/Homo_sapiens_assembly38.fasta -V merge.HC.snps.VQSR.vcf.gz --max-gaussians 4 --resource:mills,known=false,training=true,truth=true,prior=12.0 /share/home/xiehs/data/GATK/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 /share/home/xiehs/data/GATK/hg38/dbsnp_146.hg38.vcf.gz -an QD -an DP -an FS -an SOR -an ReadPosRankSum -an MQRankSum -mode INDEL -O merge.HC.snps.indel.recal --tranches-file merge.HC.snps.indel.tranches --rscript-file merge.HC.snps.indel.plots.R
- gatk ApplyVQSR -R /share/home/xiehs/data/GATK/hg38/Homo_sapiens_assembly38.fasta -V merge.HC.snps.VQSR.vcf.gz -O merge.HC.snps.indel.VQSR.vcf.gz --truth-sensitivity-filter-level 99.0 --tranches-file merge.HC.snps.indel.tranches --recal-file merge.HC.snps.indel.recal -mode INDEL
复制代码
2.2 Hard-filter
Hard-filter 硬过滤,可以根据以下标准来进行过滤,gatk 过滤的时候,snp 与 indel 是分别进行的。也可以选择 bcftools 进行简单过滤。
QualByDepth(QD)
FisherStrand (FS)
StrandOddsRatio (SOR)
RMSMappingQuality (MQ)
MappingQualityRankSumTest (MQRankSum)
ReadPosRankSumTest (ReadPosRankSum)
此脚本为硬过滤(hard-filter)的方法,主要用于不能进行 VQSR 的情况,例如非人物种,或外显子,芯片数据等。
- # 使用SelectVariants,选出SNP
- time gatk SelectVariants -select-type SNP -V merge.HC.vcf.gz -O merge.HC.vcf.snp.gz
-
- # 为SNP作硬过滤
- time gatk VariantFiltration -V merge.HC.vcf.snp.gz --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "PASS" -O merge.HC.vcf.snp.filter.gz
-
- # 使用SelectVariants,选出Indel
- time gatk SelectVariants -select-type INDEL -V merge.HC.vcf.gz -O merge.HC.indel.vcf.gz
-
- # 为Indel作过滤
- time gatk VariantFiltration -V merge.HC.indel.vcf.gz --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "PASS" -O merge.HC.indel.filter.vcf.gz
-
- # 重新合并过滤后的SNP和Indel
- time gatk MergeVcfs -I merge.HC.snp.filter.vcf.gz -I merge.HC.indel.filter.vcf.gz -O merge.HC.filter.vcf.gz
复制代码
|
|