本文主要是介绍windows ubuntu子系统，4.HaplotypeCaller和Mutect2总结，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
建立分析目录
#拷贝数据到分析目录
cp /mnt/h/data/nuodaodataRENCEL/zuoborequest/*fq.gz ./
#质量检测
ls *fq.gz | xargs fastqc -t 20 -o fastqc_row/
multiqc fastqc_row/
#制作含有样本名称的文本
ls *_1.fq.gz >1
ls *_2.fq.gz >2
paste 1 2 > merge
sed -i "s/.fq.gz//g" merge
#fastq.gz质控
cat merge | while read id ; do
fastp -i "${id}_1.fq.gz" -I "${id}_2.fq.gz" -o fastp/"${id}_clean_1.fq.gz" -O fastp/"${id}_clean_2.fq.gz" -j fastp/"${id}.json" -h fastp/"${id}.html"
done #成功
#bwa-mm
cat ../merge | while read id ; do
bwa mem -t 20 -R "@RG\tID:${id}\tLB:${id}\tPL:Illumina\tSM:${id}" /mnt/h/db/bwa.db/hg38.fa ./"${id}_clean_1.fq.gz" ./"${id}_clean_2.fq.gz" | samtools sort -@ 2 -o human/"${id}.bam"
done
#标记pcr重复
cat sample.txt | while read id ; do
java -jar /mnt/h/softwore/picard/picard.jar MarkDuplicates \
I="${id}.bam" \
O="${id}.markup.bam" \
M="${id}.markdup_metrics.txt"
done
#samtools构建索引
cat sample.txt | while read id ; do
samtools index "${id}.markup.bam" -@ 10
done
"""
#这段代码是用于进行GATK4中的BaseRecalibrator操作，用于进行碱基质量得分校正。
cat sample.txt | while read id ; do
gatk BaseRecalibrator \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-I "${id}.markup.bam" \
-known-sites /mnt/h/db/gatk/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
-known-sites /mnt/h/db/gatk/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
-O "${id}.IndelRealigner.intervals"
done#成功
# 第二步：应用校正表
cat sample.txt | while read id ; do
gatk ApplyBQSR \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-I "${id}.markup.bam" \
-bqsr ${id}.IndelRealigner.intervals \
-O "${id}.recalibrated_reads.bam"
done #成功
#单样本的变异调用
cat sample.txt | while read id ; do
gatk HaplotypeCaller \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-I "${id}.recalibrated_reads.bam" \
--dbsnp /mnt/h/db/gatk/hg38/dbsnp_146.hg38.vcf.gz \
-O "${id}.raw.vcf" \
--native-pair-hmm-threads 36 \
1>"${id}.log.HC" 2>&1
done
ATK的HaplotypeCaller命令对一个名为sample.txt的文件中的每个样本进行变异检测。
使用了以下参数：
-R：参考基因组文件的路径
-I：输入的校正后的BAM文件路径
--dbsnp：已知的单核苷酸多态性数据库（dbSNP）文件的路径
-O：输出的VCF文件路径
--native-pair-hmm-threads：使用的线程数
1>：标准输出日志的路径
2>&1：将标准错误重定向到标准输出
该命令从sample.txt文件中逐行读取样本ID，并对每个样本运行HaplotypeCaller命令。HaplotypeCaller用于发现样本中的变异位点，并生成一个原始的VCF文件。
#变异质量矫正
cat sample.txt | while read id ; do
gatk VariantRecalibrator \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-V "${id}.raw.vcf" \
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 /mnt/h/db/gatk/hg38/hapmap_3.3.hg38.vcf.gz \
-resource:omini,known=false,training=true,truth=false,prior=12.0 /mnt/h/db/gatk/hg38/1000G_omni2.5.hg38.vcf.gz \
-resource:1000G,known=false,training=true,truth=false,prior=10.0 /mnt/h/db/gatk/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 /mnt/h/db/gatk/hg38/dbsnp_146.hg38.vcf.gz \
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum -mode SNP -tranche 100.0 \
-tranche 99.9 -tranche 99.0 -tranche 95.0 -tranche 90.0 \
-O ./snprecal/"${id}.WES.snp.recal.vcf" \
--tranches-file ./snprecal/"${id}.WES.snp.tranches" \
--rscript-file ./snprecal/"${id}.WES.snp.plots.R"
done
cat sample.txt | while read id ; do ... done：这部分是一个循环，它从名为sample.txt的文件中读取每一行的内容，并将其赋给变量id。然后，循环体中的命令将针对每个id执行。
gatk VariantRecalibrator ...：这是GATK工具VariantRecalibrator的命令行调用。它用于对VCF文件进行校正和过滤，以提高SNP的质量和可信度。
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta：指定参考基因组的FASTA文件路径。
-V "${id}.raw.vcf"：指定输入VCF文件的路径，其中"${id}"将被循环中的当前id替换。
-resource:hapmap ... -resource:dbsnp ...：这些是用于校正的参考资源文件。它们提供了已知的SNP信息，帮助GATK评估和校正输入VCF文件中的SNP。
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum：指定用于评估SNP质量的注释器。这些是用于计算各种质量指标的工具。
-mode SNP：指定模式为SNP。
--tranche ...：指定不同阈值的过滤等级。
-O ./snprecal/"${id}.WES.snp.recal.vcf"：指定校正后的VCF文件的输出路径。
--tranches-file ... --rscript-file ...：指定输出的tranches文件和R脚本文件的路径。
#质控SNP
cat ../sample.txt | while read id ; do
gatk ApplyVQSR \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-V ../"${id}.raw.vcf" \
--ts-filter-level 99.0 \
--tranches-file "${id}.WES.snp.tranches" \
--recal-file "${id}.WES.snp.recal.vcf" \
-mode SNP \
-O "${id}.WES.snps.VQSR.vcf.gz"
done
#一定要注意换行符问题，路径问题。
#上个命令行的流程
使用GATK中的ApplyVQSR工具
将指定的参考基因组文件（/mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta）传递给-R参数
使用原始VCF文件（../"${id}.raw.vcf"）作为输入 (-V参数)
指定VQSR过滤的阈值为99.0% (--ts-filter-level参数)
指定tranches文件和校正文件分别为"𝑖𝑑.𝑊𝐸𝑆.𝑠𝑛𝑝.𝑡𝑟𝑎𝑛𝑐ℎ𝑒𝑠"和"id.WES.snp.tranches"和"{id}.WES.snp.recal.vcf" (--tranches-file和--recal-file参数)
指定处理的模式为SNP (-mode参数)
将处理后的VCF文件输出为"${id}.WES.snps.VQSR.vcf.gz" (-O参数)
wc 1100817-4T_L3.WES.VQSR.vcf.gz
这条命令是用来统计文件中的行数、单词数和字节数的。其中，wc 表示 word count（单词计数），1100817-4T_L3.WES.VQSR.vcf.gz 是要统计的文件名
####
Mutect2是GATK（Genome Analysis Toolkit）中的一个工具，用于检测样本之间的单核苷酸变异（SNVs）和小的插入/缺失突变（indels）。它采用了一种基于贝叶斯统计模型的方法来过滤假阳性变异，并提供了高灵敏度和高特异性的变异检测结果。
重跑bwa mem
cat merge | while read id ; do
bwa mem -t 20 -R "@RG\tID:${id}\tLB:${id}\tPL:Illumina\tSM:${id}" /mnt/h/db/bwa.db/hg38.fa ../"${id}_clean_1.fq.gz" ../"${id}_clean_2.fq.gz" | samtools sort -@ 10 -o ./"${id}.bam"
done
#最重要的是-R参数的，设计到后来的肿瘤和正常样本的分组。
#标记pcr重复
cat merge | while read id ; do
java -jar /mnt/h/softwore/picard/picard.jar MarkDuplicates \
I="${id}.bam" \
O="${id}.markup.bam" \
M="${id}.markdup_metrics.txt"
done
#建立索引
cat merge | while read id ; do
samtools index "${id}.markup.bam" -@ 20
done
#在GATK4中，局部重比对（Local Realignment）的步骤已经被整合到了BaseRecalibrator和ApplyBQSR工具中。（我搜到的，不知道对不对）
#碱基质量校准和局部重比对
cat merge | while read id ; do
gatk BaseRecalibrator \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-I "${id}.markup.bam" \
-known-sites /mnt/h/db/gatk/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
-known-sites /mnt/h/db/gatk/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
-O "${id}.IndelRealigner.intervals"
done
cat merge | while read id ; do
gatk ApplyBQSR \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-I "${id}.markup.bam" \
-bqsr "${id}.IndelRealigner.intervals" \
-O "${id}.recalibrated_reads.bam"
done
#召唤变异
gatk Mutect2 \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-I 1100817-8N_L2.recalibrated_reads.bam \
-I 1100817-4T_L3.recalibrated_reads.bam \
-tumor 1100817-4T_L3 \
-normal 1100817-8N_L2 \
--germline-resource /mnt/h/db/gatk/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
-O mutect2.vcf.gz \
-bamout mutect2.bam
FilterMutectCalls 是 GATK 工具的一部分，用于过滤 Mutect2 工具所生成的变异结果（VCF 文件）。以下是 FilterMutectCalls 工具的用法：
gatk FilterMutectCalls \
-V input.vcf.gz \
-O output.vcf.gz
其中，-V 参数指定输入的 VCF 文件路径，-O 参数指定输出的 VCF 文件路径。此外，还可以使用以下参数对变异结果进行过滤：
--contamination-table：污染表格路径，包含样本的污染比例信息。
--tumor-segmentation：肿瘤分割路径，包含突变在肿瘤中的分布信息。
--initial-tumor-lod：初始肿瘤 LOD 临界值，低于该值的变异将被过滤。
--tumor-lod：肿瘤 LOD 临界值，低于该值的变异将被过滤。
--normal-artifact-lod：正常样本伪造 LOD 临界值，高于该值的变异将被过滤。
--somatic-artifact-lod：体细胞伪造 LOD 临界值，高于该值的变异将被过滤。
--cluster-window-size：聚类窗口大小，对变异位置进行聚类。
--cluster-size：聚类大小，超出聚类大小的变异将被过滤。
--max-events-in-region：区域内的最大变异数量，超出该值的变异将被过滤
"""
#过滤变异
gatk FilterMutectCalls \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-V mutect2.vcf.gz \
-O mutect2_filtered.vcf.gz
#注释，使用snpEff
官网下载
java -jar snpEff.jar -h
java -jar snpEff.jar databases | grep -i human
可以看到与人类相关的数据库如下：
GRCh37.87：使用转录本的人类基因组GRCh37版本。
GRCh37.p13：使用RefSeq转录本的人类基因组GRCh37版本。
GRCh38.mane.0.95.ensembl：使用MANE转录本v0.95和Ensembl ID的人类基因组GRCh38版本。
GRCh38.mane.0.95.refseq：使用MANE转录本v0.95和RefSeq ID的人类基因组GRCh38版本。
GRCh38.mane.1.0.ensembl：使用MANE转录本v1.0和Ensembl ID的人类基因组GRCh38版本。
GRCh38.mane.1.0.refseq：使用MANE转录本v1.0和RefSeq ID的人类基因组GRCh38版本。
GRCh38.mane.1.2.ensembl：使用MANE转录本v1.2和Ensembl ID的人类基因组GRCh38版本。
GRCh38.mane.1.2.refseq：使用MANE转录本v1.2和RefSeq ID的人类基因组GRCh38版本。
GRCh38.p13：使用RefSeq转录本的人类基因组GRCh38版本。
GRCh38.p14：使用RefSeq转录本的人类基因组GRCh38版本。
java -jar snpEff.jar databases
java -jar /mnt/h/softwore/snpeff/snpEff/snpEff.jar download -v GRCh38.p14 #下载数据库
SnpEff是一个用于注释遗传变异的工具，以下是使用SnpEff进行变异注释的命令代码示例：
java -jar snpEff.jar -v -c snpEff.config -i vcf -o vcf GRCh37.75 input.vcf > output.ann.vcf
这个命令中的参数解释如下：
-jar snpEff.jar：指定运行SnpEff的JAR文件。
-v：启用详细输出模式，会打印出注释的详细信息。
-c snpEff.config：指定SnpEff的配置文件，其中包含数据库和其他设置的信息。
-i vcf：指定输入文件的格式为VCF。
-o vcf：指定输出文件的格式为VCF。
GRCh37.75：指定要使用的基因组版本和数据库。
input.vcf：指定要注释的输入文件。
> output.ann.vcf：将注释后的结果重定向到名为output.ann.vcf的输出文件中。
#注释
java -jar /mnt/h/softwore/snpeff/snpEff/snpEff.jar \
-v -c snpEff.config -i vcf -o vcf GRCh38.p14 mutect2_filtered.vcf.gz > human.ann/mutect2_filtered.ann.vcf
这篇关于windows ubuntu子系统，4.HaplotypeCaller和Mutect2总结的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
