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

生信菜鸟团

发布于 2021-07-29 11:28:02

2.8K0

发布于 2021-07-29 11:28:02

如果有读者仔细看过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官方文档》

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-07-18，如有侵权请联系 cloudcommunity@tencent.com 删除

spark

本文分享自生信菜鸟团微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

spark

登录后参与评论

0 条评论

热度