ANNOVAR注释

使用 ANNOVAR 软件进行基因注释。

创建注释文件

基因注释前需要自己先创建用于注释的文件。

下载 fasta 和 gff3 文件

• 如果没有 gff3 文件，也可以下载 gtf 文件（作者说 Ensembl 的 gff3 有时转换不成功，建议用 gtf ）。
• fasta 文件需要解压缩，gff3 文件不需要

示例如下

mkdir atdb
cd atdb
wget ftp://ftp.ensemblgenomes.org/pub/release-27/plants/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.27.dna.genome.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/release-27/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.27.gtf.gz

处理fasta 和 gff3 文件

结果文件命名就不要带版本了（例如pig10.2），直接用物种名(pig)

这样注释用的 annovar 命令 -buildver 参数就可以直接写物种名。

gff3 文件转格式

gff3ToGenePred 是 UCSC 开发的一款命令行工具，作用是将 GFF3 格式的基因结构注释文件转换为 GenePred 格式（也就是你命令中输出的 pig_refGene.txt）。下载网址：https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/

输入文件可以是压缩或者未压缩格式（我的脚本可以，可以测试一下，如果不支持压缩格式就解压）。

输出文件名称建议就是物种名称加 “_refGene.txt”。如果是 gtf 文件，则改为 gtfToGenePred 脚本。

gff3ToGenePred Sus_scrofa.Sscrofa10.2.89.gff3.gz pig_refGene.txt

作用是将 GFF3 格式的基因结构注释文件转换为 GenePred 格式（也就是你命令中输出的 pig_refGene.txt）。

含义是：读取压缩的 GFF3 文件，解析其中的基因/转录本/外显子/CDS 层级关系，并输出一个每一行代表一个转录本结构的表格文件（基本上都是 mRNA）。这种 GenePred 格式（常被称为 refGene 格式）是 UCSC Genome Browser 和 ANNOVAR 等工具内部常用的基因结构存储格式。

解析输入的 GFF3 文件（支持直接读 .gz压缩包），根据 GFF3 第 9 列的 ID和 Parent标签把基因（gene）、转录本（mRNA/mRNA）、外显子（exon）和 CDS 重新组装起来，最终输出一个以 Tab 分隔的 pig_refGene.txt文件，每行代表一个转录本的完整结构。

它是如何工作的（关键逻辑）

1. 寻找转录本锚点：它主要在 GFF3 中寻找 mRNA、transcript或 RNA类型的行作为每个转录本的“根”。
2. 归集子特征：把具有相同 Parent=转录本ID的 exon和 CDS行归到对应的转录本下，并按坐标排序。
3. 计算边界：自动算出整个转录本的起点/终点（txStart/txEnd），以及编码区 CDS 的起点/终点（cdsStart/cdsEnd）。
4. 命名提取：默认会把 GFF3 里转录本特征的 ID属性作为输出文件第 1 列的 name，如果有 Name属性则作为第 12 列的 geneName（也可以用 -useName参数强制用 Name 作为主名）。

输出格式

输出文件是 16 列、Tab 分隔的文本（扩展版 GenePred），含义如下：

1. name：转录本标识符（通常来自 GFF3 的 ID=）
2. chrom：染色体名称
3. strand：链方向（+或 -）
4. txStart：转录本起始位置（0-based，即 GFF3 的 start - 1）
5. txEnd：转录本终止位置（1-based，同 GFF3 的 end）
6. cdsStart：CDS 编码区起始（0-based，若无 CDS 则与 txStart 相同）
7. cdsEnd：CDS 编码区终止（1-based，若无 CDS 则与 txEnd 相同）
8. exonCount：外显子个数
9. exonStarts：所有外显子起始位置列表（0-based，逗号分隔，末尾有逗号）
10. exonEnds：所有外显子终止位置列表（1-based，逗号分隔，末尾有逗号）
11. score：分数（通常填 0）
12. name2：基因名称/别名（通常来自 GFF3 的 Name=或 gene_name=）
13. cdsStartStat：CDS 起始状态（cmpl完整、incmpl不完整、none无）
14. cdsEndStat：CDS 终止状态（同上）
15. exonFrames：外显子读码框偏移（0,1,2 或 -1）
16. extra：有时会有保留列

这种格式把原本分散在多行的 GFF3 外显子信息“打包”到了一行里，特别适合像 ANNOVAR 这类注释工具快速索引“某个坐标落在哪个转录本的哪根外显子上”，而不用再去实时解析复杂的 Parent 层级关系。

提取fasta序列

