群体遗传研究——变异注释–snpEff和ANNOVAR

ht • 2023-11-17 09:11 • 杂文

snpEff

snpEff4.0软件之后默认下载安装所有动物的库（ensembl的库）

查询对应物种的数据库

java -jar /home/sll/software/snpEff/snpEff.jar databases | grep "Bos"

但是ensembl的库有个缺点就是注释出来的基因是ensemblid，我们也不认识它是啥，因此，需要我们下载NCBI参考基因组进行建库

下载参考基因组和gtf注释文件就不多说了

1、安装，构建数据库

#会产生一个snpEff目录 所有的程序都在这里面
下载 https://pcingola.github.io/SnpEff/
unzip snpEff_latest_core.zip

2、配置自己的基因组和注释文件

以牛（cattle）参考基因组为例：打开snpEFF文件夹下的snpEff.contig，在Third party databases下面增加新的物种信息

#-------------------------------------------------------------------------------
# Third party databases    （原来的）
#-------------------------------------------------------------------------------
#sheep  genome, version cattlev1.2    (新加的)
cattlev1.2.genome : cattle     （新加的）
# Databases for human GRCh38 (hg38)    （原来的）

3.注释文件为 gff 格式

# 在SnpEff目录下，创建目录 data, data/cattlev1.2, data/genomes
mkdir data && cd data
mkdir cattlev1.2
mkdir genomes
# 将 gff3 (.gff) 文件放入sheepv1.0目录下，并改名为 genes.gff
# 将基因组序列文件（.fasta）放入 genomes 目录下，并改名为 cattlev1.2.fa
nohup  java -jar /home/sll/software/snpEff.jar build -gff3 /
                                                     -v cattlev1.2 & # 在 snpEff 目录下，执行命令

4.注释文件为 gtf 格式

# 如果注释文件为.gtf，可参考 gff 中的步骤，要将注释文件重新命名为 genes.gtf
nohup java -jar /home/sll/software/snpEff.jar build -gtf22 /
                                                    -v cattlev1.2 /
                                                    -noCheckCds /
                                                    -noCheckProtein &

# -noCheckCds 不检查CDS区
# -noCheckProtein 不检查蛋白质
#如果不加上述两个参数，snpEff软件默认检查，则需要我们提供物种相应的CDS及蛋白质的fa文件

