GATK Germline Best Practice学习

2024-03-28 21:18
文章标签 学习 practice best gatk germline

本文主要是介绍GATK Germline Best Practice学习，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

数据是sporadic的慢病case-control的组合。想用GATK germline best practice的方法进行突变的分析。这里主要参考GATK Germline best practice的教程。1 这里用的是GATK3.7的版本，目前已经出到GATK3.8。最近4.0也发布了。
部分步骤后续补完。。。

Map to Reference

bwa mem -t 8 -M -R '@RG\tID:${name}\tLB:${name}\tPL:ILLUMINA\tPM:X10\tSM:${name}' ${INDEX} ${RAW_DATA}/${name}_1.fastq ${RAW_DATA}/${name}_2.fastq > ${WORKING_DIR}/2018rerun/processed_bam/${name}.sam
$java -Xmx20g -jar $PICARD SortSam SORT_ORDER=coordinate INPUT=${WORKING_DIR}/2018rerun/processed_bam/${name}.sam OUTPUT=${WORKING_DIR}/2018rerun/processed_bam/${name}.bam

Basic Statistics

samtools flagstat ${WORKING_DIR}/2018rerun/processed_bam/${name}.bam > ${WORKING_DIR}/2018rerun/processed_bam/${name}.flagstat &
samtools stats ${WORKING_DIR}/2018rerun/processed_bam/${name}.bam > ${WORKING_DIR}/2018rerun/processed_bam/${name}.stats &

Mark Duplicates

$java -Xmx20g -jar $PICARD MarkDuplicates INPUT=${WORKING_DIR}/2018rerun/processed_bam/${name}.bam OUTPUT=${WORKING_DIR}/2018rerun/processed_bam/${name}_marked.bam METRICS_FILE=${WORKING_DIR}/2018rerun/processed_bam/${name}.metrics
samtools index ${WORKING_DIR}/2018rerun/processed_bam/${name}_marked.bam

Base Recalibration

$java -Xmx10g -jar $gatk_jar -T BaseRecalibrator -R $INDEX -I ${WORKING_DIR}/2018rerun/processed_bam/${name}_marked.bam -knownSites $DBSNP -knownSites $mills_and_1000G_gold_standard_indel -knownSites $G1000_phase1_indel -o ${WORKING_DIR}/2018rerun/processed_bam/${name}.before.table
$java -Xmx10g -jar $gatk_jar -T BaseRecalibrator -R $INDEX -I ${WORKING_DIR}/2018rerun/processed_bam/${name}_marked.bam -knownSites $DBSNP -knownSites $mills_and_1000G_gold_standard_indel -knownSites $G1000_phase1_indel -BQSR ${WORKING_DIR}/2018rerun/processed_bam/${name}.before.table -o ${WORKING_DIR}/2018rerun/processed_bam/${name}.after.table
$java -Xmx10g -jar $gatk_jar -T AnalyzeCovariates -R $INDEX -before ${WORKING_DIR}/2018rerun/processed_bam/${name}.before.table -after ${WORKING_DIR}/2018rerun/processed_bam/${name}.after.table -plots ${WORKING_DIR}/2018rerun/processed_bam/${name}.recalibration_plots.pdf
$java -Xmx10g -jar $gatk_jar -T PrintReads -R $INDEX -I ${WORKING_DIR}/2018rerun/processed_bam/${name}_marked.bam -BQSR ${WORKING_DIR}/2018rerun/processed_bam/${name}.after.table -o ${WORKING_DIR}/2018rerun/processed_bam/${name}.recalibrated.bam

Haplotypecaller

java -Xmx10g -jar GenomeAnalysisTK.jar\-T HaplotypeCaller \
 -R GATK_boundle/hg19/genome.fa \
  -I bam/${i}.recalibrated.bam \
  --emitRefConfidence GVCF \
  -o GVCF/${i}.raw_variants.g.vcf

Joint-Call Cohort

工具：GenotypeGVCFs
在这步，会mege所有的单样本GVCF，通过GenotypeGVCFs实现。这会生成一组joint-called SNP和indel，准备后续的筛选。这种队列的分析使得突变检测更为灵敏，即使是那些比较难检测的位点。同时也会提供基因型的矩阵，包含在所有样本里所有感兴趣的位点。
这一步跑的比较快，也可以重跑如果有新的样本进来，因此可以解决N+1问题。

/work/share/software/jre1.8.0_111/bin/java -Xmx30g -jar GenomeAnalysisTK.jar -T GenotypeGVCFs -R GATK_boundle/hg19/genome.fa --variant input.list -o GenotypeGVCF.vcf

分别提取SNV和Indel.

工具：SelectVariants

