GATK RNA-Seq Snps Indel 分析

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

Gatk RNA -Seq 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 -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

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

名称（点击下载）	备注
SliverWorkspace 2.0.295368	提供运行控制平台/社区版
GATK Germline SNP/Indel V1.2	分析流程文件，可以一键导入分析平台(点击查看操作) 不想复制shell的，可以使用平台一键导入流程，当然reference文件和软件还需要自己下载和安装
ucsc.hg19.gtf.tar.xz	ucsc.h19.gtf ucsc.hg19.gtf.bed 从ucsc.hg19.gtf中列数据中生成的bed文件 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" ]; thenmkdir -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`;
dorm -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`;
dovcfs="-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都非常占用时间。也难怪市面上各种加速方案了。