BSA分析（六）变异结果注释

❝

「因业务拓展，想组建一个数据分析团队（目前已有RNA-Seq、Chip-Seq、重测序与群体遗传、基因家族、比较基因组、宏基因组、微生物多样性16s/18s/ITS方向专业人员），欢迎有各种数据分析基础的朋友加入我们！——Bioinfor 生信云」

❞

变异注释

annovar 和 SnpEff 都是比较常用的 SNP 和 INDEL 注释软件。

annovar

ANNOVAR是一个高效的Perl 命令注释工具，可实现对来自下一代测序 （NGS） 数据的遗传变异进行高通量功能注释和过滤，支持包括VCF在内的多种输入和输出文件格式。

主要包含三种不同的注释方法：gene-based, region-based 和filter-based。基于基因的注释（Gene-based Annotation）揭示variant与已知基因直接的关系以及对其产生的功能性影响，基于区域的注释（Region-based Annotation）揭示variant 与不同基因组特定段的关系，例如：它是否落在转录因子结合区域等，基于筛选的注释（Filter-based Annotation）则分析变异位点是否位于指定的数据库中，比如dbSNP, 1000G，ESP 6500等数据库，计算SIFT/PolyPhen/LRT/MutationTaster/MutationAssessor等

• 数据：基因组文件（genome）、基因组注释文件（gtf或gff）、VCF文件、INDEL文件
• 软件：annovar,gtfToGenePred（gtf 文件做输入）,gff3ToGenePred(gff3 文件做输入)

# step1 转gtf或gff成refGene格式并提取转录本的序列
genome.gtf
genome.fasta

## 转gtf成refGene格式
gtfToGenePred -genePredExt genome.gtf genome_refGene.txt

##如果提供的为gff文件，使用下面命令
gff3ToGenePred genome.gff genome_refGene.txt

## 提取序列文件

/home/software/annovar/retrieve_seq_from_fasta.pl --format refGene --seqfile genome.fasta genome_refGene.txt --out genome_refGeneMrna.fa

# step2 annovar注释

## 转vcf成表格格式（软件需求）
/home/software/annovar/convert2annovar.pl \
-format vcf4 \ #输入文件格式
-allsample \ #allsample模式
-withfreq \ #输出alt频率
genome.vcf \ #输入文件
> genome.annovar.input #输出文件

## 进行变异注释 如果需要所有信息，添加-separate
/home/software/annovar/annotate_variation.pl  \
-geneanno \ #基于基因进行注释
--neargene 2000 \ #确定gene上下游范围
-buildver genome \ #基因组数据库名称
-dbtype refGene \ #数据库类型
-outfile all.anno \ #输出文件前缀
-exonsort \ # 结果按exon进行排序
genome.annovar.input \ # 输入文件
./ # 输入数据库路径

# 输出文件为：
## all.anno.exonic_variant_function
## all.anno.variant_function

# 基本统计
cut -f 1 all.anno.variant_function |sort | uniq -c > all.anno.variant_function.stat
cut -f 2 all.anno.exonic_variant_function |sort |uniq -c > all.anno.exonic_variant_function.stat

「variant_function注释结果，优先级排序为 exonic = splicing > ncRNA > UTR5/UTR3 > intron > upstream/downstream > intergenic」

「INDEL注释直接将脚本内相应snp文件替换成indel文件即可」

SnpEff

SnpEff 软件使用时需要进行基因组数据库构建，这个软件已经预先构建好了超过 2500 个物种的数据库, 如果使用的基因组版本在预先构建的数据库里存在，可以直接下载使用。

如果研究的物种没有预先建好的数据库，可以自己手动来构建。

# 注释
java -jar snpEff.jar -c ./snpEff.config \ #指定数据库
-ud 2000 \ #指定上下游范围
-csvStats test.csv \ #生成表格的统计文件
-htmlStats test.html \ #生成网页的统计文件
-o vcf \ #输出文件vcf格式
ric \ #数据库名称
genome.vcf > genome.ann.vcf

## 输出文件 genome.ann.vcf

snpEff 变 异 类 型 说 明 https://pcingola.github.io/SnpEff/se_inputoutput/#effect-prediction-details

对区域进行注释

BSA 分析得到的结果是基因组区域，需要查看区域内有哪些基因，可以使用 bedtools 进行查看。

# 准备gff或gtf文件，以及需要查看的区域bed文件
genome.gtf
genome_region.bed

# 进行区域内基因查找
bedtools intersect -wo -a genome_region.bed -b genome.gtf | awk '$6=="gene"' > genome_region.anno

欢迎加入学习交流群