原创

GATK RNA-Seq Snps Indel 分析

这是GATK Best Practice系列学习文章中的一篇,本文尝试使用:

Gatk RNA -Seq Germline spns-indels Pipeline 来分析鼻咽癌(NPT)

分析流程如下:

GATK版本的是这样的

数据

从NCBI上下载转录组数据,访问链接为:

https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP058243&o=acc_s%3Aa

第一个样本的数据下载链接如下:

Location

Name Link

NCBI

https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-5/SRR2016932/SRR2016932.1

NCBI

https://sra-downloadb.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-6/SRR2016932/SRR2016932.1

  #下载完成后获得文件为SRR2016932.1,需要转换为fastq.gz格式,这里用到NCBI sratookit工具,请自行安装
  fastq-dump -gzip --split-3 SRR2016932.1
  #得到两个fastq.gz文件
  SRR2016932_1.fastq.gz
  SRR2016932_2.fastq.gz
  #我们将文件重命名为以符合习惯
  mv SRR2016932_1.fastq.gz SRR2016932_R1.fastq.gz
  mv SRR2016932_2.fastq.gz SRR2016932_R2.fastq.gz

用到的分析系统及分析流程文件

名称 (点击下载)

备注

提供运行控制平台/社区版

分析流程文件,可以一键导入分析平台(点击查看操作) 不想复制shell的,可以使用平台一键导入流程,当然reference文件和软件还需要自己下载和安装

ucsc.h19.gtf<br />ucsc.hg19.gtf.bed 从ucsc.hg19.gtf中列数据中生成的bed文件<br />ucsc.hg19.gtf.interval_list 使用gatk IntervalToBed工具从ucsc.hg19.gtf.bed转换来的interval文件

分析结果

流程用到的变量(程序、reference文件和数据库、数值)

变量名

变量值

类型

sn

SRR2016932

字符

tools.fastqc

/opt/FastQC/fastqc

程序

tools.STAR

/opt/ref/STAR-2.7.3a/bin/Linux_x86_64_static/STAR

程序

tools.sambamba

/opt/ref/sambamba-0.7.0-linux-static

程序

tools.gatk

/opt/ref/gatk-4.1.4.1/gatk

程序

envis.read_length

100 (测序读长)

数值

envis.threads

32 (并发线程数)

数值

envis.scatter

10 (HaplotypeCaller并行数)

数值

refs.hum

/opt/ref/hg19/ucsc.hg19.fa

文件

refs.gtf

/opt/ref/RNA/ucsc.hg19.gtf

文件

refs.bed

/opt/ref/RNA/ucsc.hg19.gtf.bed

文件

refs.interval

/opt/ref/RNA/ucsc.hg19.gtf.interval_list

文件

refs.1000G

/opt/ref/hg19/1000G_phase1.indels.hg19.vcf

文件

refs.mills

/opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.vcf

文件

refs.dbsnp

/opt/ref/hg19/dbsnp_138.hg19.vcf

文件

注意:refs文件中的基因组参考序列和gtf文件以及几个vcf文件必须为同一版本,参考序列和相应的GTF文件必须为同一个网站的同一个版本,否则分析过程中会出现各种错误。很多文章推荐使用ensembl的版本,本文使用的是ucsc.hg19版本,因为之前ref文件和参考序列已经有了,只是增加了一个GTF文件,是从ucsc网站生成下载的,链接为:http://genome.ucsc.edu/cgi-bin/hgTables

分析流程:

  • 01 INPUT 输入文件
  • 02 之前用fastqc软件看过,adapter已经去掉了,这里直接开始align,Star Align
#创建一个变量star_length,值为read_length的值-1.
let star_length=${envis.read_length}-1
#如果该star_length前缀的目录不存在,说明没有创建过索引,首先创建索引
if [ ! -d "/opt/ref/RNA/$star_length-1-pass-index" ]; then
 	mkdir -p /opt/ref/RNA/$star_length-1-pass-index
    ${tools.STAR} --runMode genomeGenerate \
    --runThreadN ${envis.threads} \
    --genomeDir  /opt/ref/RNA/$star_length-1-pass-index \
    --genomeFastaFiles ${refs.hum} \
    --sjdbGTFfile      ${refs.gtf} \
    --sjdbOverhang     $star_length
fi
#使用创建的索引,执行align,为了节省空间设置了参数 --outSAMtype BAM Unsorted,也测试过sort类型,不能直接markdup
${tools.STAR} \
    --genomeDir   /opt/ref/RNA/$star_length-1-pass-index \
	--runThreadN  ${envis.threads}  \
    --readFilesIn ${data}/RNA/${sn}_R1.fastq.gz ${data}/RNA/${sn}_R2.fastq.gz \
    --readFilesCommand zcat \
    --sjdbOverhang $star_length \
    --twopassMode Basic \
    --outSAMtype BAM  Unsorted \
    --outFileNamePrefix  ${result}/${sn}.
  • 03 Sort Bam:使用了sambamba替换了samtools
  • 04-Mark duplicatate:使用了sambamba替换了gatk picard,重命名创建的索引与gatk命名一致。
  • 05 SplitNCigarReads:将落在外显子上的reads分离出来,取出N错误碱基,去除内含子区域的reads。这一步太慢了,占用整个流程一半以上运行时间,不知道有没有办法提高速度。
  • 06- 这一步添加sn样本编号等信息,前面sort如果使用samtools因为没有sn信息会报错。
  • 07- BaseRecalibrator 碱基质量校正第一步
  • 08- ApplyBQSR 应用碱基校正
  • 09- 并行运行HaplotypeCaller
