GATK RNA-Seq Snps Indel 分析

  #下载完成后获得文件为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

#创建一个变量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}.

#清理拆分生成的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