windows ubuntu 子系统，肿瘤全外篇，3. gatk中的BaseRecalibrator，HaplotypeCaller，ApplyVQSR

2中，我们对测序数据进行了比对，bam排序，标记重复和建立索引。这次我们就直接可以进入gatk流程了。

1.碱基质量校正。

  BaseRecalibrator 是 GATK（Genome Analysis Toolkit）工具集中的一个工具，用于校正测序数据中碱基质量的偏差，从而提高变异检测的准确性。

        在测序数据分析中，由于测序平台、化学试剂、PCR等因素的影响，碱基质量可能存在偏差，例如一些碱基可能比其他碱基更容易出错。BaseRecalibrator 通过分析测序数据中的已知变异位点，统计每个碱基的观察频率和预期频率，然后计算出修正因子，以校正碱基质量的偏差。

  BaseRecalibrator 的工作流程大致如下：

  1. 收集变异位点信息： 从已知的变异位点（如 dbSNP 数据库中的位点）中提取出来，并根据这些位点的信息统计每个碱基的观察频率和质量分数。
  2. 建立模型： 基于收集到的信息，建立一个模型，用于描述碱基质量与观察频率之间的关系。
  3. 计算校正因子： 使用建立的模型，计算每个碱基的校正因子，用于修正碱基质量的偏差。
  4. 生成校正后的 BAM 文件： 应用计算得到的校正因子，对测序数据中的碱基质量进行校正，生成校正后的 BAM 文件。

1.1生成校正表

gatk BaseRecalibrator \
        -R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
        -I ./L1.markdup.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 L1.IndelRealigner.intervals

gatk BaseRecalibrator: 这是运行 GATK 工具集中 BaseRecalibrator 工具的命令。

-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta: -R 参数指定参考基因组的 FASTA 文件路径。这里使用了 Homo sapiens 的 hg38 版本作为参考基因组。

-I ./L1.markdup.bam: -I 参数指定输入的 BAM 文件，其中 L1.markdup.bam 是标记去重后的测序数据文件。

-known-sites /mnt/h/db/gatk/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz: -known-sites 参数指定已知的变异位点文件，这里使用了 1000 Genomes 项目中高置信度的单核苷酸多态性（SNPs）的 VCF 文件。

-known-sites /mnt/h/db/gatk/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz: -known-sites 参数再次指定已知的变异位点文件，这次使用了 Mills & 1000 Genomes 项目中的金标准插入/缺失（Indels）的 VCF 文件。

-O L1.IndelRealigner.intervals: -O 参数指定输出文件的路径和文件名，这里生成的文件名为 L1.IndelRealigner.intervals，用于存储间插片段重对齐所需的间隔信息。

1.2 应用校正表

gatk ApplyBQSR \
        -R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
        -I ./L1.markdup.bam \
        -bqsr L1.IndelRealigner.intervals \
        -O L1_recalibrated_reads.bam

gatk ApplyBQSR: 这是运行 GATK 工具集中 ApplyBQSR 工具的命令。

-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta: -R 参数指定参考基因组的 FASTA 文件路径，这里使用了 Homo sapiens 的 hg38 版本作为参考基因组。

-I ./L1.markdup.bam: -I 参数指定输入的 BAM 文件，其中 949743-T_L2_1.markdup.bam 是标记去重后的测序数据文件。

-bqsr L1.IndelRealigner.intervals: -bqsr 参数指定之前生成的校正因子文件，这里使用了L1.IndelRealigner.intervals 文件。

-O L1_recalibrated_reads.bam: -O 参数指定输出文件的路径和文件名，这里生成的文件名为 L1_recalibrated_reads.bam，用于存储校正后的测序数据。

2.单样本变异调用

