满分室间质评之GATK Somatic SNV+Indels+CNV+SV（上）

最新推荐文章于 2024-06-09 09:58:28 发布
最新推荐文章于 2024-06-09 09:58:28 发布
阅读量3.3k 收藏 10
点赞数 6
分类专栏： 生信分析流程 生信分析
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_39900139/article/details/107711305
版权

卫计委在2017年，2019年，2020年（还没有答案）提供标准数据用于肿瘤生信分析的室间质评。这样预知结果的数据自然是不能放过了，本文尝试参考GATK Best Practice：Somatic SNVs + Indels ，Cnvkit，Manta的pipeline来完成满分流程分析，也可以使用标准数据反向判断GATK Mutect2的实际准确度,算法优劣。

注：本文仅用于学习，距离真正的临床应用还有相当大距离，欢迎大佬批评指正

**1. 分析流程概览如下：
在这里插入图片描述

2. 本文用到的分析系统及分析流程文件

名称 （点击下载）备注
SliverWorkspace 2.0.277363提供运行控制平台/个人版
FFPE SNV CNV SV V1.workflow分析流程文件，可以一键导入分析平台(点击查看操作)
当然可以参照图片中运行脚本，shell里运行，效果也是一样
最终结果过滤脚本（python2.7 )及编译版本Illumina_pt2.bed 等用到的bed，intelval等文件
SnvAnnotationFilter.py SNV过滤脚本
CnvAnnotationFilter.py CNV过滤脚本
SvAnnotationFilter.py SV 过滤脚本
QcProcessor.py 获取整体QC数据的脚本
report_template.docx 分析报告模板
分析结果（pipeline结果与标准答案）result.zip pipeline结果与标准答案

3. 本文用到的原始文件，用fastqc查看质量状态是clean data，Q值均高于30，这里就不需要去接头和QC了。

百度网盘：链接: https://pan.baidu.com/s/1t9R14aQNP6Xk4tq1wFWJpg 提取码: u5yp

文件名（按照需要有调整）文件大小MD5
B1701_R1.fq.gz4.85G07d3cdccee41dbb3adf5d2e04ab28e5b
B1701_R2.fq.gz4.77Gc2aa4a8ab784c77423e821b9f7fb00a7
B1701NC_R1.fq.gz3.04G4fc21ad05f9ca8dc93d2749b8369891b
B1701NC_R2.fq.gz3.11Gbc64784f2591a27ceede1727136888b9

4. 本文用到的环境变量（目录/程序/文件/数值/字符）reference文件和数据库体积过于庞大请自行下载安装（如:ftp.broadinstitute.org/Annovar等等）

名称数值类型
sn1701字符
data/opt/data目录
result/opt/result目录
tools.java/opt/jdk1.8.0_162/bin/java程序
tools.bgzip/opt/tabix-0.2.6/bgzip程序
tools.bwa/opt/bwa/bwa程序
tools.cnvkit/usr/local/bin/cnvkit.py程序
tools.fastqc/opt/FastQC/fastqc程序
tools.gatk/opt/ref/gatk-4.1.3.0/gatk程序
tools.manta/opt/manta/dist/bin/configManta.py程序
tools.samtools/opt/samtools/samtools程序
tools.tabix/opt/tabix-0.2.6/tabix程序
tools.annovar/opt/ref/annovar/table_annovar.pl程序
tools.convertor/opt/ref/annovar/convert2annovar.pl程序
cutoff.cnv_max0.5数值
cutoff.cnv_min-0.5数值
cutoff.cnvdep1000数值
cutoff.depth500数值
cutoff.event2数值
cutoff.nvaf0.005数值
cutoff.sv200数值
cutoff.TLOD16数值
cutoff.vaf0.01数值
envis.scatter8数值
envis.threads32数值
envis.memory32G字符
refs.1000G/opt/ref/1000G_phase1.indels.hg19.vcf文件
refs.af_only/opt/ref/af-only-gnomad.raw.sites.hg19.vcf.gz文件
refs.bed/opt/ref/projects/Illumina_pt2.bed文件
refs.interval/opt/ref/projects/Illumina_pt2.interval_list文件
refs.dbsnp/opt/ref/dbsnp_138.hg19.vcf文件
refs.dict/opt/ref/ucsc.hg19.dict文件
refs.gene/opt/ref/hg19_refGene.txt文件
refs.hum/opt/ref/ucsc.hg19.fa文件
refs.mills/opt/ref/Mills_and_1000G_gold_standard.indels.hg19.vcf文件
refs.small.exac/opt/ref/small_exac_common_3_b37.vcf文件

5. 详细流程分析具体如下（不习惯图形软件的使用shell变量+脚本运行也是一样的）

GATK提供的标准流程从Normal/Tumor的fastq文件开始 map>reorder>sort>mark duplicate>recalibrator>apply BQSR

  • 流程输入文件

