本章介绍如何通过RNA-seq找到可能的DNA上的Single Nucleotide Variance (SNV)和Insert/Deletion (INDEL)，我们将本示例使用STAR将RNA测序数据比对到参考基因组后，使用GATK(4.0以上版本)进行SNV/INDEL 的检测，最后使用ANNOVAR对这些SNV进行注释。
1) Software introduction
1a) STAR
目前，可以用于将RNA测序数据（reads）比对到参考基因组的软件有：Bowtie、TopHat、HISAT、STAR等。
STAR在运行时候占用机器的内存较大，一般可达到20~30G，因此需要根据可用内存控制同时运行的STAR任务数量。
参考文献： Alexander Dobin , et al. STAR: ultrafast universal RNA-seq aligner Bioinformatics . 2012. 29(1): 15-21.
1b) GATK
GATK是Broad Institute开发的一款用于检测变异（SNV/INDEL）的软件，拥有较高的引用率（已有上万次引用）。
参考文献：Aaron McKenna , et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Research . 2010. 20: 1297-1303.
2) Running steps
Copy docker load -i ~/Downloads/bioinfo_snv.tar.gz
docker run -dt --name=snv -v ~/Downloads/data:/data gangxu/snv:2.0
docker exec -it snv bash 2a) Alignment
GATK要求输入的SAM/BAM文件中有Read Group信息，因此我们需要在--outSAMattrRGline中填写相应的Read Group信息，其中ID为必填项。Read Group格式详见SAM格式
Copy echo alignment start ` date `
mkdir -p /home/test/output/1.STAR_alignment/SRR5714908
/home/test/STAR-2.7.1a/bin/Linux_x86_64_static/STAR \
--genomeDir /home/test/Homo_sapiens_GRCh38_ch1_STAR_Index \
--runThreadN 1 --readFilesCommand "gunzip -c" --readFilesIn /home/test/chr1.fq.gz --outSAMtype BAM SortedByCoordinate --outSAMattrRGline ID:1 LB:exoRBase PL:ILLUMINA PU:unit1 SM:SRR5714908 --outFileNamePrefix /home/test/output/1.STAR_alignment/SRR5714908.
echo alignment end ` date `
#echo alignment start `date`
#source /BioII/lulab_b/containers/singularity/wrappers/bashrc
#STAR \
#--genomeDir /BioII/lulab_b/chenyinghui/project/Docker/SNP/reference/Homo_sapiens_GRCh38_ch1_STAR_Index \
#--runThreadN 1 \
#--readFilesIn /BioII/lulab_b/chenyinghui/project/Docker/SNP/chr1.fq \
#--readFilesCommand "gunzip -c" \
#--outSAMtype BAM SortedByCoordinate \
#--outSAMattrRGline ID:1 LB:exoRBase PL:ILLUMINA PU:unit1 SM:SRR5714908 \
#--outFileNamePrefix /BioII/lulab_b/chenyinghui/project/Docker/SNP/output/1.STAR_alignment/SRR5714908.
#echo alignment end `date` 2b) MarkDuplicates
在建库的PCR过程中会形成一些重复的DNA片段序列，这些重复序列被称为PCR duplicates。另外，在测序仪进行光学测量时候，也会形成一些光学重复，optical duplicates。如果变异位点位于这些重复的序列中，可能导致变异频率偏高，因此需要对重复序列进行标记，使得后续变异检测软件可以识别这些重复序列。
我们使用GATK包含的MarkDuplicates功能模块进行重复序列标记。（该功能原属于Picard软件，后被GATK收录）
Copy echo 2 .MarkDuplicates start ` date `
mkdir -p /home/test/output/2.MarkDuplicates/
/gatk/gatk MarkDuplicates \
--java-options '-Xmx2G' \
--INPUT /home/test/output/1.STAR_alignment/SRR5714908.Aligned.sortedByCoord.out.bam \
--OUTPUT /home/test/output/2.MarkDuplicates/SRR5714908.sorted.MarkDup.bam \
--CREATE_INDEX true \
--METRICS_FILE /home/test/output/2.MarkDuplicates/SRR5714908.sorted.MarkDup.bam.metrix \
--VALIDATION_STRINGENCY SILENT
echo 2 .MarkDuplicates end ` date `
#echo 2.MarkDuplicates start `date`
#/BioII/lulab_b/chenyinghui/software/GATK/gatk-4.1.3.0/gatk MarkDuplicates \
#--java-options '-Xmx2G'
#--INPUT /BioII/lulab_b/chenyinghui/project/Docker/SNP/output/1.STAR_alignment/SRR5714908.Aligned.sortedByCoord.out.bam \
#--OUTPUT /BioII/lulab_b/chenyinghui/project/Docker/SNP/output/2.MarkDuplicates/SRR5714908.sorted.MarkDup.bam \
#--CREATE_INDEX true \
#--METRICS_FILE /BioII/lulab_b/chenyinghui/project/Docker/SNP/output/2.MarkDuplicates/SRR5714908.sorted.MarkDup.bam.metrix \
#--VALIDATION_STRINGENCY SILENT
#echo 2.MarkDuplicates end `date` 2c) SplitNCigarReads
SAM/BAM文件的第6列为CIGAR表达式，用来表示该序列各个位置的碱基的比对情况。
由于RNA测序得到的reads中的一部分可能会跨越不同的exon，对于这些reads的CIGAR表达式中会出现N值。使用GATK中的SplitNCigarReads功能将这种reads切分为k+1个reads（K为CIGAR中的N的数量）
Copy echo 3 .SplitNCigarReads start ` date `
/gatk/gatk SplitNCigarReads \
--java-options '-Xmx2G' \
-R /home/test/reference/Homo_sapiens.GRCh38.ch1.fa \
-I /home/test/output/2.MarkDuplicates/SRR5714908.sorted.MarkDup.bam \
-O /home/test/output/3.SplitNCigarReads/SRR5714908.sorted.MarkDup.SplitNCigar.bam
echo 3 .SplitNCigarReads end ` date `
#echo 3.SplitNCigarReads start `date`
#/BioII/lulab_b/chenyinghui/software/GATK/gatk-4.1.3.0/gatk SplitNCigarReads \
#--java-options '-Xmx2G'
#-R /BioII/lulab_b/chenyinghui/project/Docker/SNP/reference/Homo_sapiens.GRCh38.ch1.faa \
#-I /BioII/lulab_b/chenyinghui/project/Docker/SNP/output/2.MarkDuplicates/SRR5714908.sorted.MarkDup.bam \
#-O /BioII/lulab_b/chenyinghui/project/Docker/SNP/output/3.SplitNCigarReads/SRR5714908.sorted.MarkDup.SplitNCigar.bam
#echo 3.SplitNCigarReads end `date` 2d) HaplotypeCaller
该步骤是正式使用GATK进行变异检测的步骤。
Copy echo 4 .HaplotypeCaller start ` date `
mkdir /home/test/output/4.HaplotypeCaller
/gatk/gatk HaplotypeCaller \
--java-options '-Xmx2G' \
--input /home/test/output/3.SplitNCigarReads/SRR5714908.sorted.MarkDup.SplitNCigar.bam \
--output /home/test/output/4.HaplotypeCaller/SRR5714908.raw.vcf.gz \
--reference /home/test/reference/Homo_sapiens.GRCh38.ch1.fa \
-dont-use-soft-clipped-bases \
--standard-min-confidence-threshold-for-calling 20
echo 4 .HaplotypeCaller end ` date `
#echo 4.HaplotypeCaller start `date`
#/BioII/lulab_b/chenyinghui/software/GATK/gatk-4.1.3.0/gatk HaplotypeCaller \
#--java-options '-Xmx2G' \
#--input /BioII/lulab_b/chenyinghui/project/Docker/SNP/output/3.SplitNCigarReads/SRR5714908.sorted.MarkDup.SplitNCigar.bam \
#--output /BioII/lulab_b/chenyinghui/project/Docker/SNP/output/4.HaplotypeCaller/SRR5714908.raw.vcf.gz \
#--reference /BioII/lulab_b/chenyinghui/project/Docker/SNP/reference/Homo_sapiens.GRCh38.ch1.fa \
#-dont-use-soft-clipped-bases \
#--standard-min-confidence-threshold-for-calling 20
#echo 4.HaplotypeCaller end `date` 2e) VariantFiltration
我们可以根据变异的聚集程度、变异的链偏好性、变异的平均质量水平、位点测序深度等指标进行过滤。
Copy echo 5 .VariantFiltration start ` date `
mkdir /home/test/output/5.VariantFiltration
/gatk/gatk VariantFiltration \
--java-options '-Xmx1G' \
--variant /home/test/output/4.HaplotypeCaller/SRR5714908.raw.vcf.gz \
--output /home/test/output/5.VariantFiltration/SRR5714908.filtered.vcf.gz \
--reference /home/test/reference/Homo_sapiens.GRCh38.ch1.fa \
--window 35 \
--cluster 3 \
--filter-name 'FS' \
--filter 'FS > 30.0' \
--filter-name 'QD' \
--filter 'QD < 2.0' \
--filter-name 'DP' \
--filter 'DP < 5'
zcat /home/test/output/5.VariantFiltration/SRR5714908.filtered.vcf.gz | awk -F '\t' '{if ($0 ~ "#" || $7 == "PASS") print $0 }' - >/home/test/output/5.VariantFiltration/SRR5714908.filtered.clean.vcf
echo 5 .VariantFiltration end ` date `
#echo 5.VariantFiltration start `date`
#/BioII/lulab_b/chenyinghui/software/GATK/gatk-4.1.3.0/gatk VariantFiltration \
#--java-options '-Xmx1G' \
#--variant /BioII/lulab_b/chenyinghui/project/Docker/SNP/output/4.HaplotypeCaller/SRR5714908.raw.vcf.gz \
#--output /BioII/lulab_b/chenyinghui/project/Docker/SNP/output/5.VariantFiltration/SRR5714908.filtered.vcf.gz \
#--reference /BioII/lulab_b/chenyinghui/project/Docker/SNP/reference/Homo_sapiens.GRCh38.ch1.fa \
#--window 35 \
#--cluster 3 \
#--filter-name 'FS' \
#--filter 'FS > 30.0' \
#--filter-name 'QD' \
#--filter 'QD < 2.0' \
#--filter-name 'DP' \
#--filter 'DP < 5'
#zcat /BioII/lulab_b/chenyinghui/project/Docker/SNP/output/5.VariantFiltration/SRR5714908.filtered.vcf.gz | awk -F '\t' '{if ($0 ~ "#" || $7 == "PASS") print $0 }' - #>/BioII/lulab_b/chenyinghui/project/Docker/SNP/output/5.VariantFiltration/SRR5714908.filtered.clean.vcf
#echo 5.VariantFiltration end `date` --window 35 \ --cluster3 : （变异的聚集程度）if there are more than 3 variants cluster in 35 bp window, these variants will be filtered.
--filter-name 'FS' \ --filter 'FS > 30.0': （变异的链偏好性（strand bias） ）Phred-scaled p-value （-log10(p-value) ） using Fisher's exact test to detect strand bias. Phred-score closer to 0 means there is a lower chance of there being bias. Higher FS values therefore indicate more bias.
--filter-name 'QD': （变异的平均质量水平）Variant Confidence/Quality by Depth.
--filter-name 'DP'：（位点的测序深度）read depth that cover this site.
--filter-name中的过滤参数的具体含义可以参考GATK生成的VCF文件中开头的注释部分。用户可以根据VCF文件中出现的其他指标对变异进行更多样的过滤筛选。
值得指出的是，满足用户所设置的过滤表达式（如平均质量QD低于2: --filter 'QD < 2.0'）的变异才是我们需要过滤的变异。这些需要被过滤的“不合格”变异仍然会被保留在VCF文件中，但是在VCF第6列 QUAL中会被标注过滤的原因（平均质量QD太低，则标记为QD），通过筛选的、合格的变异位点会被标记PASS。 我们可以用awk等命令去除VCF中不合格变异，保留合格变异。
2f) Annotation
我们可以使用软件以及数据库对得到的变异进行注释，可以获得注释信息包括：变异的位置（位于哪个基因？ 位于exon/intron/UTR ?）、在人群（例如东亚人群）中的频率，临床意义（Pathogenic/Benign）等等。这些注释信息可以帮助研究人员对变异的重要性作出判断。
ANNOVAR将注释分为gene-based annotation、filter-based annotation、Region-based annotation等，在下方的--operation参数中分别对应g与f。
可以用于ANNOVAR注释的公共数据库以及对应的注释类型可以参考ANNOVAR的官网 中的解释。
下面脚本展示的是使用COSMIC、ExAC、ClinVAR数据库对变异进行注释的过程。另外，gnomAD与dbSNP也是常用的注释数据库，推荐大家今后可以自行尝试使用。
Copy echo 6 .Annotation start ` date `
mkdir /home/test/output/6.Annotation
perl /home/test/annovar/table_annovar.pl \
/home/test/output/5.VariantFiltration/SRR5714908.filtered.clean.vcf \
/home/test/annovar/Annovar_database \
--buildver hg38 \
--outfile /home/test/output/6.Annotation/SRR5714908.annotated.variants \
--remove \
--nastring . \
--vcfinput \
--thread 2 \
--protocol ensGene,cosmic70,exac03,clinvar_20190305 \
--operation g,f,f,f
echo 6 .Annotation end ` date `
#echo 6.Annotation start `date`
#perl /BioII/lulab_b/chenyinghui/software/annovar/annovar/table_annovar.pl \
#/BioII/lulab_b/chenyinghui/project/Docker/SNP/output/5.VariantFiltration/SRR5714908.filtered.clean.vcf \ #需要注释的VCF文件
#/BioII/lulab_b/chenyinghui/database/Homo_sapiens/annovar/ \ #注释使用的数据库所处的文件夹路径
#--buildver hg38 \ #比对所用的人类参考基因组的版本，决定了变异坐标是基于何种版本参考基因组
#--outfile /BioII/lulab_b/chenyinghui/project/Docker/SNP/output/6.Annotation/SRR5714908.annotated.variants \ #输出文件路径与文件名前缀
#--remove \ #自动删除中间文件
#--nastring . \ #注释空缺值用.代替
#--vcfinput \ #输入文件为VCF格式
#--thread 2 \
#--protocol ensGene,cosmic70,gnomad211_genome,exac03,clinvar_20190305,avsnp150 \ #使用的数据库
#--operation g,f,f,f,f,f \ #数据库对应的注释类型
#
#echo 6.Annotation end `date` 运行上述脚本之后，最后在输出目录下会有三个文件，后缀名分别是：“.avinput”、“_multianno.txt”、“_multianno.vcf”。
程序先将输入文件从vcf格式转换为“.avinput”格式，然后再对“.avinput”格式保存的变异进行注释。
最后，程序将注释之后的结果的以2种不同格式保存为“_multianno.txt”与“_multianno.vcf”。
“_multianno.txt”中各列以“tab”键分隔，将“_multianno.txt”后缀改为“_multianno.xls”即可在Windows系统中打开。
“_multianno.vcf”则是以vcf格式保存变异以及注释信息。
3) Utilities
3a) ANNOVAR
本示例中使用ANNOVAR进行变异位点信息注释。ANNOVAR是一款优秀的变异注释软件，注释速度快，且可以免费使用。用户可以选择下载公共数据库进行注释，也可以用自己制作的数据库文件（ANNOVAR接受BED/VCF格式）进行注释。
可以使用ANNOVAR提供的Perl脚本下载数据库，如下：
文件已经下载好了。
Copy mkdir /home/test/annovar/Annovar_database
perl /home/test/annovar/annotate_variation.pl \
-buildver hg38 \
-downdb \
-webfrom annovar \
refGene \
/home/test/annovar/Annovar_database 