HaplotypeCaller是GATK（Genome Analysis Toolkit）中的一个重要工具，用于对测序数据进行变异（SNP和Indel）的检测和分析。它采用了全局的分析方法，利用多个样本的信息来识别可能的变异位点，并生成高质量的变异呼叫。用这个方法不用提前局部重比对。
主要特点：

  1. 基于组装的变异检测：HaplotypeCaller通过对碱基片段进行组装来识别变异位点，而不是简单地对每个碱基位置进行独立的分析。这种方法可以更准确地区分变异和测序错误。

  2. 局部重比对：在变异检测之前，HaplotypeCaller会进行局部重比对，以更好地处理插入缺失（Indels）。

  3. 多样本模式：HaplotypeCaller能够同时处理多个样本的测序数据，以提高变异检测的准确性。

  4. 基于贝叶斯方法的变异呼叫：HaplotypeCaller使用贝叶斯方法来计算每个可能的变异的概率，并在确定最终呼叫时考虑这些概率。

  5. 局部重比对：首先，HaplotypeCaller会对比对过的测序数据进行局部重比对，以纠正插入缺失导致的比对错误，并提高变异检测的准确性。

  6. 组装候选片段：接下来，HaplotypeCaller会从测序数据中组装出可能的等位基因片段（haplotypes），这些片段是由候选变异引起的，或者是由于测序错误而引入的。

  7. 变异检测：然后，对每个组装出的候选片段进行变异检测，通过比较片段与参考基因组之间的差异来识别SNP和Indel。

  8. 变异呼叫：最后，使用贝叶斯方法计算每个候选变异的概率，并基于这些概率进行变异呼叫，以确定最终的变异集合。

gatk  HaplotypeCaller \
     -R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta  \
     -I L1_recalibrated_reads.bam \
      --dbsnp /mnt/h/db/gatk/hg38/dbsnp_146.hg38.vcf.gz \
      -O L1_raw.vcf \
      1>L1_log.HC 2>&1

-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta: -R 参数指定参考基因组的 FASTA 文件路径，这里使用了 Homo sapiens 的 hg38 版本作为参考基因组。

-I L1_recalibrated_reads.bam: -I 参数指定输入的 BAM 文件，其中 L1_recalibrated_reads.bam 是经过校正后的测序数据文件。

--dbsnp /mnt/h/db/gatk/hg38/dbsnp_146.hg38.vcf.gz: --dbsnp 参数指定已知的变异位点文件，这里使用了 dbSNP 数据库中的 hg38 版本的 VCF 文件。

-O L1_raw.vcf: -O 参数指定输出文件的路径和文件名，这里生成的文件名为 L1_raw.vcf，用于存储未经过过滤的原始变异结果。

1>L1_log.HC 2>&1: 这部分是将标准输出和标准错误输出重定向到文件 L1_log.HC 中，以便记录日志信息。

        GATK HaplotypeCaller之后，需要进行质控，可以用GATK Variant Quality Score Recalibration（VQSR）。VQSR的原理是利用自身数据和已知变异位点集之间的重叠，通过高斯混合模型（GMM）构建一个分类器来对变异数据进行打分，从而评估每个位点的可靠性。具体而言，VQSR使用已知变异位点集作为训练集，通过统计特征（如变异质量、测序深度等）和已知位点的质量得分，构建一个GMM模型。然后，该模型会根据自身数据的特征和已知位点的质量得分，为每个变异位点分配一个新的质量得分。最终，VQSR根据这个新的质量得分，对变异位点进行过滤和分类，以确定哪些位点是高可靠度的变异。通过这一方法，VQSR可以帮助识别那些在自身数据中可能存在假阳性或假阴性的变异位点，并提供更可靠的变异集合，从而提高变异检测的准确性和可信度。

3.变异质控

接下来的命令需要 R: conda insatll R

3.1生成校正模型

gatk VariantRecalibrator \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-V L1_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/L1_WES.snp.recal.vcf \
--tranches-file ./snprecal/L1_WES.snp.tranches \
--rscript-file ./snprecal/L_WES.snp.plots.R

-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta：指定参考基因组的FASTA文件路径。

-V L1_raw.vcf：指定原始的VCF文件路径，这是待校正的文件。

-resource:hapmap,known=false,training=true,truth=true,prior=15.0 /mnt/h/db/gatk/hg38/hapmap_3.3.hg38.vcf.gz：指定用于训练和校正的参考资源。这里包括hapmap、omini、1000G和dbsnp四个资源。每个资源都有不同的权重（prior），hapmap和omini的truth值为true，而1000G和dbsnp的known值为true。这些资源通常用于帮助确定变异的真实性和可靠性。

