RNA-seq 检测变异之 GATK 最佳实践流程

RNA-seq 序列比对

对 RNA-seq 产出的数据进行变异检测分析,与常规重测序的主要区别就在序列比对这一步,因为 RNA-seq 的数据是来自转录本的,比对到参考基因组需要跨越转录剪切位点,所以 RNA-seq 进行变异检测的重点就在于跨剪切位点的精确序列比对。

文献 systematic evaluation of spliced alignment programs for RNA-seq data 中对 RNA-seq 数据常用的 11 款比对软件进行了详细测试,包括 STAR 2-pass,而 GATK 对 RNA-seq 数据变异检测的最佳实践流程中选用了 STAR 2-pass 这一方法进行比对,STAR 发表的文章至今已被引用 1900 余次,这款软件的比对速度很快,也是 ENCODE 项目的御用比对软件。

STAR 2-pass 模式需要进行两次序列比对,建立两次参考基因组索引。它的思路是第一次建参考基因组索引之后进行初步的序列比对,根据初步比对结果得到该样本所有的剪切位点信息,包括参考基因组注释 GTF 中已知的剪切位点和比对时新发现的剪切位点,然后利用第一次比对得到的剪切位点信息重新对参考基因组建立索引,然后进行第二次的序列比对,这样可以得到更精确的比对结果。

这里使用了一个测试数据演示流程

  • 第一次对参考基因组建索引:
# star 1-pass indexSTAR --runThreadN 8 --runMode genomeGenerate \        --genomeDir ./star_index/ \        --genomeFastaFiles ./genome/chrX.fa \        --sjdbGTFfile ./genes/chrX.gtf
  • 然后进行第一次序列比对:
#star 1-pass alignSTAR --runThreadN 8 --genomeDir ./star_index/ \        --readFilesIn ./samples/ERR188044_chrX_1.fastq.gz ./samples/ERR188044_chrX_2.fastq.gz \        --readFilesCommand zcat \        --outFileNamePrefix ./star_1pass/ERR188044
  • 之后根据第一次比对得到的所有剪切位点,重新对参考基因组建立索引:
# star 2-pass indexSTAR --runThreadN 8 --runMode genomeGenerate \        --genomeDir ./star_index_2pass/ \        --genomeFastaFiles ./genome/chrX.fa \        --sjdbFileChrStartEnd ./star_1pass/ERR188044SJ.out.tab
  • 再进行 STAR 二次序列比对:
# star 2-pass alignSTAR --runThreadN 8 --genomeDir ./star_index_2pass/ \        --readFilesIn ./samples/ERR188044_chrX_1.fastq.gz ./samples/ERR188044_chrX_2.fastq.gz \        --readFilesCommand zcat \        --outFileNamePrefix ./star_2pass/ERR188044

由于后面要用 GATK 进行 call 变异,还需要对比对结果 SAM 文件进行一些处理,这些都可以用 picard 来做,包括 SAM 头文件添加 @RG 标签,SAM 文件排序并转 BAM 格式,然后标记 duplicate:

# picard Add read groups, sort, mark duplicates, and create indexjava -jar picard.jar AddOrReplaceReadGroups \        I=./star_2pass/ERR188044Aligned.out.sam \        O=./star_2pass/ERR188044_rg_added_sorted.bam \        SO=coordinate \        RGID=ERR188044 \        RGLB=rna \        RGPL=illumina \        RGPU=hiseq \        RGSM=ERR188044 java -jar picard.jar MarkDuplicates \        I=./star_2pass/ERR188044_rg_added_sorted.bam \        O=./star_2pass/ERR188044_dedup.bam  \        CREATE_INDEX=true \        VALIDATION_STRINGENCY=SILENT \        M=./star_2pass/ERR188044_dedup.metrics

到此序列比对就完成了。

使用 GATK 进行变异检测

感觉 GATK 里面的工具都很慢(相对于其他的软件特别慢!),都是单线程在跑,有的虽然可以设置为多线程但是感觉没啥速度上的提升,所以要有点耐心……

由于 STAR 软件使用的 MAPQ 标准与 GATK 不同,而且比对时会有 reads 的片段落到内含子区间,需要进行一步 MAPQ 同步和 reads 剪切,使用 GATK 专为 RNA-seq 应用开发的工具 SplitNCigarReads 进行操作,它会将落在内含子区间的 reads 片段直接切除,并对 MAPQ 进行调整。

DNA 测序的重测序应用中也有序列比对软件的 MAPQ 与 GATK 无法直接对接的情况,需要进行调整。

# samtools faidx chrX.fa# samtools dict chrX.fajava -jar GenomeAnalysisTK.jar -T SplitNCigarReads \        -R ./genome/chrX.fa \        -I ./star_2pass/ERR188044_dedup.bam \        -o ./star_2pass/ERR188044_dedup_split.bam \        -rf ReassignOneMappingQuality \        -RMQF 255 \        -RMQT 60 \        -U ALLOW_N_CIGAR_READS

之后就是可选的 Indel Realignment,对已知的 indel 区域附近的 reads 重新比对,可以稍微提高 indel 检测的真阳性率,如果时间紧张也可以不做,这一步影响很轻微 :)