5.开始注释 (提取的要注释的文件中染色体应为NC格式且文件以tab键分隔，提取三列即可：染色体号、起始位置、终止位置，若为SNP位点则第三列用1代替，且位于snpEff目录下操作命令）

（1）若文件为按照窗口计算Fst后输出的文件，则提取为bed文件格式，bed格式(取染色体号，起始位置和结束位置和Fst 的值)

# 计算Fst后输出表格的格式: 1.CHROM（染色体号） 2.BIN_START （开头）3.BIN_END（结尾） 4.N_VARIANTS（SNP数） 5.WEIGHTED_FST（加权Fst） 6.MEAN_FST（平均Fst）
awk  -F '/t'  '{print $1"/t"$2"/t"$3"/t"$5}'   XXX.windowed.weir.sorted.1%.fst   > XXX.windowed.weir.sorted.1%.bed
#如果bed文件中染色体为1.2.3格式，需替换染色体
nohup java -Xmx4g -jar /home/sll/software/snpEff.jar cattlev1.2 -i bed XXX.windowed.weir.sorted.1%.bed-chr.bed   > XXX.windowed.weir.sorted.1%.bed-chr.anno.bed &

（2）若文件为按照位点计算（例如：fst按位点计算、重测序SNP数据、重测序INDEL数据）但重测序使用中，得到的注释文件不理想

# 用awk提取fst中的snp 位置信息，并用vcftools对应回vcf文件中
nohup awk   '{print $1,$4,$4+1}'  test.map > test-awk.bed  &
# 注释
nohup java -Xmx4g -jar /home/sll/software/snpEff.jar cattlev1.2 -i bed test-awk.bed   > test-awk.anno.bed &

（3）若文件为vcf文件 (重测序SNP数据、重测序INDEL数据)

比如只关注CDS中的注释信息，不考虑上游、下游、UTR、基因间区等信息
可以选择以下参数简化输出
-no-downstream
-no-upstream
-no-utr
-no-intergenic
-no-intron

nohup java -jar snpEff.jar ann -no-utr /
                               -no-downstream /
                               -no-upstream /
                               -no-intergenic sheepv1.0 test.vcf.gz /
                               > snpeff.vcf &
# 此处可将变异信息删除做一个小的vcf文件，只保留到vcf的info信息即可，则最后出现的文件小一些

ANNOVAR

1、下载地址：

http://www.openbioinformatics.org/annovar/download/0wgxR2rIVP/annovar.latest.tar.gz

2、解压：

tar -zxvf annovar.latest.tar.gz

3、构建注释数据库：

# 下载注释gtf和fa文件:建议在Ensembl或UCSC下载

4、下载安装gtfToGenePred工具

wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64.v369/gtfToGenePred

chmod +x gtfToGenePred

5 用 gtfToGenePred 工具将 GTF file 转换 GenePred file

/home/sll/software/gtfToGenePred -genePredExt GCF_002263795.1_ARS-UCD1.2_genomic.gtf ARS-UCD1.2_refGene.txt

NCBI下载的第一次可能会能报错：使用下面脚本解决
位置：/home/sll/software/annovar/cattle
replacegtf.py （记得改内容）

/home/sll/software/gtfToGenePred -genePredExt GCF_002263795.1_ARS-UCD1.2_genomic_replace.gtf ARS-UCD1.2_refGene.txt

6 生成转录组信息文件

perl /home/sll/software/annovar/retrieve_seq_from_fasta.pl --format refGene /
                                                           --seqfile GCF_002263795.1_ARS-UCD1.2_genomic.fna ARS-UCD1.2_refGene.txt /
                                                           -outfile ARS-UCD1.2_refGeneMrna.fa

# -format指定gene definition file格式
# -seqfile 指定基因组序列文件名称
# -outfile 指定输出mRNA序列文件的名称

建库完成

注释示例：

1、vcf转为适宜的格式

perl /home/sll/software/annovar/convert2annovar.pl -format vcf4old QC.JPBC-geno005-maf003.bed.vcf /
                                                   -outfile QC.JPBC-geno005-maf003.avinput

#关于-format vcf4,并没有保留全部位点：
#WARNING to old ANNOVAR users: this program no longer does line-to-line conversion for multi-sample VCF files. If you want to include all variants in output, use '-format vcf4old' or use '-format vcf4 -allsample -withfreq' instead.

2、annotate_variation注释

NC染色体
第一种，使用annotate_variation.pl ，三种类型分开：

perl /home/sll/software/annovar/annotate_variation.pl -out jpbc.anno /
                                                      -dbtype refGene /
                                                      -buildver cattle1.2 QC.JPBC-geno005-maf003.avinput /
                                                      /home/sll/software/annovar/Bos/

#目前以上自建库只能基于gene注释
# -geneanno  表示使用基于基因的注释 一般是默认的
# -dbtype refGene  表示使用"refGene"类型的数据库
# -out jpbc.anno  表示输出以jpbc.anno为前缀的结果文件

基于基因的注释运行后会产生三个文件：variant_function，exonic_variant_function，.log

variant_function 注释所有变异所在基因及位置。第1列为变异所在的类型，如外显子等，第2列是对应的基因名(若有多个基因名用，隔开)，第3-7列为输入的那5列主要信息，剩余为注释信息。如果变异找到多种注释，ANNOVAR将根据优先权重进行比较，取最优的表示，可以使用-seperate参数列出该变异所有注释。

image.png
exonic_variant_function 详细注释外显子区域的变异功能、类型、氨基酸改变等。第1列为.variant_function文件中该变异所在行号，第2列为变异功能性后果，如外显子改变导致的氨基酸变化，阅读框移码，无义突变，终止突变等，第3列包括基因名称、转录识别标志和相应的转录本的序列变化，第4-9列为输入文件内容。

image.png

查看有多少种gene variant

# exonic_variant_function文件
cat x.exonic_variant_function |cut -f2|cut -d" " -f1|sort |uniq -c |sort -nr

# variant_function文件
cat x.variant_function |cut -f1|cut -d" " -f1|sort |uniq -c |sort -nr

# 基于基因
annotate_variation.pl -geneanno /
                      -dbtype refGene /
                      -buildver hg19 /
                      example/ex1.avinput /
                      humandb/
#基于区域
annotate_variation.pl -regionanno /
                      -dbtype cytoBand -buildver hg19 /
                      example/ex1.avinput /
                      humandb/ 
#基于筛选
annotate_variation.pl -filter /
                      -dbtype exac03 /
                      -buildver hg19 /
                      example/ex1.avinput /
                      humandb/

第二种，基于table_annovar.pl，直接注释三种类型（如果是自建库，则推荐使用第一种，因为只能基于基因注释）：

table_annovar.pl是ANNOVAR多个脚本的封装，可以一次性完成三种类型的注释

perl /home/sll/software/annovar/table_annovar.pl QC.JPBC-geno005-maf003.avinput /
                                             /home/sll/software/annovar/cattle/ /
                                             -buildver ARS-UCD1.2 /
                                             -out jpbc.anno /
                                             -remove /
                                             -protocol refGene /
                                             -operation g,r,f /
                                             -nastring NA /
                                             -csvout

# -buildver ARS-UCD1.2 表示使用的参考基因组版本为ARS-UCD1.2
# -out jpbc.anno 指定输出文件前缀为jpbc.anno
# -remove 表示删除中间文件
# -protocol 后跟注释来源数据库名称，每个protocal名称或注释类型之间只有一个逗号，并且没有空白
# -operation 后跟指定的注释类型，和protocol指定的数据库顺序是一致的，g代表gene-based、r代表region-based、f代表filter-based
# -nastring NA 表示用NA替代缺省值
# -csvout 表示最后输出.csv文件

Ensembl基因ID转为NCBI基因ID

1. 写如下文件，储存为GeneName2Symble.pl文件。

除了前三个文件名需要更改为自己的，分别为注释文件，输入文件名，输出文件名，其他的都一致即可。将该.pl文件与里面的注释文件，文本文件放在同一文件夹内。
链接：https://www.jianshu.com/p/c5010d5fc997

##用法perl GeneName2Symble.pl
use strict;
use warnings;

my $gtfFile="Homo_sapiens.GRCh38.98.chr.gtf";
my $expFile="tpm.txt";
my $outFile="tpm_symble.txt";

my %hash=();
open(RF,"$gtfFile") or die $!;
while(my $line=)
{
    chomp($line);
    if($line=~/gene_id /"(.+?)/"/;.+gene_name "(.+?)"/;.+gene_biotype /"(.+?)/"/;/)
    {
        $hash{$1}=$2;
    }
}
close(RF);

open(RF,"$expFile") or die $!;
open(WF,">$outFile") or die $!;
while(my $line=)
{
    if($.==1)
    {
        print WF $line;
        next;
    }
    chomp($line);
    my @arr=split(//t/,$line);
    $arr[0]=~s/(.+)/..+/$1/g;
    if(exists $hash{$arr[0]})
    {
        $arr[0]=$hash{$arr[0]};
        print WF join("/t",@arr) . "/n";
    }
}
close(WF); 
close(RF)

2、输入perl GeneName2Symble.pl运行即可

版权声明：
作者：ht
链接：https://www.techfm.club/p/86136.html
来源：TechFM
文章版权归作者所有，未经允许请勿转载。

THE END

download GitHub

二维码

汉字美学读书摘抄

< <上一篇

《<史记>研读》146：晋世家（十二）

下一篇>>

搜索内容

群体遗传研究——变异注释–snpEff和ANNOVAR

snpEff

snpEff4.0软件之后默认下载安装所有动物的库（ensembl的库）

1、安装，构建数据库

2、配置自己的基因组和注释文件

3.注释文件为 gff 格式

4.注释文件为 gtf 格式

5.开始注释 (提取的要注释的文件中染色体应为NC格式且文件以tab键分隔，提取三列即可：染色体号、起始位置、终止位置，若为SNP位点则第三列用1代替，且位于snpEff目录下操作命令）

ANNOVAR

1、下载地址：

2、解压：

3、构建注释数据库：

4、下载安装gtfToGenePred工具

5 用 gtfToGenePred 工具将 GTF file 转换 GenePred file

6 生成转录组信息文件

注释示例：

1、vcf转为适宜的格式

2、annotate_variation注释

查看有多少种gene variant

Ensembl基因ID转为NCBI基因ID

1. 写如下文件，储存为GeneName2Symble.pl文件。

2、输入perl GeneName2Symble.pl运行即可

取消回复

共有 0 条评论

Ads

群体遗传研究——变异注释–snpEff和ANNOVAR

snpEff

snpEff4.0软件之后默认下载安装所有动物的库（ensembl的库）

1、安装 ，构建数据库

2、配置自己的基因组和注释文件

3.注释文件为 gff 格式

4.注释文件为 gtf 格式

5.开始注释 (提取的要注释的文件中染色体应为NC格式且文件以tab键分隔，提取三列即可：染色体号、起始位置、终止位置，若为SNP位点则第三列用1代替，且位于snpEff目录下操作命令）

ANNOVAR

1、下载地址：

2、解压：

3、构建注释数据库：

4、下载安装gtfToGenePred工具

5 用 gtfToGenePred 工具将 GTF file 转换 GenePred file

6 生成转录组信息文件

注释示例：

1、vcf转为适宜的格式

2、annotate_variation注释

查看有多少种gene variant

Ensembl基因ID转为NCBI基因ID

1. 写如下文件，储存为GeneName2Symble.pl文件。

2、输入perl GeneName2Symble.pl运行即可

取消回复

共有 0 条评论

Ads

1、安装，构建数据库