起始输入文件，使用变量${sn}表示样本编号

分析流程输入文件，这里使用变量${sn}表示样本编号，对室间质评文件名做了调整。

  • Tumor 比对，管道操作给samtools，直接输出bam格式文件。
    在这里插入图片描述

  • Reorder Bam
    GATK文档中这样的

Not to be confused with SortSam which sorts a SAM or BAM file with a valid sequence dictionary, ReorderSam reorders reads in a SAM/BAM file to match the contig ordering in a provided reference file, as determined by exact name matching of contigs. Reads mapped to contigs absent in the new reference are dropped. Runs substantially faster if the input is an indexed BAM file.

这一步省略其实对最终结果影响不大。
在这里插入图片描述

  • Sort

  • MarkDuplicate 标记重复
    在这里插入图片描述

  • 重新校正碱基质量值第一步，BaseRecalibrator：计算所有需要重校正的reads和特征值，然后把这些信息输出为校准表文件
    在这里插入图片描述

  • 重新校正碱基质量值第二步，ApplyBQSR：用第一步得到的校准表文件，重新调整BAM文件中的碱基质量值，并使用这个新的质量值重新输出一个新的BAM文件。

在这里插入图片描述

  • Normal 数据对着Tumor的步骤重复一遍。序列比对
    在这里插入图片描述

  • Normal Reorder
    在这里插入图片描述

  • Normal Sort
    在这里插入图片描述

  • Normal MarkDuplicate

在这里插入图片描述

  • Normal BaseRecalibrator

  • Normal ApplyBQSR
    在这里插入图片描述

  • Tumor GetPileupSummaries
    在这里插入图片描述

  • Normal GetPileupSummaries
    在这里插入图片描述

  • CalculateContamination
    在这里插入图片描述

  • Mutect2
    在这里插入图片描述

  • FilterMutectCalls 使用GATK提供的过滤器，过滤得到的突变
    在这里插入图片描述

  • 将过滤后的文件转换为Annovar注释所需要的格式
    在这里插入图片描述

  • 使用Annovar对结果注释
    在这里插入图片描述

table_annovar.pl ${result}/${sn}.annovar \
  /opt/ref/annovar/humandb \
  -buildver hg19 \
  -out ${result}/${sn} \
  -remove \
  -protocol refGene,avsnp150,esp6500siv2_all,1000g2015aug_all,clinvar_20190305,cosmic89_coding \
  -operation g,f,f,f,f,f \
  -nastring NAN

  • 使用自己写的脚本对注释后的结果过滤，比如按照室间质评要求，过滤掉突变频率低于1%，测序深度低于500的突变。对GATK某些过滤器过滤掉的结果进行保留和排除，后面使用IGV进行人工筛选。最终输出的结果为，${sn}.result.SNV.xls（其实是个csv文件，扩展名改为.xls是便于使用excel打开，很多人都这么干）
    在这里插入图片描述

  • 使用Manta Call SV，第一步，生成分析脚本文件runWorkflow.py
    在这里插入图片描述

  • 运行runWorkflow.py获取SV突变结果
    在这里插入图片描述

  • 使用自己过滤的脚本文件，对Manta获取的SV过滤。如根据SOMATICSCORE分数过滤，根据hg19_refGene.txt提供文件，计算突变基因等等。
    在这里插入图片描述

  • 使用CnvKIt，获取CNV突变
    在这里插入图片描述
    在这里插入图片描述

  • 使用CnvKit画图
    在这里插入图片描述

  • 使用py脚本文件，对CnvKit输出结果过滤。同样根据hg19_refGene.txt文件匹配基因，以及发生拷贝数变异的区域的外显子区域等。
    在这里插入图片描述

  • 获取整体QC数据，insertSize，depth等

  • Samtools获取碱基测序深度
    在这里插入图片描述
  • Samtools flagstat统计比对信息
    在这里插入图片描述
  • 使用py脚本文件对以上获取的数据做统一处理，获取整体QC状态结果。
    在这里插入图片描述

最终结果：
在这里插入图片描述

最终结果及各个命令/任务运行时间如下：
在这里插入图片描述

整个分析过程耗时 3h 58min 29s，比较耗时，这个版本为了逻辑清晰一些，多数任务都是串行运行，没有很好利用计算资源。

后续文章会对整个流程进行优化（主要是并行化），最终分析时间在1h 13min左右（提升约3倍）。

GATK 输出结果中SNV&INDEL的准确度问题，经过反复试验，不论如何设置过滤参数，最终的结果始终会有假阴性问题，这是GATK（4.1.3.0）中个别过滤器的问题，目前的补救措施是将部分GATK过滤器过滤掉的结果仍然包含在最终结果中，再使用IGV工具人工过滤（官方文档也是这么推荐的），判断该结果是否可靠。

