windows ubuntu子系统，4.HaplotypeCaller和Mutect2总结

2024-05-06 20:36
文章标签 ubuntu 总结 windows 子系统 haplotypecaller mutect2

本文主要是介绍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总结的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/965308。

相关文章

JavaSE正则表达式用法总结大全

JavaSE正则表达式用法总结大全

《JavaSE正则表达式用法总结大全》正则表达式就是由一些特定的字符组成,代表的是一个规则,：本文主要介绍JavaSE正则表达式用法的相关资料,文中通过代码介绍的非常详细,需要的朋友可以参考下... 目录常用的正则表达式匹配符正则表China编程达式常用的类Pattern类Matcher类PatternSynta
阅读更多...
基于Python开发Windows屏幕控制工具

基于Python开发Windows屏幕控制工具

《基于Python开发Windows屏幕控制工具》在数字化办公时代,屏幕管理已成为提升工作效率和保护眼睛健康的重要环节,本文将分享一个基于Python和PySide6开发的Windows屏幕控制工具,... 目录概述功能亮点界面展示实现步骤详解1. 环境准备2. 亮度控制模块3. 息屏功能实现4. 息屏时间
阅读更多...
在Windows上使用qemu安装ubuntu24.04服务器的详细指南

在Windows上使用qemu安装ubuntu24.04服务器的详细指南

《在Windows上使用qemu安装ubuntu24.04服务器的详细指南》本文介绍了在Windows上使用QEMU安装Ubuntu24.04的全流程：安装QEMU、准备ISO镜像、创建虚拟磁盘、配置... 目录1. 安装QEMU环境2. 准备Ubuntu 24.04镜像3. 启动QEMU安装Ubuntu4
阅读更多...
Windows下C++使用SQLitede的操作过程

Windows下C++使用SQLitede的操作过程

《Windows下C++使用SQLitede的操作过程》本文介绍了Windows下C++使用SQLite的安装配置、CppSQLite库封装优势、核心功能（如数据库连接、事务管理）、跨平台支持及性能优... 目录Windows下C++使用SQLite1、安装2、代码示例CppSQLite：C++轻松操作SQ
阅读更多...
基于Python实现一个Windows Tree命令工具

基于Python实现一个Windows Tree命令工具

《基于Python实现一个WindowsTree命令工具》今天想要在Windows平台的CMD命令终端窗口中使用像Linux下的tree命令,打印一下目录结构层级树,然而还真有tree命令,但是发现... 目录引言实现代码使用说明可用选项示例用法功能特点添加到环境变量方法一：创建批处理文件并添加到PATH1
阅读更多...
SQL中JOIN操作的条件使用总结与实践

SQL中JOIN操作的条件使用总结与实践

《SQL中JOIN操作的条件使用总结与实践》在SQL查询中,JOIN操作是多表关联的核心工具,本文将从原理,场景和最佳实践三个方面总结JOIN条件的使用规则,希望可以帮助开发者精准控制查询逻辑... 目录一、ON与WHERE的本质区别二、场景化条件使用规则三、最佳实践建议1.优先使用ON条件2.WHERE用
阅读更多...
Windows的CMD窗口如何查看并杀死nginx进程

Windows的CMD窗口如何查看并杀死nginx进程

《Windows的CMD窗口如何查看并杀死nginx进程》：本文主要介绍Windows的CMD窗口如何查看并杀死nginx进程问题,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地... 目录Windows的CMD窗口查看并杀死nginx进程开启nginx查看nginx进程停止nginx服务
阅读更多...
Ubuntu设置程序开机自启动的操作步骤

Ubuntu设置程序开机自启动的操作步骤

《Ubuntu设置程序开机自启动的操作步骤》在部署程序到边缘端时,我们总希望可以通电即启动我们写好的程序,本篇博客用以记录如何在ubuntu开机执行某条命令或者某个可执行程序,需要的朋友可以参考下... 目录1、概述2、图形界面设置3、设置为Systemd服务1、概述测试环境：Ubuntu22.04 带图
阅读更多...
Nginx Location映射规则总结归纳与最佳实践

Nginx Location映射规则总结归纳与最佳实践

《NginxLocation映射规则总结归纳与最佳实践》Nginx的location指令是配置请求路由的核心机制,其匹配规则直接影响请求的处理流程,下面给大家介绍NginxLocation映射规则... 目录一、Location匹配规则与优先级1. 匹配模式2. 优先级顺序3. 匹配示例二、Proxy_pa
阅读更多...
Android学习总结之Java和kotlin区别超详细分析

Android学习总结之Java和kotlin区别超详细分析

《Android学习总结之Java和kotlin区别超详细分析》Java和Kotlin都是用于Android开发的编程语言,它们各自具有独特的特点和优势,：本文主要介绍Android学习总结之Ja... 目录一、空安全机制真题 1：Kotlin 如何解决 Java 的 NullPointerExceptio
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2