这里要改的就是输入输出文件名称，输入是参考基因组的fasta文件和上一步的 pig_refGene.txt ，输出文件是 pig_refGeneMrna.fa 。

这里用的软件就是 annovar 的一个脚本。

perl /mnt/data/zhouziwen/bin/annovar/retrieve_seq_from_fasta.pl --format refGene --seqfile Sus_scrofa.Sscrofa10.2.dna.toplevel.fa pig_refGene.txt --outfile pig_refGeneMrna.fa

这个脚本的作用是根据你指定的“转录本区间文件”，从参考基因组 FASTA 中提取对应的 DNA 序列，并输出为 FASTA 文件。具体到这个例子，就是

✅ 从猪的参考基因组 FASTA 中，根据 pig_refGene.txt里记录的每个转录本（mRNA）的外显子坐标，拼接出每个转录本的 cDNA / mRNA 序列

参数说明：

参数	含义
--format refGene	输入文件是 refGene 格式（UCSC GenePred）
--seqfile ...fa	单文件参考基因组 FASTA（不是按 chr 拆分）
pig_refGene.txt	转录本区间定义文件
--outfile	输出的 mRNA 序列 FASTA

脚本会：

• 按外显子坐标从 FASTA 中取序列
• 正链直接拼接
• 负链反向互补
• 最终得到 mRNA 序列

输出文件：FASTA 格式，每个转录本一条序列

>NM_001001 Comment: this sequence (leftmost exon at chr1:1000) is generated by ANNOVAR on ...
ATGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGC
TAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGC

内容	含义
>NM_001001	转录本 ID
序列	拼接后的 mRNA / cDNA 序列

检查两个文件的转录本数目是否相同（二者必须相同）

wc -l pig_refGene.txt

grep ">" pig_refGeneMrna.fa | wc -l

提取基因ID与名称的对应关系

从 gff3 文件中提取基因ID与名称的对应关系，这是用于 annovar 注释完之后的分析。因为查找基因功能的时候往往需要使用的是基因名称，因此在 annovar 注释到基因之后需要转换成基因名称。

原理：自动跳过注释行并仅提取类型为 gene的记录，从第九列属性中解析出 ID和 Name，自动去除 Ensembl 风格的 gene:前缀，最终生成一个两列（基因ID \t基因名）的映射文件，并在结束时打印提取到的基因总数。

python pick_gene_name.py Sus_scrofa.Sscrofa10.2.89.chr.gff3 geneIDmap.txt

annovar 注释

官方文档

github 文档 和 它的官网的文档应该是一致的，只是 github 可以下载下来。

github文档网址：https://github.com/WGLab/doc-ANNOVAR

官网文档网址：http://annovar.openbioinformatics.org/en/latest/user-guide/download/

下载软件

从 here 下载，需要注册

脚本说明

• table_annovar.pl : 目前的核心脚本，可以一次调用多个数据库（如 gnomAD、dbSNP、dbNSFP）进行多功能注释，基于参数 -operation详解： g= 基因注释 (gene-based) r= 区域注释 (region-based) f= 过滤注释 (filter-based) 。其输入文件就是 VCF 文件按

  但是例子给的是人的，动物的不太会，目前我还是用专门做基因注释的 annotate_variation.pl 。

输入文件

这里主要介绍 ANNOVAR 可接受的文件输入格式，以及如何利用 convert2annovar.pl工具将其他格式的变异文件转换为 ANNOVAR 可识别的格式。

输入文件格式概述

• VCF 文件：table_annovar.pl程序可以直接接受 VCF 文件作为输入（需加 -vcfinput参数）。目前 VCF 已成为大多数研究者使用的黄金标准格式。
• ANNOVAR 输入格式 (.avinput)：annotate_variation.pl等程序需要一个简单的基于文本的格式（称为 ANNOVAR input format）。
  • 文件中每一行对应一个变异位点。
  • 每一行的前 5 列（空格或 Tab 分隔）分别代表：染色体 (Chromosome)、起始位置 (Start)、终止位置 (End)、参考碱基 (Reference Allele)、变异碱基 (Alternative Allele)。
  • 第 5 列之后可以添加任意额外的列，这些列会在输出中被原样打印。
  • 如果参考碱基信息未知，可用 0填充。
  • 插入或缺失可用 -表示空碱基（例如 TC到 -代表缺失，-到 AT代表插入）。
  • 默认使用 1-based 坐标系统。
  • 如果 ANNOVAR 遇到无效输入行，会将其写入 <outfile>.invalid_input文件。

示例 (ex1.avinput)：