/work/share/software/jre1.8.0_111/bin/java -Xmx30g -jar GenomeAnalysisTK.jar -T SelectVariants -R GATK_boundle/hg19/genome.fa --variant GenotypeGVCF.vcf -o GenotypeGVCF.snp.vcf --selectTypeToInclude SNP
/work/share/software/jre1.8.0_111/bin/java -Xmx30g -jar GenomeAnalysisTK.jar -T SelectVariants -R GATK_boundle/hg19/genome.fa --variant GenotypeGVCF.vcf -o GenotypeGVCF.indel.vcf --selectTypeToInclude INDEL

Filter Variants by Variant (Quality Score) Recalibration 2

工具：VariantRecalibrator 3, ApplyRecalibration 4
这一步可以尽量不失去真的variants，同时也意味着要对raw callset进行筛选。Variant quality score recalibration(VQSR)主要用了机器学习方法，根据突变的annotation profile去找真的突变，并赋一个VQSLOD值给每个突变（比QUAL准）。方法共有2步，第一步基于training variants建立一个模型，然后基于模型给每个variant call一个well-calibrated probability。第二步基于这个variant quality score去筛raw call set，得到a subset of calls with desired level of quality, fine-tuned to balance specificity and sensitivity.
Resource Datasets Used in two steps of VQSR:
基于human genomes，依照GATK Best Pracitce。人类的training, truth和known resource datasets都在GATK resource bundle可以下到。对于非人类，需要自己去生成包括至少truth和training resource datasets（方法在2）。
SNP的resources:
-true sites traning resource: HapMap
-True sites training resource: Omni
-Non-true sites training resource: 1000G
-Known sites resource, not used in training: dbSNP
Indel的resources:
-True sites training resource: Mills
-Known sites resource, not used in training: dbSNP

SNP
step1 VariantRecalibrator
java -Xmx30g -jar GenomeAnalysisTK.jar \
-T VariantRecalibrator \
-R GATK_boundle/hg19/genome.fa \
-input ../GenotypeGVCF.snp.vcf \
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 ~/tools/GATK/GATK_bundle/hapmap_3.3_hg19_pop_stratified_af.vcf \
-resource:omni,known=false,training=true,truth=false,prior=12.0 ~/tools/GATK/GATK_bundle/1000G_omni2.5.hg19.sites.vcf \
-resource:1000G,known=false,training=true,truth=false,prior=10.0 ~/tools/GATK/GATK_bundle/1000G_phase1.snps.high_confidence.hg19.sites.vcf \
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 ~/tools/GATK/GATK_bundle/dbsnp_138.hg19.vcf \
-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an InbreedingCoeff \
-mode SNP \
-recalFile output.recal \
-tranchesFile output.tranches \
-rscriptFile output.plots.R
step2 ApplyRecalibration
java -Xmx30g -jar GenomeAnalysisTK.jar \
-T ApplyRecalibration \
-R GATK_boundle/hg19/genome.fa \
-input ../GenotypeGVCF.snp.vcf \
-tranchesFile output.tranches \
-recalFile output.recal \
-o snp.filtered.vcf \
--ts_filter_level 99.5 \
-mode SNP
Indel
step1 VariantRecalibrator
java -Xmx30g -jar GenomeAnalysisTK.jar \
-T VariantRecalibrator \
-R GATK_boundle/hg19/genome.fa \
-input ../GenotypeGVCF.indel.vcf \
-resource:mills,known=false,training=true,truth=true,prior=12.0 ~/tools/GATK/GATK_bundle/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 ~/tools/GATK/GATK_bundle/dbsnp_138.hg19.vcf \
-an QD -an DP -an FS -an SOR -an ReadPosRankSum -an MQRankSum -an InbreedingCoeff \
-mode INDEL \
-recalFile output.recal \
-tranchesFile output.tranches \
-rscriptFile output.plots.R
step2 ApplyRecalibration
java -Xmx30g -jar GenomeAnalysisTK.jar \
-T ApplyRecalibration \
-R GATK_boundle/hg19/genome.fa \
-input ../GenotypeGVCF.indel.vcf \
-tranchesFile output.tranches \
-recalFile output.recal \
-o indel.filtered.vcf \
--ts_filter_level 99.0 \
-mode INDEL

过滤recalibrated VCF

5, 6

这篇关于GATK Germline Best Practice学习的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


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

相关文章

HarmonyOS学习(七)——UI（五）常用布局总结

HarmonyOS学习(七)——UI（五）常用布局总结

自适应布局 1.1、线性布局（LinearLayout） 通过线性容器Row和Column实现线性布局。Column容器内的子组件按照垂直方向排列，Row组件中的子组件按照水平方向排列。 属性说明space通过space参数设置主轴上子组件的间距，达到各子组件在排列上的等间距效果alignItems设置子组件在交叉轴上的对齐方式，且在各类尺寸屏幕上表现一致，其中交叉轴为垂直时，取值为Vert
阅读更多...
Ilya-AI分享的他在OpenAI学习到的15个提示工程技巧

Ilya-AI分享的他在OpenAI学习到的15个提示工程技巧