weixin_39900139
关注
  • 6
    点赞
  • 踩
  • 10
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
专栏目录
09 拷贝数变异分析（GATK流程）
qq_23435961的博客
01-08 4708
09 拷贝数变异分析（GATK流程） 我们已经分析了 Somatic mutations，并进行了注释和可视化，接下来我们进行拷贝数变异的分析。 这里我们还是先从 GATK 的 somatic cnv 的最佳实践开始 拷 贝数变异（copy number variations, CNVs）是属于基因组结构变异（structural variation, SV），是指 DNA 片 段长度在 1Kb-3Mb 的基因组结构变异。我们首先从 GATK 的 CNV 流程开始 CNV 的分析。 首先是流程图，先从 b
GenomicsDB：C ++中的高性能数据存储，用于导入，查询和转换带有CC ++ JavaSpark绑定的变体数据。 在gatk4中使用
02-05
主开发 GenomicsDB（最初来自 ）建立在和基于图块的阵列存储系统的基础上，用于导入，查询和转换变体数据。 变异数据本质上是稀疏的（相对于整个基因组稀疏），并且使用稀疏数组数据存储非常适合存储此类数据。 ...
开箱即用版本 满分室间质评之GATK Somatic SNV+Indel+CNV+SV（20240430更新）
SliverWorkspace
10-07 445
使用校准表对bam碱基质量校准，因为这一步gatk效率感人，所以同时计算insertsize，拆分interval list（后续mutect2并行运行需要），运行cnvkit batch，运行samtools depth计算测序深度，samtools flagstat 统计mapping比例及质量。编写脚本处理cnvkit输出，计算cnv基因，exon位置，gain/lost，cn数。编写脚本处理manta的输出，获取最终sv输出结果，起始位置，基因、频率等。使用cnvkit提供工具输出分布图和热图。
gatk-sv
03-19
GATK-SV 用于Illumina短读全基因组测序（WGS）数据的结构变异发现管道。 目录 原始呼叫者和证据收集 批量质量控制 -gCNV模型创建 批量证据合并，BAF生成和深度调用者 网站集群 网站指标 过滤 跨批次站点合并 基因分型 基因型优化（可选） 跨批次集成，复杂事件解决和VCF清理 下游过滤 注释 质量控制和可视化 附加模块-Mosaic和de novo 部署和执行： 一个帐户。 支持（WDL）的流执行系统，可以： （v36或更高版本）。强烈建议使用专用服务器。 或 （请注意，此平台尚不支持预配置的GATK-SV工作流程） 推荐： 。由于许可限制，我们无法为此算法提供公共docker映像或参考面板VCF。 推荐：用于与专用Cromwell服务器进行交互的 。 推荐：用于验证WDL / json文件的 。 数据： Illumina短读全基因组CRAM或BAM，与带有hg38
探秘GATK-SV：结构变异发现的利器
最新发布
gitblog_00081的博客
06-09 422
探秘GATK-SV：结构变异发现的利器 项目地址:https://gitcode.com/broadinstitute/gatk-sv 在基因组研究领域，结构变异（SV）的准确检测是揭示遗传多样性与疾病关联的关键一环。GATK-SV 是一个专为Illumina短读测序（WGS）数据设计的结构变异发现管道，它提供了一套全面且高效的方法来处理大规模样本集中的复杂遗传信息。 项目简介 GATK-SV 使...
GATK Somatic 流程学习1
coolsen133的专栏
01-11 819
1. 主要是参考链接https://software.broadinstitute.org/gatk/documentation/article?id=9183#1 的步骤，阅读了解相应的介绍 2. https://console.cloud.google.com/storage/browser/gatk-best-practices/somatic-b37/?project=broad-dsd...
满分室间质评之GATK Somatic SNV+Indel+CNV+SV（下）性能优化
SliverWorkspace
08-04 3019
我们接上文：满分室间质评之GATK Somatic SNV+Indel+CNV+SV一文中实现了对于卫计委室间质评数据分析以及与满分结果的匹配。本文将着重解决，保证最终结果一致的情况下，如何优化分析性能（并行化），如何将分析时间从 3h 59m 53s缩短至 1h 10m 38s。 优化的方向：实际运行GATK4.X的工具如Mutect2时，发现其运行效率相当低，从CPU占用率，内存占用，硬盘I/O都占用很低，起初自己DIY时候，将要分析的bed/interval_list文件按照染色体编号拆分（不太确定
使用gatk检测WES数据中的cnv
Q.1的博客
04-26 6258
gatk的cnv流程对环境依赖较高，需要调用许多python包，推荐在dockerhub里找官方镜像，或者用conda来配置环境。 1、dockerhub 在本地的docker环境中直接拉取镜像，如果没有root权限就用conda安装。 docker pull broadinstitute/gatk:4.1.6.0 2、conda 先下载一个miniconda或者anaconda，然后下载好ga...
CallVars：CallVars是一种自动的NGS工作流程，可将配对末端的样品FastQ直接带到感兴趣的少数临床相关变体（SNPs + Indels）
02-23
该工作流基于准则，即针对单个样品的种系短变异发现（SNP + Indel），并使用基因组分析工具包（GATK）v4.1.9.0。 CallVars配置为使用“ config.yaml”文件（此存储库中随附）中列出的参数运行。您可以更改config...
atlas_gatk4:在Atlas HPC上测试gatk4模块
04-07
Atlas HPC上的GATK4模块 在Atlas HPC上测试gatk4模块 示例数据集 重用GATK nextflow管道中的数据集（ ） wget https://iastate.box.com/shared/static/wt85l6s4nw4kycm2bo0gpgjq752osatu.gz tar -xf wt85l6s4nw4...
gatk:GATK版本4及更高版本的官方代码存储库
05-11
使用Google Cloud Storage上的输入来运行GATK4 在本地运行GATK4 Spark工具 在Spark集群上运行GATK4 Spark工具 在Google Cloud Dataproc上运行GATK4 Spark工具 使用R生成图 Bash的GATK标签完成
FishingCNV:外显子组测序数据中的拷贝数变异（CNV）检测-开源
04-29
FishingCNV是由麦吉尔大学开发的软件工具，是用于全面分析高通量外显子组测序数据中稀有拷贝数变异的工具。 输入是由Genome Analysis ToolKit（GATK）生成的标准coverage文件，而输出是包含推定CNV的文件。 该程序有2个不同的发行版* GUI版本（FishingCNV_X.XXzip）*命令行版本（FishingCNV_X.X_pipeline）浏览我们的文件以查找不同的发行版。
GATK使用方法详细介绍
07-31
GATK软件使用方法详细描述，好的资源不容错过！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！
GenomeAnalysisTK（GATK3.8.0）便携版，无需编译直接使用
04-29
包含GenomeAnalysisTK.jar和README.md文件。 Linux使用 java -jar GenomeAnalysisTK.jar -h 可查看帮助文档，直接运行GATK3.8.0软件。
gatk-cwl-generator:从GATK文档生成CWL文件
04-27
gatk-cwl生成器 从生成文件 安装 一，安装模块 git clone https://github.com/wtsi-hgi/gatk-cwl-generator cd gatk-cwl-generator python setup.py install 您可能还需要安装来运行生成的CWL文件 要求 的Python 3 ...
GATK HaplotypeCaller SNP Calling 自动化流程
03-13
使用GATK进行 SNP Callling的自动化流程。基于GATK3.3以上版本的HaplotypeCaller标准流程进行，测试脚本以sra文件为最初输入，将读取文件见内的所有sra文件（视为同一个样本的数据），进行SNP Calling产生GVCF，然后...
GATK变异检测
weixin_45606903的博客
08-13 9173
变异检测 GATK 变异检测 简述：GATK(全称The Genome Analysis Toolkit) ,是基于二代测序的分析软件，主要功能是寻找变异位点。 数据类型：Illumina数据 软件版本：Gatk 4.1.8.1、fastp 0.20.0、bwa 0.7.17、samtools 1.9 测试数据： $ ref = /home/jingke.meng/database/hg19_BWA /hg19.fa $ fq1 = /home/jingke.meng/GATK_passway/Illu
GATK4 Mutect2和HaplotypeCaller的不同
tangerine1222的博客
04-11 1148
也就是说，如果病例中的一个变异点在匹配的对照组或人群资源中也是变异点，如dbSNP，即使变异的等位基因与对照组或资源不同，也会从体细胞调用中扣除。该工具不考虑与正常人小组或与匹配的正常对照的明确的变异部位共享的变异部位。根据与上述类似的逻辑，并假设常见的变异位点是双联的，任何在病例样本中出现多个变异等位基因的位点都是可疑的。其中，不在生殖系资源中的等位基因的频率（--af-of-alleles-not-in-resource）在生殖系变体先验中定义，Mutect2在计算变体是生殖系的可能性时使用。
GATK4最佳实践-体细胞突变的检测与识别
庐州月光的博客
05-28 3155
欢迎关注"生信修炼手册"！分析体细胞突变时，通常采用tumor_vs_nomal 的实验设计。在检测时，由于同时会检测出生殖细胞突变和体细胞突变，需要做的就是去除生殖细胞...
spark gatk
01-24
Spark GATK是指在Spark集群上运行GATK（Genome Analysis Toolkit）工具。GATK是一个广泛使用的工具集，用于分析基因组数据。它提供了一系列的工具，用于变异检测、变异过滤、变异注释等。 在Spark集群上运行GATK...
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值