1   948921   948921   T   C   comments: rs15842, a SNP in 5' UTR of ISG15
1   13211293 13211294 TC  -   comments: rs59770105, a 2-bp deletion
1   11403596 11403596 -   AT  comments: rs35561142, a 2-bp insertion
13  20797176 21105944 0   -   comments: a 342kb deletion encompassing GJB6

输入文件要求

• 第5列（ALT）不能同时包含碱基和 -（这个我测试过了）
• 第4列 （REF）也不能同时包含碱基和 - 。

格式转换工具 convert2annovar.pl

该脚本可将其他“基因型调用”格式转换为 ANNOVAR 的 .avinput格式。目前支持 Samtools pileup、Illumina CASAVA、SOLiD GFF、Complete Genomics、SOAPsnp、MAQ 和 VCF 等格式。如今 VCF 最常用，其他格式已较少使用。

转换 VCF 文件

官方案例使用 -format vcf4参数（实际我使用 vcf4old 参数）。

输入VCF文件可以是 vcf.gz 压缩格式。

# 基本转换（默认只处理第一个样本）
perl convert2annovar.pl -format vcf4 example/ex2.vcf > ex2.avinput

# 指定输出文件
perl convert2annovar.pl -format vcf4 example/ex2.vcf -outfile ex2.avinput

# 转换所有样本（生成多个文件，如 ex2.NA00001.avinput, ex2.NA00002.avinput）
perl convert2annovar.pl -format vcf4 example/ex2.vcf -outfile ex2 -allsample

# 包含 VCF 的 INFO 信息
perl convert2annovar.pl -format vcf4 example/ex2.vcf -outfile ex2 -includeinfo

# 包含杂合性、质量、深度信息
perl convert2annovar.pl -format vcf4 example/ex2.vcf -outfile ex2 -withzyg

# 打印基于文件内所有样本的等位基因频率（-withfreq）
perl convert2annovar.pl -format vcf4 example/ex2.vcf -outfile ex2 -withfreq -includeinfo

# 保留 VCF 头信息
perl convert2annovar.pl -format vcf4 example/ex2.vcf -outfile ex2 -comment

注意：如果某个位点有多个可变等位基因（Multiple Alleles），输出会产生多行，每行对应一个变异。

vcf4 & vcf4old

我看之前的脚本有的使用 -format vcf4old ，这里和 -format vcf4 的区别就是输入文件 VCF 版本问题。

✅ -format vcf4用于“新版 VCF（4.x）”，支持多样本、多 ALT、基因型字段

✅ -format vcf4old用于“老版 VCF（4.0 之前 / 简化版）”，只处理单样本、单列 ALT

在 misc/faq.md 文件中，作者进行了解释。

1. vcf4old：过时的全量转换模式

该模式保留了 2013 年之前的旧逻辑。它的核心特点是**“见行就转”**：无论 VCF 中的特定样本在该位点是否有突变（即使是 0/0纯合参考），只要 VCF 文件中有这一行记录，它就会被转换为输出文件的一行。这种模式在现代多样本联合 Calling（Joint Calling）场景下会产生大量无意义的冗余数据，因此已被视为 Obsolete（过时），不建议新用户使用。

2. vcf4：智能的现代转换模式*

这是目前的默认推荐模式。它理解“变异”的定义：对于一个特定样本，只有当基因型不是 0/0时，才算作变异。默认情况下，它只处理 VCF 中的第一个样本，并且只输出该样本真正发生突变的位点，这大大减少了无效数据。此外，vcf4支持通过参数扩展功能，能够优雅地处理多样本数据。

3. 参数组合与行为对比

下表详细列出了两种模式在处理多样本 VCF 时的行为差异：

模式 / 参数	行为逻辑	适用场景
vcf4(默认)	仅处理第一个样本。仅输出该样本基因型非 0/0的位点。	单样本 VCF，或只关心 VCF 文件中的第一个样本。
vcf4 -allsample	处理所有样本。为每个样本生成一个独立的 .avinput输出文件。	多样本 VCF，需要将每个人的变异数据拆分开进行独立分析。
vcf4 -allsample -withfreq	全量保留 + 计算频率。保留 VCF 中的所有行（包括 0/0），并在输出中新增一列，计算该变异在所有样本中的等位基因频率。	群体遗传学分析，需要知道每个位点在整个队列中的频率。
vcf4old	全量保留。无论样本基因型如何，只要 VCF 有记录就输出。	极不推荐。仅用于兼容非常古老、依赖旧输出格式的特定流程。

总结: vcf4 会丢失无变异位点，因此还是使用 vcf4old

基因注释