Ilya（不是本人，claude AI）在社交媒体上分享了他在OpenAI学习到的15个Prompt撰写技巧。 以下是详细的内容： 提示精确化:在编写提示时,力求表达清晰准确。清楚地阐述任务需求和概念定义至关重要。例:不用"分析文本",而用"判断这段话的情感倾向:积极、消极还是中性"。 快速迭代:善于快速连续调整提示。熟练的提示工程师能够灵活地进行多轮优化。例:从"总结文章"到"用
阅读更多...
【前端学习】AntV G6-08 深入图形与图形分组、自定义节点、节点动画（下）

【前端学习】AntV G6-08 深入图形与图形分组、自定义节点、节点动画（下）

【课程链接】 AntV G6：深入图形与图形分组、自定义节点、节点动画（下）_哔哩哔哩_bilibili 本章十吾老师讲解了一个复杂的自定义节点中，应该怎样去计算和绘制图形，如何给一个图形制作不间断的动画，以及在鼠标事件之后产生动画。（有点难，需要好好理解） <!DOCTYPE html><html><head><meta charset="UTF-8"><title>06
阅读更多...
学习hash总结

学习hash总结

2014/1/29/   最近刚开始学hash，名字很陌生，但是hash的思想却很熟悉，以前早就做过此类的题，但是不知道这就是hash思想而已，说白了hash就是一个映射，往往灵活利用数组的下标来实现算法，hash的作用：1、判重；2、统计次数；
阅读更多...
零基础学习Redis(10) -- zset类型命令使用

零基础学习Redis(10) -- zset类型命令使用

zset是有序集合，内部除了存储元素外，还会存储一个score，存储在zset中的元素会按照score的大小升序排列，不同元素的score可以重复，score相同的元素会按照元素的字典序排列。 1. zset常用命令 1.1 zadd  zadd key [NX | XX] [GT | LT]   [CH] [INCR] score member [score member ...]
阅读更多...
【机器学习】高斯过程的基本概念和应用领域以及在python中的实例

【机器学习】高斯过程的基本概念和应用领域以及在python中的实例

引言 高斯过程（Gaussian Process，简称GP）是一种概率模型，用于描述一组随机变量的联合概率分布，其中任何一个有限维度的子集都具有高斯分布 文章目录 引言一、高斯过程1.1 基本定义1.1.1 随机过程1.1.2 高斯分布 1.2 高斯过程的特性1.2.1 联合高斯性1.2.2 均值函数1.2.3 协方差函数（或核函数） 1.3 核函数1.4 高斯过程回归（Gauss
阅读更多...
【学习笔记】 陈强-机器学习-Python-Ch15 人工神经网络（1）sklearn

【学习笔记】 陈强-机器学习-Python-Ch15 人工神经网络（1）sklearn

系列文章目录 监督学习：参数方法 【学习笔记】 陈强-机器学习-Python-Ch4 线性回归 【学习笔记】 陈强-机器学习-Python-Ch5 逻辑回归 【课后题练习】 陈强-机器学习-Python-Ch5 逻辑回归（SAheart.csv） 【学习笔记】 陈强-机器学习-Python-Ch6 多项逻辑回归 【学习笔记 及 课后题练习】 陈强-机器学习-Python-Ch7 判别分析 【学
阅读更多...
系统架构师考试学习笔记第三篇——架构设计高级知识（20）通信系统架构设计理论与实践

系统架构师考试学习笔记第三篇——架构设计高级知识（20）通信系统架构设计理论与实践

本章知识考点：         第20课时主要学习通信系统架构设计的理论和工作中的实践。根据新版考试大纲,本课时知识点会涉及案例分析题(25分),而在历年考试中,案例题对该部分内容的考查并不多,虽在综合知识选择题目中经常考查,但分值也不高。本课时内容侧重于对知识点的记忆和理解,按照以往的出题规律,通信系统架构设计基础知识点多来源于教材内的基础网络设备、网络架构和教材外最新时事热点技术。本课时知识
阅读更多...
线性代数|机器学习-P36在图中找聚类

线性代数|机器学习-P36在图中找聚类

文章目录 1. 常见图结构2. 谱聚类 感觉后面几节课的内容跨越太大，需要补充太多的知识点，教授讲得内容跨越较大，一般一节课的内容是书本上的一章节内容，所以看视频比较吃力，需要先预习课本内容后才能够很好的理解教授讲解的知识点。 1. 常见图结构 假设我们有如下图结构： Adjacency Matrix:行和列表示的是节点的位置，A[i,j]表示的第 i 个节点和第 j 个
阅读更多...
Node.js学习记录（二）

Node.js学习记录（二）

目录 一、express 1、初识express 2、安装express 3、创建并启动web服务器 4、监听 GET&POST 请求、响应内容给客户端 5、获取URL中携带的查询参数 6、获取URL中动态参数 7、静态资源托管 二、工具nodemon 三、express路由 1、express中路由 2、路由的匹配 3、路由模块化 4、路由模块添加前缀 四、中间件
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2