-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum：指定用于校正的注释标签（annotations）。这些标签包括深度（DP）、质量/深度比值（QD）、分段比例（FS）、strand偏差（SOR）、读取位置秩和（ReadPosRankSum）以及映射质量秩和（MQRankSum）等。

-mode SNP：指定模式为SNP，表示这个命令用于对SNP进行校正。

-tranche：指定用于切分数据的阈值。这里有四个不同的阈值：100.0、99.9、99.0和95.0。这些阈值用于确定变异的可靠性等级，例如，阈值100.0表示最高可靠性。

-O ./snprecal/L1_WES.snp.recal.vcf：指定输出的校正后的VCF文件路径。

--tranches-file ./snprecal/L1_WES.snp.tranches：指定输出的阈值文件路径，其中包含了根据给定阈值进行的变异分级。

--rscript-file ./snprecal/L_WES.snp.plots.R：指定输出的R脚本文件路径，用于绘制校正过程中生成的各种图表。

        这个命令通过对原始的VCF文件进行训练和校正，生成了一个校正后的VCF文件，其中包含了经过过滤和分级的SNP数据，以及相应的阈值文件和R脚本用于可视化分析。

3.2 运用校正模型

gatk ApplyVQSR \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-V ../L_1_raw.vcf \
--ts-filter-level 99.0 \
--tranches-file L1_WES.snp.tranches \
--recal-file L1_WES.snp.recal.vcf \
-mode SNP \
-O 19P0126636WES.snps.VQSR.vcf.gz

-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta：指定参考基因组的FASTA文件路径。

-V ../L1_raw.vcf：指定原始的VCF文件路径，这是待应用VQSR的文件。

--ts-filter-level 99.0：指定过滤阈值的水平，这里设定为99.0，表示将满足99.0%可靠性阈值的变异保留下来。

--tranches-file L1_WES.snp.tranches：指定之前生成的阈值文件路径，用于对变异进行过滤。

--recal-file L1_WES.snp.recal.vcf：指定之前生成的校正文件路径，用于对变异进行校正。

-mode SNP：指定模式为SNP，表示这个命令用于对SNP数据进行VQSR。

-O 19P0126636WES.snps.VQSR.vcf.gz：指定输出的VCF文件路径，这是应用VQSR后的结果文件，经过过滤和校正的SNP数据。

        这个命令将会根据之前建立的校正模型，对原始的VCF文件中的SNP进行过滤和校正，并生成一个经过VQSR处理的结果文件。

         这就生成了一个VCF文件，流程跑完了，但是感觉还是啥也不会，先不说了，还是有空补充下基础知识吧。跑流程的过程中遇到了无数大坑，一坑更比一坑坑，大家想学，还是要自己跑一下。

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.mzph.cn/pingmian/141.shtml

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈email:809451989@qq.com，一经查实，立即删除！

相关文章

02 - Git 之命令 + 删除 + 版本控制 + 分支 + 标签 + 忽略文件 + 版本号

02 - Git 之命令 + 删除 + 版本控制 + 分支 + 标签 + 忽略文件 + 版本号

1 Git相关概念 1.1 以下所谈三个区，文件并不只是简单地在三个区转移，而是以复制副本的方式转移 使用 Git 管理的项目，拥有三个区域，分别是 Working area工作区（亦称为 工作树Working Tree）、stage area …
阅读更多...
[渗透测试学习] Monitored-HackTheBox

[渗透测试学习] Monitored-HackTheBox

Monitored-HackTheBox 信息搜集 nmap扫描一下端口 nmap -sV -sC -v --min-rate 1000 10.10.11.248扫描结果如下 PORT STATE SERVICE VERSION 22/tcp open ssh OpenSSH 8.4p1 Debian 5+deb11u3 (protocol 2.0) 80/tcp open http Apache httpd …
阅读更多...
视觉位置识别与多模态导航规划

视觉位置识别与多模态导航规划

