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