搜索
查看: 423|回复: 2

GATK RNA-Seq Snps Indel 分析

[复制链接]

15

主题

20

帖子

357

积分

中级会员

Rank: 3Rank: 3

积分
357
发表于 2020-9-21 17:33:06 | 显示全部楼层 |阅读模式
本帖最后由 豆浆包子 于 2020-10-7 16:11 编辑

这是GATK Best Practice系列学习文章中的一篇,本文尝试使用:Gatk RNA -Seq spns-indels Pipeline 来分析鼻咽癌(NPT)分析流程如下:
GATK版本的是这样的
数据
从NCBI上下载转录组数据,访问链接为:
第一个样本的数据下载链接如下:
[td]
Location
Name Link
NCBIhttps://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-5/SRR2016932/SRR2016932.1
NCBIhttps://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 -O ./
  #得到两个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用到的分析系统及分析流程文件[td]
名称 (点击下载)
备注
SliverWorkspace 2.0.295368提供运行控制平台/社区版
GATK Germline SNP/Indel V1.0分析流程文件,可以一键导入分析平台(点击查看操作) 不想复制shell的,可以使用平台一键导入流程,当然reference文件和软件还需要自己下载和安装
ucsc.hg19.gtf.tar.xzucsc.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文件和数据库、数值)[td]
变量名
变量值
类型
snSRR2016932字符
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_length100     (测序读长)数值
envis.threads32          (并发线程数)数值
envis.scatter10          (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都非常占用时间。也难怪市面上各种加速方案了。





回复

使用道具 举报

6

主题

72

帖子

1319

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1319
发表于 2020-10-6 12:21:21 | 显示全部楼层
我勒个去,这么长, 你居然都排版起来了,赞!!!
还是发邮件给我吧,我给你检查看看,然后发在生信技能树公众号哈!
我的邮箱是:jmzeng1314@163.com
回复 支持 反对

使用道具 举报

15

主题

20

帖子

357

积分

中级会员

Rank: 3Rank: 3

积分
357
 楼主| 发表于 2020-10-6 19:51:33 | 显示全部楼层
bioinfotrainee 发表于 2020-10-6 12:21
我勒个去,这么长, 你居然都排版起来了,赞!!!
还是发邮件给我吧,我给你检查看看,然后发在生信技能树 ...

我今天刚看到。
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2020-10-26 21:29 , Processed in 0.028115 second(s), 29 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.