# 可选步骤 IndelRealignjava -jar GenomeAnalysisTK.jar -T RealignerTargetCreator \        -R ./genome/chrX.fa \        -I ./star_2pass/ERR188044_dedup_split.bam \        -o ./star_2pass/ERR188044_realign_interval.list \        -known Mills_and_1000G_gold_standard.indels.hg19.sites.vcf java -jar GenomeAnalysisTK.jar -T IndelRealigner \        -R ./genome/chrX.fa \        -I ./star_2pass/ERR188044_dedup_split.bam \        -known Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \        -o ./star_2pass/ERR188044_realign.bam \        -targetIntervals ./star_2pass/ERR188044_realign_interval.list

然后还是可选的 BQSR,这一步操作主要是针对测序质量不太好的数据,进行碱基质量再校准,如果对自己的测序数据质量足够自信可以省略,2500 之后 Hiseq 仪器的数据质量都挺不错的,可以根据 FastQC 结果来决定。这一步省了又能节省时间 :)

# 可选步骤 BQSRjava -jar GenomeAnalysisTK.jar \        -T BaseRecalibrator \        -R ./genome/chrX.fa \        -I ./star_2pass/ERR188044_realign.bam \        -knownSites 1000G_phase1.snps.high_confidence.hg19.sites.vcf \        -knownSites Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \        -o ./star_2pass/ERR188044_recal_data.tablejava -jar GenomeAnalysisTK.jar  \        -T PrintReads \        -R ./genome/chrX.fa \        -I ./star_2pass/ERR188044_realign.bam \        -BQSR ./star_2pass/ERR188044_recal_data.table \        -o ./star_2pass/ERR188044_BQSR.bam

比如下面的数据就可以放心的省略这两步了:

现在终于可以进行变异检测了,GATK 官网说 HC 表现比 UC 好,所以这里用 HC 进行变异检测:

java -jar GenomeAnalysisTK.jar -T HaplotypeCaller \        -R ./genome/chrX.fa \        -I ./star_2pass/ERR188044_dedup_split.bam \        -dontUseSoftClippedBases \        -stand_call_conf 20.0 \        -o ./star_2pass/ERR188044.vcf

call 完变异之后再进行过滤:

java -jar GenomeAnalysisTK.jar \        -T VariantFiltration \        -R ./genome/chrX.fa \        -V ./star_2pass/ERR188044.vcf \        -window 35 \        -cluster 3 \        -filterName FS -filter "FS > 30.0" \        -filterName QD -filter "QD < 2.0" \        -o ./star_2pass/ERR188044_filtered.vcf

然后就拿到变异检测结果了,可以用 ANNOVAR 或 SnpEff 或 VEP 进行注释,根据自己的需要进行筛选了。

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-06-05

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏生信宝典

分子对接简明教程 (三)

Docking非原生配体 在前面的例子中,AutoDock Vina能把配体构象调整到几乎原生的构象,验证了这一预测方法的准确度。下面,我们尝试docking另...

6339
来自专栏生信技能树

你以为的可能不是你以为的

最近生信技能树管理员小朋友XZG跟我炫耀他植物的简化基因组的gvcf模式,两百个测序数据,我一直没用过这个gvcf功能,因为的确没有需求。癌症研究,关注的主要是...

1204
来自专栏编程微刊

【前端统计图】echarts实现简单柱状图项目地址下载:

1814

基因组分析工具包:Apache Spark

自2000年人类基因组计划(Human Genome Project)产生人类基因组首份草案序列以来,测序成本从几乎每个基因组的1亿美元左右急剧下降到今天的约1...

4866
来自专栏yw的数据分析

安装使用pyclone进行克隆演化推断

pyclone介绍   可以根据多个样品突变的allele frequency 和 copy number,推断出有该突变的细胞克隆所占的比例(cellular...

5315
来自专栏生信宝典

分子对接简明教程 (一)

分子对接(Molecular Docking)理论 所谓分子对接就是两个或多个分子之间通过几何匹配和能量匹配相互识别找到最佳匹配模式的过程。分子对接对酶学研究和...

58610
来自专栏生信技能树

【直播】我的基因组58:用R包SNPRelate来对我的基因型跟hapmap计划数据比较

hapmap计划的人群分布结果和千人基因组计划的分布结果来分析是一样的!【直播】我的基因组55:简单的PCA分析千人基因组的人群分布 这两个计划里面收集的样本的...

4596
来自专栏机器之心

这是一份你们需要的Windows版深度学习软件安装指南

选自Github 机器之心编译 参与:蒋思源、刘晓坤 本文从最基本的依赖项开始,依次配置了 VS 2015、Anaconda 4.4.0、CUDA 8.0.61...

4118
来自专栏生信宝典

不是原配也可以-对接非原生配体

Docking非原生配体 在前面的例子中,AutoDock Vina能把配体构象调整到几乎原生的构象,验证了这一预测方法的准确度。下面,我们尝试docking另...

2408
来自专栏生信技能树

使用sequenza软件判定肿瘤纯度

发表该软件的文章当时使用了10 个 ovarian serous carcinomas (OVCA) 和 20 个clear-cell renal cell c...

2223

扫码关注云+社区