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 条评论
登录 后参与评论

相关文章

来自专栏杨建荣的学习笔记

海量数据迁移之分区并行抽取(r2笔记53天)

在之前的章节中分享过一些数据迁移中并行抽取的细节,比如一个表T 很大,有500G的数据,如果开启并行抽取,默认数据库中并行的最大值为64,那么生成的dump文件...

2688
来自专栏CSDN技术头条

极具参考价值的MySQL性能调优技巧

摘要:针对购物旺季网站流量会对数据库造成的压力,作者给出了MySQL性能调优的一些技巧,这些技巧极具参考价值,通过这些调优,可以有效避免因为流量过大造成服务器宕...

1926
来自专栏Kubernetes

Kubernetes 1.8抢占式调

Author: xidianwangtao@gmail.com 阅读本博文前,建议先阅读解析Kubernetes 1.8中的基于Pod优先级的抢占式调度...

39313
来自专栏java架构师

Hadoop学习19--推测式执行

  所谓推测式执行,就是计算框架判断,如果有一个task执行的过慢,则会启动备份任务,最终使用原任务+备份任务中执行较快task的结果。产生原因一般是程序bug...

2649
来自专栏沃趣科技

Oracle压缩黑科技(三):OLTP压缩

原文链接:https://www.red-gate.com/simple-talk/sql/oracle/compression-in-oracle-part-...

3127
来自专栏瓜大三哥

时序分析中的基本概念和术语

1.建立保持时间 ? 2.四种时序路径 ? 第一类时序路径:从设备A的时钟到FPGA的第一级寄存器的数据输入端口 第二类时序路径:两个同步原件之间的路径,...

2058
来自专栏SDNLAB

SDNLAB技术分享(一):ODL的Service Function Chaining入门和Demo

在网络通信过程中,包含各式各样的网络服务功能。既可以包含传统的像防火墙,NAT等功能,也有包含特定的网络应用功能(Service Function)。将特定的网...

3669
来自专栏Spark学习技巧

戳破 | hive on spark 调优点

微信交流群里有人问浪尖hive on spark如何调优,当时浪尖时间忙没时间回答,这里就给出一篇文章详细聊聊。强调一下资源设置调优,这个强经验性质的,这里给出...

703
来自专栏精讲JAVA

key / value 数据库的选型

这个项目有很多 key/value 数据(约 100 GB)需要使用,使用时基本是只读的,偶尔更新时才会批量导入,且可以忍受短暂的停机导入。我一想 TiKV 和...

803
来自专栏企鹅号快讯

如何在不会导致服务器宕机的情况下,用 PHP 读取大文件

英文:Christopher Pitt ,译文:oschina www.oschina.net/translate/performant-reading-big...

1789

扫描关注云+社区