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

2021-07-29 11:28:02 浏览数 (1)

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

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

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

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

代码语言：javascript复制
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添加的套件进行转换。

代码语言：javascript复制
gatk SplitNCigarReads 
-R ~/reference/linux/STAR/STAR_GRCh38_genecode_v33/ref_genome.fa 
-I ./SRR11178348_dedup.bam 
-O ./SRR11178348_dedup_split.bam

3、进行碱基质量重校正

代码语言：javascript复制
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

代码语言：javascript复制
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进行过滤

代码语言：javascript复制
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文件

代码语言：javascript复制
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官方文档》

spark na

0 人点赞