专栏首页生信菜鸟团RNA-seq数据分析完全指北-10:gatk找突变

RNA-seq数据分析完全指北-10:gatk找突变

如果有读者仔细看过RNA-seq结题报告,就会发现在定量分析以外通常还会有SNP和INDEL分析。目前,对人类测序数据找突变最常用的软件是GATK,除了速度慢以外,没有其他明显缺点(可以通过部署Spark提高速度;当然,如果有钱,可以购买Sentieon,快了15-20倍)。

和WES不同,RNA-seq对于外显子区域的覆盖度极度不均一,并且由于其数据来自转录本,跨越剪切位点的对比方式也使其与常规全外显子测序大不相同。对于RNA-seq与WES进行变异检测的异同,黄树嘉老师在《RNA-Seq能替代WES完成外显子的变异检测吗?》一文中有详细解释。

本文后续内容主要来自GATK官网,但是这一部分的分析需要对WES检测变异具有一定的基础。如果跟不上下面的代码,墙裂推荐生信技能树的《肿瘤外显子专栏》,对于变异检测进行一个系统的了解之后再来follow。

1、向bam文件添加头文件并去重

gatk AddOrReplaceReadGroups \
-I ./SRR11178348Aligned.sortedByCoord.out.bam \
-O ./SRR11178348_rg_added_sorted.bam \
-ID SRR11178348 \
-LB rna \
-PL illumina \
-PU hiseq \
-SM SRR11178348 && \
gatk MarkDuplicates \
-I ./SRR11178348_rg_added_sorted.bam \
-O ./SRR11178348_dedup.bam  \
-M ./SRR11178348_dedup.metrics \
--CREATE_INDEX true \
--VALIDATION_STRINGENCY SILENT

对WES流程有了解的同学回想bwa-men流程,就会发现在bwa对比的那一步已经添加了头文件,在此我们向STAR对比的bam文件追加这一部分。

2、转换MAPQ

由于STAR的MAPQ和GATK的MAPQ并不一致,所以需要使用GATK专门为RNA-seq添加的套件进行转换。

gatk SplitNCigarReads \
-R ~/reference/linux/STAR/STAR_GRCh38_genecode_v33/ref_genome.fa \
-I ./SRR11178348_dedup.bam \
-O ./SRR11178348_dedup_split.bam

3、进行碱基质量重校正

gatk BaseRecalibrator \
-R ~/reference/linux/STAR/STAR_GRCh38_genecode_v33/ref_genome.fa \
-I ./SRR11178348_dedup_split.bam \
--known-sites ~/reference/linux/gatk/GRCh38/1000G_phase1.snps.high_confidence.hg38.vcf \
--known-sites ~/reference/linux/gatk/GRCh38/dbsnp_146.hg38.vcf \
--known-sites ~/reference/linux/gatk/GRCh38/Mills_and_1000G_gold_standard.indels.hg38.vcf \
-O ./SRR11178348_recal_data.table && \
gatk ApplyBQSR \
-R ~/reference/linux/STAR/STAR_GRCh38_genecode_v33/ref_genome.fa \
-I ./SRR11178348_dedup_split.bam \
-bqsr ./SRR11178348_recal_data.table \
-O ./SRR11178348_BQSR.bam

4、call mutations

gatk HaplotypeCaller \
-R ~/reference/linux/STAR/STAR_GRCh38_genecode_v33/ref_genome.fa \
-I ./SRR11178348_BQSR.bam \
-O ./SRR11178348.vcf \
-L ~/reference/linux/gatk/GRCh38/GRCh38.interval_list

5、对mutations进行过滤