这是最重要的一章：网址为 https://annovar.openbioinformatics.org/en/latest/user-guide/gene/，下载的 github 中为 gene.md (这个貌似内容更多)

本文档主要介绍 ANNOVAR 最核心的功能之一：基于基因的变异注释。即判断一个变异（SNP/Indel）落在基因组的什么位置（外显子、内含子等），以及是否引起氨基酸改变，并详细说明其输出文件和各类参数。

注释前准备

在进行基因注释前，必须将对应的基因定义文件（如 refGene.txt）和转录本 FASTA 文件（如 refGeneMrna.fa）下载到 humandb/ 目录中（你自己的某个目录中）。

命令及参数说明

文档压根没有一个完整的参数说明和解释。

我一般使用的注释命令如下

annotate_variation.pl -buildver LS -geneanno -outfile Anno SNP.annovar /work/home/ach0yyfeio/reference/LS_na

命令拆解：

annotate_variation.pl \
  -buildver LS \
  -geneanno \
  -outfile Anno \
  SNP.annovar \
  /work/home/ach0yyfeio/reference/LS_na

参数解释：这里 -buildver 后面接的就是注释文件的前缀，输入文件为 SNP.annovar ，输出前缀为 Anno ，最后是注释文件的文件夹的路径。

参数	含义
annotate_variation.pl	ANNOVAR 的主注释脚本
-buildver LS	指定参考基因组版本为 LS （ANNOVAR 会在数据库目录中寻找 LS_refGene.txt、LS_refGeneMrna.fa等文件）
-geneanno	只做基因注释 （不跑区域注释、过滤注释等）
-outfile Anno	输出文件前缀 会生成 Anno.variant_function和 Anno.exonic_variant_function
SNP.annovar	输入文件 ANNOVAR 专用格式（.avinput），包含染色体、位置、Ref、Alt
/work/home/ach0yyfeio/reference/LS_na	ANNOVAR 数据库目录 必须包含 LS 物种的基因注释文件

其他参数与说明（除了第一个，其他没看懂）

• –separate：默认只输出最高优先级的注释（如 exonic 覆盖了 intronic，只报 exonic）。加此参数会列出该变异所有符合的功能类别（每行一个）。
• –hgvs：让 cDNA 和蛋白变化注释符合 HGVS 命名标准（如 c.122C>T 而非 c.C122T）。
• –transcript_function：默认第 2 列输出基因名，加此参数可输出转录本名（如 NM_ 或 ENST_ 编号）。
• -exonsort：当同一变异影响多个转录本时，强制按外显子序号排序输出，保证结果绝对一致（避免随机排序）。
• –aamatrixfile：可在错义变异后附加氨基酸替代矩阵分数（如 Grantham score）。

基因定义系统（Database / -dbtype）

ANNOVAR 支持多种基因集，只需下载对应文件并使用 -dbtype 指定（不懂其他的基因集如何使用，但是自己生成的注释文件就是使用默认的 refGene，所以这个应该不设定）：

• refGene：RefSeq 基因（经典，较少基因，质量高）。
• refGeneWithVer：带版本号的 RefSeq（如 NM_005101.4），推荐用于精确追踪。
• knownGene：UCSC Known Gene（较多，包含部分预测基因）。
• ensGene：Ensembl/GENCODE 基因（最全面，包含很多 lncRNA 等，hg38 建议使用 GENCODE）。
• ccdsGene：CCDS 基因（严格保守的核心基因集，输出只有 CCDS ID，需自行转换基因名）。

注意：不同基因集对同一变异的注释可能不同（如某个变异在一个系统里是 synonymous，在另一个里因不同转录本可能是 nonsynonymous），ANNOVAR 允许用户灵活切换对比。

输出文件详解（核心）

当运行基因注释（无论是通过 table_annovar.pl还是底层命令）时，主要生成两个关键文件：

*variant_function（所有变异的注释）

• 格式：在原始输入行前增加 2 列。

