6.5.SNV/INDEL
本章介绍如何通过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）的软件，拥有较高的引用率（已有上万次引用）。
GATK Forum （GATK开发人员与在该论坛回答用户疑问）
参考文献：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
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 bash2a) Alignment
GATK要求输入的SAM/BAM文件中有Read Group信息，因此我们需要在--outSAMattrRGline中填写相应的Read Group信息，其中ID为必填项。Read Group格式详见SAM格式
2b) MarkDuplicates
在建库的PCR过程中会形成一些重复的DNA片段序列，这些重复序列被称为PCR duplicates。另外，在测序仪进行光学测量时候，也会形成一些光学重复，optical duplicates。如果变异位点位于这些重复的序列中，可能导致变异频率偏高，因此需要对重复序列进行标记，使得后续变异检测软件可以识别这些重复序列。
我们使用GATK包含的MarkDuplicates功能模块进行重复序列标记。（该功能原属于Picard软件，后被GATK收录）
2c) SplitNCigarReads
SAM/BAM文件的第6列为CIGAR表达式，用来表示该序列各个位置的碱基的比对情况。
由于RNA测序得到的reads中的一部分可能会跨越不同的exon，对于这些reads的CIGAR表达式中会出现N值。使用GATK中的SplitNCigarReads功能将这种reads切分为k+1个reads（K为CIGAR中的N的数量）
2d) HaplotypeCaller
该步骤是正式使用GATK进行变异检测的步骤。
2e) VariantFiltration
我们可以根据变异的聚集程度、变异的链偏好性、变异的平均质量水平、位点测序深度等指标进行过滤。
--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也是常用的注释数据库，推荐大家今后可以自行尝试使用。
运行上述脚本之后，最后在输出目录下会有三个文件，后缀名分别是：“.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脚本下载数据库，如下：
文件已经下载好了。
参考文献： Wang K, et al. ANNOVAR: Functional annotation of genetic variants from next-generation sequencing data Nucleic Acids Research. 2010. 38:e164.
Last updated
Was this helpful?