gatk VariantRecalibrator \
-R ~/reference/linux/STAR/STAR_GRCh38_genecode_v33/ref_genome.fa \
-V ./SRR11178348.vcf \
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 ~/reference/linux/gatk/GRCh38/hapmap_3.3.hg38.vcf \
-resource:omini,known=false,training=true,truth=false,prior=12.0 ~/reference/linux/gatk/GRCh38/1000G_omni2.5.hg38.vcf \
-resource:1000G,known=false,training=true,truth=false,prior=10.0 ~/reference/linux/gatk/GRCh38/1000G_phase1.snps.high_confidence.hg38.vcf \
-resource:dbsnp,known=true,training=false,truth=false,prior=6.0 ~/reference/linux/gatk/GRCh38/dbsnp_146.hg38.vcf \
-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 \
--tranches-file SRR11178348.snps.tranches \
-O SRR11178348.snps.recal && \
gatk ApplyVQSR \
-R ~/reference/linux/STAR/STAR_GRCh38_genecode_v33/ref_genome.fa \
-V ./SRR11178348.vcf.gz \
--tranches-file SRR11178348.snps.tranches \
--recal-file SRR11178348.snps.recal \
-ts-filter-level 99.0 \
-mode SNP \
-O SRR11178348.snps.VQSR.vcf.gz

gatk VariantRecalibrator \
-R ~/reference/linux/STAR/STAR_GRCh38_genecode_v33/ref_genome.fa \
-V ./SRR11178348.snps.VQSR.vcf.gz \
-resource:mills,known=true,training=true,truth=true,prior=12.0 ~/reference/linux/gatk/GRCh38/Mills_and_1000G_gold_standard.indels.hg38.vcf \
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
-mode INDEL \
--max-gaussians 6 \
--tranches-file SRR11178348.indels.tranches \
-O SRR11178348.snps.indels.recal && \
gatk ApplyVQSR \
-R ~/reference/linux/STAR/STAR_GRCh38_genecode_v33/ref_genome.fa \
-V ./SRR11178348.snps.VQSR.vcf.gz \
--tranches-file SRR11178348.indels.tranches \
--recal-file SRR11178348.snps.indels.recal \
-ts-filter-level 99.0 \
-mode INDEL \
-O SRR11178348.VQSR.vcf.gz

6、使用VEP进行注释并转换为maf文件

vep --cache --species homo_sapiens \
--dir ~/reference/linux/vep/GRCh38 \
--clin_sig_allele 0 \
--af --af_gnomad \
-R ~/reference/linux/STAR/STAR_GRCh38_genecode_v33/ref_genome.fa \
-i SRR11178348.VQSR.vcf --vcf \
-a GRCh38 \
-o SRR11178348.anno.vcf

~/biosoft/vcf2maf/vcf2maf.pl \
--input-vcf SRR11178348.anno.vcf \
--output-maf SRR11178348.maf \
--normal-id SRR11178348 \
--tumor-id SRR11178348 \
-R ~/reference/linux/STAR/STAR_GRCh38_genecode_v33/ref_genome.fa \
--vep-path ~/miniconda3/envs/vep/share/ensembl-vep-104.3-0/ \
--vep-data ~/reference/linux/vep/GRCh38 \
--ncbi-build GRCh38

其实,在获得VCF文件之后,还可以使用数据库文件或自有的测序文件进行germline突变的去除,以进一步提高结果的可信度。

之后,就可以使用maftools对结果进行可视化以及其他分析啦~~

部分内容引用自下列文章

《RNA-seq 检测变异之 GATK 最佳实践流程》——生信技能树,王鹏 《RNA-Seq能替代WES完成外显子的变异检测吗? 》——碱基矿工,黄树嘉 《GATK官方文档》