• 第 1 列（变异类型）：ANNOVAR 将基因组划分为 9 种区域，并遵循严格的优先级（Precedence）（数字越小越优先）。当一个变异同时属于多个功能区时，ANNOVAR 只报“优先级最高”的那个，而不会全都报。

  exonic> splicing> ncRNA> UTR5/UTR3> intronic> upstream/downstream> intergenic。

  优先级	值 (Value)	解释 (Explanation)	序列本体 (Sequence Ontology)
  1	exonic(外显子)	变异位于编码区	exon_variant(SO:0001791)
  1	splicing(剪接)	变异位于剪接位点 2bp 范围内(可通过 -splicing_threshold修改范围)	splicing_variant(SO:0001568)
  2	ncRNA(非编码RNA)	变异重叠的基因定义中无编码注释的转录本	non_coding_transcript_variant(SO:0001619)
  3	UTR5(5’非翻译区)	变异重叠 5’ 端非翻译区	5_prime_UTR_variant(SO:0001623)
  3	UTR3(3’非翻译区)	变异重叠 3’ 端非翻译区	3_prime_UTR_variant(SO:0001624)
  4	intronic(内含子)	变异位于内含子中	intron_variant(SO:0001627)
  5	upstream(上游)	变异位于转录起始位点上游 1kb 区域	upstream_gene_variant(SO:0001631)
  5	downstream(下游)	变异位于转录终止位点下游 1kb 区域(可通过 -neargene修改范围)	downstream_gene_variant(SO:0001632)
  6	intergenic(基因间)	变异位于基因间区域	intergenic_variant(SO:0001628)

• 第 2 列（基因信息）：

  • 若落在基因区域：给出基因名（多个基因用逗号分隔，如 exonic ATG16L1,BRAF）。
  • 若落在基因间区：给出两侧最近的基因名和距离（如 intergenic UBIAD1(dist=43968),DISP3(dist=135699)）。

• 修改优先级：可用 -precedence intronic,utr5,utr3等参数自定义优先级顺序。

*exonic_variant_function（仅外显子变异的详细注释）

• 格式：仅包含外显子上的变异，在原始输入行前增加 3 列。

• 第 1 列：该变异在原始输入文件中的行号（Line #）。

• 第 2 列（功能后果）：详细描述对蛋白的影响，包括：

  • nonsynonymous SNV（错义）、synonymous SNV（同义）
  • frameshift insertion/deletion（移码插入/缺失）
  • nonframeshift insertion/deletion（非移码插入/缺失）
  • stopgain（无义，引入终止密码子）、stoploss（终止密码子丢失）
  • unknown（未知，通常因转录本缺乏完整 ORF 或存在提前终止密码子）。

  优先级如下，默认情况下是唯一的，并且严格按照“蛋白层面优先级”只输出最高的那一个。

  优先级 (Precedence)	注释类型 (Annotation)	解释 (Explanation)	序列本体 (Sequence Ontology)
  1	移码插入 (frameshift insertion)	插入一个或多个核苷酸，导致蛋白质编码序列发生移码变化。	frameshift_elongation(SO:0001909)
  2	移码缺失 (frameshift deletion)	缺失一个或多个核苷酸，导致蛋白质编码序列发生移码变化。	frameshift_truncation(SO:0001910)
  3	移码块替换 (frameshift block substitution)	替换一个或多个核苷酸的块，导致蛋白质编码序列发生移码变化。	frameshift_variant(SO:0001589)
  4	无义突变 (stopgain)	单核苷酸变异（SNV）、移码插入/缺失、非移码插入/缺失或块替换，导致在变异位点立即产生终止密码子。对于移码突变，在变异位点下游产生的终止密码子不计入“stopgain”！	stop_gained(SO:0001587)
  5	终止密码子丢失 (stoploss)	单核苷酸变异（SNV）、移码插入/缺失、非移码插入/缺失或块替换，导致在变异位点立即消除终止密码子。	stop_lost(SO:0001578)
  6	非移码插入 (nonframeshift insertion)	插入3个或3的倍数的核苷酸，不导致蛋白质编码序列发生移码变化。	inframe_insertion(SO:0001821)
  7	非移码缺失 (nonframeshift deletion)	缺失3个或3的倍数的核苷酸，不导致蛋白质编码序列发生移码变化。	inframe_deletion(SO:0001822)
  8	非移码块替换 (nonframeshift block substitution)	替换一个或多个核苷酸的块，不导致蛋白质编码序列发生移码变化。	inframe_variant(SO:0001650)
  9	非同义 SNV (nonsynonymous SNV)	单个核苷酸的改变导致氨基酸发生改变。	missense_variant(SO:0001583)
  10	同义 SNV (synonymous SNV)	单个核苷酸的改变不导致氨基酸发生改变。	synonymous_variant(SO:0001819)
  11	未知 (unknown)	功能未知（由于数据库文件中基因结构定义的多种错误）。	sequence_variant(SO:0001060)

• 第 3 列（细节）：包含 基因名:转录本ID:外显子号:cDNA变化:氨基酸变化。

  • 例：IL23R:NM_144701:exon9:c.G1142A:p.R381Q。
  • 若一个基因有多个转录本（可变剪切），会用逗号分隔列出所有情况。