最新版针对RNA-seq数据的GATK找变异流程
RNA-seq标准分析，我们已经讲解的太多了，表达矩阵到差异分析等下游生物学注释都没有啥新颖之处，融合基因和可变剪切算是出彩的地方，如果加上GATK找变异流程就更棒了，反正都使用了star软件进行序列比对拿到bam文件了。
如果你简单谷歌搜索关键词：gatk best practices pipeline rna-seq 会搜索到大量过期的教程：
2014年5月6日 - Best Practices workflow for RNAseq https://gist.github.com/PoisonAlien/c6c03539cf4b1ac41cf1
2016年10月27日 - https://www.biostars.org/p/219281/
2016年10月17日 - https://gatkforums.broadinstitute.org/gatk/discussion/3892/the-gatk-best-practices-for-variant-calling-on-rnaseq-in-full-detail
2017年3月17日 - 2017 Mar 17. doi: 10.12688/wellcomeopenres.10501.2
因为软件和数据库都是在持续更新，所以我们必须得紧跟潮流。
我们需要看最新的：Best Practices Workflows | Created 2018-01-09 | Last updated 2019-07-11 ，链接是 https://software.broadinstitute.org/gatk/best-practices/workflow?id=11164
回顾一下star的2-pass比对代码
我在介绍star-fusion：最好用的融合基因查找工具终于正式发表了 的时候，附带了比对代码，核心代码是：
start=$(date +%s.%N)
echo star start `date`
$bin_star --runThreadN 4 --genomeDir $star_index \
--twopassMode Basic --outReadsUnmapped None --chimSegmentMin 12 \
--alignIntronMax 100000 --chimSegmentReadGapMax parameter 3 \
--alignSJstitchMismatchNmax 5 -1 5 5 \
--readFilesCommand zcat --readFilesIn $fq1 $fq2 --outFileNamePrefix ${sample}_
mv ${sample}_Aligned.out.sam $sample.sam
$bin_samtools sort -o $sample.bam $sample.sam
$bin_samtools index $sample.bam
$bin_samtools flagstat $sample.bam > $sample.flagstat
# 这里需要判断上一个步骤是否成功，判断命令状态
touch ok.star.$sample.status
rm $sample.sam
echo star end `date`
dur=$(echo "$(date +%s.%N) - $start" | bc)
printf "Execution time for star : %.6f seconds" $dur
实际上就是一行命令在运行比对过程，但是呢，参数太多了，调起来很麻烦，通常如果不理解的话就不建议修改参数。如果你确实需要做star-fusion 有一个链接需要注意：https://github.com/STAR-Fusion/STAR-Fusion/issues/104 我们这里仅仅是需要star比对后的bam文件来走GATK找变异流程。
参考基因组都使用star-fusion的31G数据库文件里面的：
~/biosoft/starFusion/db/GRCh38_gencode_v31_CTAT_lib_Oct012019.plug-n-play/ctat_genome_lib_build_dir
值得注意的是，因为下载star-fusion的31G数据库文件解压后只有 ref_genome.fa ，并没有 ref_genome.dict，需要自己构建。
Picard CreateSequenceDictionary creates .dict file and samtools faidx creates a .fai file. Both are needed for GATK. 自行搜索如何构建哈。
下载最新版GATK
官网：https://software.broadinstitute.org/gatk/ 目前最新版是：4.1.4.0 我后面的教程都是基于此：
地址是：https://github.com/broadinstitute/gatk/releases/download/4.1.4.0/gatk-4.1.4.0.zip
mkdir -p ~/biosoft
cd ~/biosoft
mkdir gatk
wget https://github.com/broadinstitute/gatk/releases/download/4.1.4.0/gatk-4.1.4.0.zip
unzip gatk-4.1.4.0.zip
主要流程
首先去除PCR重复
conda activate rna
我们对bam文件取部分作为测试数据
samtools view -h ../star/SRR2016932.bam|head -100000 > tmp.sam
samtools sort -O bam -@ 4 tmp.sam -o tmp.bam
samtools index tmp.bam
conda install -y -c bioconda sambamba
sample=tmp
sambamba markdup --overflow-list-size 600000 --tmpdir='./' -r $sample.bam ${sample}_rmd.bam
然后SplitNCigarReads
这个时候需要参考基因组文件，以及配套的dict好fai文件。
gatk='/home/jmzeng/biosoft/gatk/gatk-4.1.4.0/gatk'
GENOME='/home/jmzeng/biosoft/starFusion/db/GRCh38_gencode_v31_CTAT_lib_Oct012019.plug-n-play/ctat_genome_lib_build_dir/ref_genome.fa'
$gatk CreateSequenceDictionary -R $GENOME # 最新版gatk，这个步骤非常快
$gatk SplitNCigarReads -R $GENOME \
-I ${sample}_rmd.bam \
-O ${sample}_rmd_split.bam
这样得到的bam文件，才跟DNA流程的类似，可以走gatk后续流程了。
接着走Base Quality Recalibration
gatk='/home/jmzeng/biosoft/gatk/gatk-4.1.4.0/gatk'
DBSNP=/home/jmzeng/biosoft/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz
kgSNP=/home/jmzeng/biosoft/GATK/resources/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
kgINDEL=/home/jmzeng/biosoft/GATK/resources/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
GENOME='/home/jmzeng/biosoft/starFusion/db/GRCh38_gencode_v31_CTAT_lib_Oct012019.plug-n-play/ctat_genome_lib_build_dir/ref_genome.fa'
# 因为我的star比对得到的bam文件里面没有 Read groups
# 参考：https://gatkforums.broadinstitute.org/gatk/discussion/6472/read-groups
$gatk AddOrReplaceReadGroups -I ${sample}_rmd_split.bam \
-O ${sample}_rmd_split_add.bam -LB $sample -PL ILLUMINA -PU $sample -SM $sample
$gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" BaseRecalibrator \
-I ${sample}_rmd_split_add.bam \
-R $GENOME \
-O ${sample}_recal.table --known-sites $kgSNP --known-sites $kgINDEL
$gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" ApplyBQSR \
-I ${sample}_rmd_split_add.bam -R $GENOME --output ${sample}_recal.bam -bqsr ${sample}_recal.table
最后Variant Calling
因为最后所有的样本需要合并比较，所以这里使用gvcf模式的Variant Calling
bed=/home/jmzeng/annotation/CCDS/human/exon_probe.GRCh38.gene.150bp.bed
# ftp://ftp.ncbi.nlm.nih.gov/snp/organisms/human_9606/VCF/All_20180418.vcf.gz
bam=${sample}_recal.bam
$gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" HaplotypeCaller \
-ERC GVCF -L $bed -R $GENOME -I $bam --dbsnp $DBSNP -O ${sample}_gatk.gvcf
一个样本的star比对后的bam文件，走RNA-seq数据的GATK找变异流程得到的全部文件如下：
如果是普通的vcf文件，是需要过滤一下的，官网是这样说的：We recommend specific hard filters, since VQSR and CNNScoreVariants require truth data for training that we don’t yet have for RNA.
因为我这里走的是gvcf模式，所以暂时不需要这个过滤步骤。
串起来成为shell脚本流程
只需要把我们的全部star的bam文件全路径保存为bam_path.txt 文件，就可以使用下面的脚本提交任务批量运行我们的gvcf流程啦，当然也可以控制并行的样本数量，具体需要理解shell脚本的语法了。
cat bam_path.txt |while read id
do
file=$(basename $id )
sample=${file%%.*}
if((i%$number1==$number2))
then
if [ ! -f ok.gvcf.$sample.status ]; then
time sambamba markdup --overflow-list-size 600000 --tmpdir='./' -r $id ${sample}_rmd.bam
time $gatk SplitNCigarReads -R $GENOME \
-I ${sample}_rmd.bam \
-O ${sample}_rmd_split.bam
time $gatk AddOrReplaceReadGroups -I ${sample}_rmd_split.bam \
-O ${sample}_rmd_split_add.bam -LB $sample -PL ILLUMINA -PU $sample -SM $sample
time $gatk BaseRecalibrator \
-I ${sample}_rmd_split_add.bam \
-R $GENOME \
-O ${sample}_recal.table --known-sites $kgSNP --known-sites $kgINDEL
time $gatk ApplyBQSR \
-I ${sample}_rmd_split_add.bam -R $GENOME --output ${sample}_recal.bam -bqsr ${sample}_recal.table
time $gatk HaplotypeCaller \
-ERC GVCF -L $bed -R $GENOME -I ${sample}_recal.bam --dbsnp $DBSNP -O ${sample}_gatk.gvcf
if [ $? -eq 0 ]; then
echo "succeed" $sample
touch ok.gvcf.$sample.status
rm ${sample}_rmd.bam ${sample}_rmd_split.bam ${sample}_rmd_split_add.bam
rm ${sample}_rmd.bam.bai ${sample}_rmd_split.bam.bai
else
echo "failed" $sample
fi # check gvcf status
fi ## check whether we need to process current sampel
fi # check the batch
i=$((i+1))
done
一个真正的样本测序数据走这个针对RNA-seq数据的GATK找变异流程会输出的文件非常多，而且耗时很夸张的！
下面的5个文件是需要删除的：
得到vcf文件，目前不显示了。
大家注意一下java问题
仔细看我分步骤调试代码的时候，有时候加上了 " --java-options -Xmx20G -Djava.io.tmpdir=./" 的参数，但是合并起来成为脚本的时候，又删除它了，具体细节看：虽然不知道为什么但是我可以解决这个bug，但它就是被解决了