前言 机器人感知决策是机器人移动的前提，机器人需要对周围环境实现理解，而周围环境通常由静态环境与动态环境构成。机器人在初始状态或者重启时需要确定当前所处的位置，然后根据用户的指令或意图，开展相应移动或抓取操作。通过视觉…
阅读更多...
【分治】Leetcode 排序数组

【分治】Leetcode 排序数组

题目讲解 912. 排序数组 算法讲解 我们这里使用三指针&#xff0c;将数组分成三块&#xff1a;<key 和 key 和 >key,如果当前指针指向的数字<key&#xff0c;我们就swap(nums[left]), nums[i] 。如果当前的数字key &#xff0c;就让i。如果当前的数字>key&…
阅读更多...
geolife笔记/python笔记：trackintel.io.read_geolife

geolife笔记/python笔记：trackintel.io.read_geolife

此函数解析 geolife_path 目录中可用的所有 geolife 数据 trackintel.io.read_geolife(geolife_path, print_progressFalse) 参数&#xff1a; geolife_path (str) 包含 geolife 数据的目录路径 print_progress (Bool, 默认为 False)如果设置为 True&#xff0c;则显示每个…
阅读更多...
退出 beeline

退出 beeline

退出 beeline 的命令是 !quit或 !exit 或者&#xff0c;直接来 Ctrl-D 我们下期见&#xff0c;拜拜&#xff01;
阅读更多...
java自定义排序规则

java自定义排序规则

java自定义排序规则 背景&#xff1a;开发的一个查询功能&#xff0c;主要字段&#xff1a;账号、币种等&#xff0c;业务提供了一个显示顺序&#xff0c;没有什么规则的排序&#xff0c;例如 3,1,6,2 无规则的顺序 思路&#xff1a;通过提供的顺序来指定排序规则 代码&#…
阅读更多...
Linux 内核复合页(compound page)原理分析

Linux 内核复合页(compound page)原理分析

源码基于&#xff1a;Linux5.15 约定&#xff1a; 芯片架构&#xff1a;ARM64内存架构&#xff1a;UMACONFIG_ARM64_VA_BITS&#xff1a;39CONFIG_ARM64_PAGE_SHIFT&#xff1a;12CONFIG_PGTABLE_LEVELS &#xff1a;3 1. 简介 复合页(Compound Page) 只是将两个或更多物理上…
阅读更多...
【笔试强训】DFS、优先队列、滑动窗口笔试题目！

【笔试强训】DFS、优先队列、滑动窗口笔试题目！

文章目录 1. 单词搜索2. 除 2 操作3. dd 爱框框 1. 单词搜索 题目链接 解题思路&#xff1a; DFS (深度优先遍历)&#xff0c;用一个 pos 记录要匹配单词 word 的位置&#xff0c;每次与 pos 进行匹配判断&#xff08;这样做的好处是不用把答案存下来&#xff09; 注意细节…
阅读更多...
算法：期望场景；鲁棒优化

算法：期望场景；鲁棒优化