本文分享自微信公众号 - 生信菜鸟团(bio_123456789),作者:开啤酒摊的老程

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2021-07-18

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • CNS图表复现11—RNA-seq数据可不只是表达量矩阵结果

    CNS图表复现之旅前面我们已经进行了10讲,你可以点击图表复现话题回顾。如果你感兴趣也想加入交流群,自己去:你要的rmarkdown文献图表复现全套代码来了(单...

    生信技能树jimmy
  • 用WES和RNA-Seq数据提取到的somatic SNVs不一致

    其实两者均可用于检测遗传变异,特别是在单核苷酸变异方面(SNVs)。如果大家对RNA-seq数据如何找变异位点的流程不是很清楚,可以看我们生信技能树以前的教程:

    生信技能树
  • Variant 分析阶段小结2- 变异寻找碎碎念

    写在前面:『思考问题的熊』专栏上次更新还要追溯到4月19号的 Variant 分析阶段小结1-基础碎碎念,过去接近一个月的时间里我分别经历了两次长途出差和电脑无...

    生信技能树
  • 最新版针对RNA-seq数据的GATK找变异流程

    如果你简单谷歌搜索关键词:gatk best practices pipeline rna-seq 会搜索到大量过期的教程:

    生信技能树
  • RNA-seq数据分析完全指北-03:去除奇怪的RNA

    不同于我们常见的polyA富集方法,号称全转录组测序的rRNA depletion建库对于实验的要求更高,并且在建库过程中引入的我们并不想分析的序列也更多。

    生信菜鸟团
  • RNA-seq 检测变异之 GATK 最佳实践流程

    RNA-seq 序列比对 对 RNA-seq 产出的数据进行变异检测分析,与常规重测序的主要区别就在序列比对这一步,因为 RNA-seq 的数据是来自转录本的,...

    生信技能树
  • GATK RNA-Seq Snps Indel 分析

    https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP058243&o=acc_s%3Aa

    SliverWorkspace
  • 单个肿瘤病人的外显子数据分析策略

    但是实际上肿瘤外显子队列是很烧钱的,通常来说,一个肿瘤病人需要测50X的血液加上200X的肿瘤,基本上3000块钱是跑不了的,100人的队列就是三十好几万了。而...

    生信技能树
  • 多靶点自体免疫细胞技术

    本期文章题目是 Immune recognition of somatic mutations leading to complete durable regr...

    生信技能树jimmy
  • RNA-seq数据分析完全指北-02:fastq文件质量控制

    说了这么多闲话,下面进入正题。在我们终于拿到了fastq文件之后,又该做些什么呢?

    生信菜鸟团
  • RNA-seq数据分析完全指北-05:去接头以及过滤

    一般来说,测序结果如下如所示,包括barcode和部分insert。而barcode部分在demultiplexing会被去除,剩下的就只有测到的一部分inse...

    生信菜鸟团
  • 【生信文献200篇】27 多靶点自体免疫细胞技术

    「英文标题」 Immune recognition of somatic mutations leading to complete durable regre...

    生信菜鸟团
  • RNA-seq数据分析完全指北-04:创建本地blast库分析物种组成

    啊~~~本来是半个月的专栏不知道到底过了多久才又和大家见面,其中经历不足为外人道也

    生信菜鸟团
  • 生信老司机以中心法则为主线讲解组学技术的应用和生信分析心得 - 限时免费

    海哥,中国科学院遗传与发育生物学研究所,生物信息学博士。在生信宝典出品过多部“傻瓜式”教程。

    生信宝典
  • 图形化开放式生信分析系统开发 - 3 生信分析流程的进化

    接上两篇内容,本文主要讲述工作中NGS从科研进入医学临床领域,工作中接触到生信流程,以及最终在实现的过程。

    SliverWorkspace
  • GATK4基本概念整理

    GATK 是 Genome Analysis ToolKit 的缩写,是一款从高通量测序数据中分析变异信息的软件,是目前最主流的snp calling 软件之一...

    生信修炼手册
  • 生物信息学流程框架的4个流派

    比如Nextflow、Snakemake等等,这方面的各种教程多如牛毛,我这里就不赘述了,大家根据关键词搜索即可自行学习。

    生信技能树
  • 找寻gvcf失败的原因(java啊java,配置啊配置)

    因为我只想保留recal后的bam和,call出来的gvcf文件,但是发现有些样本根本就走不通这个流程,就需要debug了。

    生信技能树
  • 转录组高级分析之融合基因

    2017年BMC文章:De novo assembly and characterization of breast cancer transcriptomes...

    生信技能树

扫码关注云+社区

领取腾讯云代金券