#清理拆分生成的interval文件夹:类似temp_0001_of_10,防止重复运行被上次的结果干扰
rm    -rf ${result}/${sn}/temp_*
#创建${sn}目录。
mkdir -p  ${result}/${sn}
#使用Gatk IntervalListTools拆分interval,拆分数量为${envis.scatter}
${tools.gatk} IntervalListTools \
    --SCATTER_COUNT=${envis.scatter} \
    --SUBDIVISION_MODE=BALANCING_WITHOUT_INTERVAL_SUBDIVISION_WITH_OVERFLOW \
    --UNIQUE=true \
    --SORT=true \
    --INPUT=${refs.interval} \
    --OUTPUT=${result}/${sn}
#循环遍历生成的interval list文件,运行HaplotypeCaller
index=1
for i in `ls ${result}/${sn}/*/scattered.interval_list`;
do
	rm -f ${result}/${sn}_$index.vcf
	${tools.gatk}  HaplotypeCaller \
        -L $i \
        -R ${refs.hum} \
        -I ${result}/${sn}_bqsr.bam \
        -O ${result}/${sn}/${sn}_$index.vcf \
        --dbsnp ${refs.dbsnp} \
        -dont-use-soft-clipped-bases \
        --standard-min-confidence-threshold-for-calling 20 &
    let index+=1
done
#等待所有HaplotypeCaller运行结束
wait
#用生成的vcf文件名拼接输入字段
vcfs=
for z in `ls ${result}/${sn}/${sn}*.vcf`;
do
    vcfs="-I $z $vcfs"
done
echo $vcfs
#合并生成的vcf文件。
${tools.gatk} MergeVcfs \
	$vcfs \
    -O ${result}/${sn}.vcf

  • 10-VariantFiltration:对于RNA-Seq,推荐使用硬过滤,不支持VQSR。
  • 11-顺便使用fastqc看下QC结果:
  • 12-最终分析结果:

运行结果

可以看到,GATK的工具一如既往的慢,HaplotypeCaller这一步通过拆分interval并行分析,最后合并vcf,速度从1个小时以上降到了9分钟。剩下的几步,SplitNCigarReads,BaseRecalibrator,ApplyBQSR都非常占用时间。也难怪市面上各种加速方案了。

原创声明,本文系作者授权云+社区发表,未经许可,不得转载。

如有侵权,请联系 yunjia_community@tencent.com 删除。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 使用Gatk Germline spns-indels Pipeline分析遗传病(耳聋)

    这次没有拿到遗传病的室间质评的数据,直接从NCBI上找一些数据来分析。NCBI上搜索deaf,点击第一条搜索结果,最后几经跳转找到数据下载页面:https://...

    SliverWorkspace
  • 满分室间质评之GATK Somatic SNV+Indel+CNV+SV(下)性能优化

    #此处是原先Manta分析SV的步骤一,生成runWorkflow.py,因为这一不步速度很快,所以串行执行 rm -f ${result}/${sn}/r...

    SliverWorkspace
  • 转录组RNA-Seq使用docker+bioconda搭建分析环境

    近期学习转录组分析,从ncbi下载数据,转成fastq,STAR/hisat2 map到基因组上,使用featureCount拿到表达矩阵文件挺顺利的,就是到了...

    SliverWorkspace
  • 机器学习(二)--------单变量线性回归(Linear Regression with One Variable)

    这就是代价函数了 在两个参数下 真实值与求出的值的差的平方和 除以2m 其实就是求误差的平均数

    实时计算
  • iOS底层原理总结 - 探寻KVO本质

    xx_Cc
  • 马斯克挖地道公司获1.2亿美元融资,估值逼近独角兽

    据彭博社报道,挖地道为主业的“无聊公司”(Boring),获得1.2亿美元资金新融资。

    量子位
  • 上传高通量测序原始文件

    在我们发表高通量测序文章之前通常要上传测序数据到GEO数据库,现总结流程如下。 注册账户、填写MetaSheet 在NCBI GEO官网注册一个账号,然后登陆。...

    生信宝典
  • java——内存中的数组

    2、arr=new int[5],这段代码执行完成后,系统会在堆内存中为该数组分配一块内存空间,并初始化值0,如下图所示:

    说故事的五公子
  • 测序文章数据上传找哪里

    在我们发表高通量测序文章之前通常要上传测序数据到GEO数据库,现总结流程如下。 注册账户、填写MetaSheet 在NCBI GEO官网注册一个账号,然后登陆。...

    生信宝典
  • github极简指南

    入生信的坑已经3年多了,但开始github的旅程才一年多,起初主要是为了建立bioconductor中文社区而学习的,现在也在自己的github上面分享了不少代...

    生信技能树

扫码关注云+社区

领取腾讯云代金券