部分代码 for i1:T stst[D_DGk(i)*min_P_DG<P_DGk(i)<D_DGk(i)*max_P_DG]; end for i2:T indicatorD_DGk(i)-D_DGk(i-1); rangei:min(T,iT_up-1); st st[D_DGk(range)>indicator]; end for i2:T indicatorD_DGk(i-1)-D_DGk(i); rangei:min(T…
阅读更多...
时序分解 | Matlab实现TVF-EMD时变滤波器的经验模态分解信号分量可视化

时序分解 | Matlab实现TVF-EMD时变滤波器的经验模态分解信号分量可视化

时序分解 | Matlab实现TVF-EMD时变滤波器的经验模态分解信号分量可视化 目录 时序分解 | Matlab实现TVF-EMD时变滤波器的经验模态分解信号分量可视化效果一览基本介绍程序设计参考资料 效果一览 基本介绍 Matlab实现TVF-EMD(时变滤波器的经验模态分解)可直接替换 Matlab语言 1.…
阅读更多...
IDEA远程调试debug

IDEA远程调试debug

IDEA远程调试debug jar包启动脚本配置IDEA配置 通俗的说&#xff1a;本地有代码&#xff0c;服务器项目出现问题&#xff0c;环境的中间件配置不同&#xff0c;用idea远程调试&#xff0c;能快速定位问题&#xff0c;解决问题。 jar包启动脚本配置 jdk5-8写法 java -Xdebug -…
阅读更多...
xftp、xshell连不上虚拟机解决方法

xftp、xshell连不上虚拟机解决方法

一、检查连接虚拟机ip看是否正确 查看虚拟机系统 IP ifconfig 二、检查虚拟机防火墙是否关闭 查看防火墙状态(ubuntu) sudo ufw status 关闭防火墙 sudo ufw disable 查看防火墙状态(centos) systemctl status firewalld.service 关闭防火墙 systemctl stop firewalld.se…
阅读更多...
基于Python dlib的实时人脸识别，附源码

基于Python dlib的实时人脸识别，附源码

博主介绍&#xff1a;✌IT徐师兄、7年大厂程序员经历。全网粉丝15W、csdn博客专家、掘金/华为云//InfoQ等平台优质作者、专注于Java技术领域和毕业项目实战✌ &#x1f345;文末获取源码联系&#x1f345; &#x1f447;&#x1f3fb; 精彩专栏推荐订阅&#x1f447;&#x1f3…
阅读更多...
【STM32】西南交大嵌入式系统设计实验：环境配置

【STM32】西南交大嵌入式系统设计实验：环境配置

把走过的坑记录一下&#xff0c;希望后来人避坑 No ST-Link device detected.问题解决 如果跟着指导书出现这个问题&#xff1a; 直接跳过这一步不用再更新固件&#xff0c;后面直接创建项目写程序就行了。 在keil里配置成用DAP_link即可。 详细的可以看这篇文章&#xff1a…
阅读更多...
30. 【Android教程】吐司提示：Toast 的使用方法

30. 【Android教程】吐司提示：Toast 的使用方法

在使用 Android 手机的时候&#xff0c;有没有遇到过如图中这种类型的消息提示&#xff1f; 这个在 Android 中被称为 Toast&#xff0c;用来短暂的展示一些简短的提示信息。相比弹窗来讲它对用户的打扰更小&#xff0c;在提示一段时间之后会自动消失&#xff0c;通常用来提示当…
阅读更多...
【Python-基础】字符串合集

【Python-基础】字符串合集

字符串格式化 f # 例如: # f{train_path}/{f}: 将train_path字符串和f字符串结合 # f{root}.csv:将root字符串和.csv字符串结合判断字符串是否以…结尾 root.endswith(".csv") # True未待完续…
阅读更多...
java 溯本求源之基础（十七）之Monitoring--jstatd

java 溯本求源之基础（十七）之Monitoring--jstatd

目录 1.简介 2.jstatd 命令概述 2.1功能和工作原理 2.2 使用场景 2.3 运行条件 3.命令行选项解析 3.1常用选项 4. 实际应用示例 4.1 内部 RMI 注册表配置 4.1.1 基础示例&#xff1a; 4.1.2 外部 RMI 注册表配置 4.1.3启用 RMI 日志记录 5.总结 1.简介 jstatd 是 Jav…
阅读更多...
Abp中ef操作新增重复子级数据问题

Abp中ef操作新增重复子级数据问题

在偶然开发中&#xff0c;导入的多条数据中&#xff0c;可能都存在同一个字段生成主外键关联子级数据的逻辑&#xff0c;此时循环去生成子级数据&#xff0c;会导致重复添加子级数据&#xff0c;有点绕吧&#xff0c;那就那实例说吧 如下&#xff1a;实现导入两条论文数据&…
阅读更多...
【贪心算法】发饼干问题详解python

【贪心算法】发饼干问题详解python

作者介绍&#xff1a;10年大厂数据\经营分析经验&#xff0c;现任大厂数据部门负责人。 会一些的技术&#xff1a;数据分析、算法、SQL、大数据相关、python 欢迎加入社区&#xff1a;码上找工作 作者专栏每日更新&#xff1a; LeetCode解锁1000题: 打怪升级之旅 python数据分析…
阅读更多...
最新文章

Copyright @ 2022~2023