搜索
查看: 4594|回复: 7

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

[复制链接]

5

主题

37

帖子

485

积分

中级会员

Rank: 3Rank: 3

积分
485
发表于 2017-6-2 17:03:03 | 显示全部楼层 |阅读模式
本帖最后由 W.Peng 于 2017-6-2 17:56 编辑

原文链接:http://www.jianshu.com/p/b400dc7c5eea
原文用 markdown 排版的,格式清晰一些。论坛里用 markdown here 插件写帖子,但是代码部分格式还是没搞定

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 中已知的剪切位点和比对时新发现的剪切位点,然后利用第一次比对得到的剪切位点信息重新对参考基因组建立索引,然后进行第二次的序列比对,这样可以得到更精确的比对结果。

这里使用了一个测试数据演示流程,第一次对参考基因组建索引:
[Bash shell] 纯文本查看 复制代码
STAR —runThreadN 8 —runMode genomeGenerate \
        —genomeDir ./star_index/ \
        —genomeFastaFiles ./genome/chrX.fa \
        —sjdbGTFfile ./genes/chrX.gtf
然后进行第一次序列比对:
[Bash shell] 纯文本查看 复制代码
STAR —runThreadN 8 —genomeDir ./star_index/ \
        —readFilesIn ./samples/ERR188044_chrX_1.fastq.gz ./samples/ERR188044_chrX_2.fastq.gz \
        —readFilesCommand zcat \
        —outFileNamePrefix ./star_1pass/ERR188044
之后根据第一次比对得到的所有剪切位点,重新对参考基因组建立索引:
[Bash shell] 纯文本查看 复制代码
STAR —runThreadN 8 —runMode genomeGenerate \
        —genomeDir ./star_index_2pass/ \
        —genomeFastaFiles ./genome/chrX.fa \
        —sjdbFileChrStartEnd ./star_1pass/ERR188044SJ.out.tab
再进行 STAR 二次序列比对:
[Bash shell] 纯文本查看 复制代码
STAR —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:
[Bash shell] 纯文本查看 复制代码
java -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

[Bash shell] 纯文本查看 复制代码
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 无法直接对接的情况,需要进行调整。
[Bash shell] 纯文本查看 复制代码
java -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 检测的真阳性率,如果时间紧张也可以不做,这一步影响很轻微
[Bash shell] 纯文本查看 复制代码
java -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

[Bash shell] 纯文本查看 复制代码
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 结果来决定。这一步省了又能节省时间
[Bash shell] 纯文本查看 复制代码
java -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.table

[Bash shell] 纯文本查看 复制代码
java -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 进行变异检测:
[Bash shell] 纯文本查看 复制代码
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 完变异之后再进行过滤:
[Bash shell] 纯文本查看 复制代码
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 进行注释,根据自己的需要进行筛选了。





上一篇:VCF文件里面比较重要的头 tag-欢迎补充
下一篇:R黑魔法合集
回复

使用道具 举报

5

主题

37

帖子

485

积分

中级会员

Rank: 3Rank: 3

积分
485
 楼主| 发表于 2017-6-2 17:27:49 | 显示全部楼层
代码部分是我手动改的格式
回复 支持 反对

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-6-2 18:52:10 | 显示全部楼层
W.Peng 发表于 2017-6-2 17:27
代码部分是我手动改的格式

辛苦了,这个代码的问题我也还在寻找解决方案,你这个帖子非常优秀
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

1

主题

43

帖子

461

积分

中级会员

Rank: 3Rank: 3

积分
461
发表于 2017-6-3 12:04:18 | 显示全部楼层
不错的帖子 值得学习!
人生若只如初见!
回复 支持 反对

使用道具 举报

1

主题

11

帖子

132

积分

注册会员

Rank: 2

积分
132
发表于 2017-8-7 21:06:01 | 显示全部楼层
有个问题想请问下,我们做RNAseq一般都会测重复,那在这种情况下我们的重复应该怎么处理呢?,是不是DNA重测中不同的lane的处理一样?
回复 支持 反对

使用道具 举报

5

主题

37

帖子

485

积分

中级会员

Rank: 3Rank: 3

积分
485
 楼主| 发表于 2017-8-9 21:31:57 | 显示全部楼层
qiling5733 发表于 2017-8-7 21:06
有个问题想请问下,我们做RNAseq一般都会测重复,那在这种情况下我们的重复应该怎么处理呢?,是不是DNA重 ...

我觉得变异检测通常都是针对单个样本的吧,多个同类样本也需要每个单独进行变异检测,再做后续分析。
回复 支持 反对

使用道具 举报

1

主题

11

帖子

132

积分

注册会员

Rank: 2

积分
132
发表于 2017-8-14 09:48:31 | 显示全部楼层
W.Peng 发表于 2017-8-9 21:31
我觉得变异检测通常都是针对单个样本的吧,多个同类样本也需要每个单独进行变异检测,再做后续分析。 ...

多谢解答
回复 支持 反对

使用道具 举报

2

主题

41

帖子

385

积分

中级会员

Rank: 3Rank: 3

积分
385
发表于 2017-8-28 21:29:07 | 显示全部楼层
精彩的帖子,谢谢分享~
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-7-18 09:32 , Processed in 0.040351 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.