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信息即可，则最后出现的文件小一些
脚本TratesnpEffOutfile.py
用于对snpEff软件注释后的.anno.bed文件进行提取，提取后的每个变异占一行，方便查看对应的选择信号值与变异的关系
#!/usr/bin/env python
# -*- encoding: utf-8 -*-
'''
@File : TratesnpEffOutfile.py
@Time : 2024/05/29 09:20:14
@Author : Lulu Shi
@Mails : crazzy_rabbit@163.com
@link : https://github.com/Crazzy-Rabbit
'''
import click
@click.command()
@click.option('--infile', type=click.File('r'), help="the out .anno.bed file that annotated by snpEff")
@click.option('--outfile', type=click.File('w'), help='out file name for you extract')
def main(infile, outfile):
"""
用于对snpEff软件注释后的.anno.bed文件进行提取，\n
提取后的每个变异占一行，方便查看对应的选择信号值与变异的关系
"""
for line in infile:
if line.startswith('#'):
new_line = '\t'.join(line.split('|'))
outfile.write('\t'.join(new_line.split(';')))
else:
parts = line.strip().split(';')
if len(parts) >= 2:
new_parts1 = '\t'.join(parts[1].split('|'))
outfile.write(parts[0] + '\t' + new_parts1 + "\n")
for content in parts[2:]:
new_content = '\t'.join(content.split('|'))
outfile.write(parts[0] + '\t' + new_content +'\n')
print(f"处理完成，结果已保存")
if __name__ == '__main__':
main()
处理后效果如下
# SnpEff version 5.0 (build 2020-08-09 21:23), by Pablo Cingolani
# Command line: SnpEff sheepv1.0 -i bed all.chart.ihs.out.50k.windows.5%.sort.chr.bed
# Chromo Start End Name Effect Gene BioType Score
NC_040258.1 80900001 80950001 0.175676 Upstream:29974 Transcript:XM_015097123.2:protein_coding Gene:FUT8:protein_coding
NC_040258.1 80900001 80950001 0.175676 Upstream:28742 Transcript:XM_012181844.3:protein_coding Gene:FUT8:protein_coding
NC_040258.1 80900001 80950001 0.175676 Exon:1:12:RETAINED Transcript:rna-XM_015097123.2 Gene:exon-XM_015097123.2-1:null
NC_040258.1 80900001 80950001 0.175676 Upstream:29974 Transcript:rna-XM_015097123.2 Gene:exon-XM_015097123.2-1:null
NC_040258.1 80900001 80950001 0.175676 Transcript:XM_012181844.3:protein_coding Gene:FUT8:protein_coding
NC_040258.1 80900001 80950001 0.175676 Transcript:XM_015097123.2:protein_coding Gene:FUT8:protein_coding
NC_040258.1 80900001 80950001 0.175676 Exon:1:12:RETAINED Transcript:rna-XM_012181844.3 Gene:exon-XM_012181844.3-1:null
NC_040258.1 80900001 80950001 0.175676 Upstream:28742 Transcript:rna-XM_012181844.3 Gene:exon-XM_012181844.3-1:null
NC_040258.1 80900001 80950001 0.175676 Intergenic:LOC101111101-GENE_exon-XM_012181844.3-1
NC_040260.1 77300001 77350001 0.175676 Intergenic:GENE_exon-XM_015097839.2-1-LOC114116467
NC_040277.1 39800001 39850001 0.175772 Exon:5:6:RETAINED Transcript:rna-XM_012106001.3 Gene:exon-XM_012106001.3-1:null
NC_040277.1 39800001 39850001 0.175772 Transcript:XM_004021795.4:protein_coding Gene:GOLGA7:protein_coding
NC_040277.1 39800001 39850001 0.175772 Exon:6:8:RETAINED Transcript:rna-XM_015104630.2 Gene:exon-XM_015104630.2-1:null
NC_040277.1 39800001 39850001 0.175772 Exon:3:8:RETAINED Transcript:rna-XM_015104630.2 Gene:exon-XM_015104630.2-1:null
NC_040277.1 39800001 39850001 0.175772 Exon:6:6:RETAINED Transcript:rna-XM_004021795.4 Gene:exon-XM_004021795.4-1:null
NC_040277.1 39800001 39850001 0.175772 Transcript:XM_015104630.2:protein_coding Gene:GINS4:protein_coding
NC_040277.1 39800001 39850001 0.175772 Downstream:0 Transcript:rna-XM_015104630.2 Gene:exon-XM_015104630.2-1:null
NC_040277.1 39800001 39850001 0.175772 Exon:7:8:RETAINED Transcript:rna-XM_015104630.2 Gene:exon-XM_015104630.2-1:null
NC_040277.1 39800001 39850001 0.175772 Upstream:27365 Transcript:rna-XM_015104630.2 Gene:exon-XM_015104630.2-1:null
NC_040277.1 39800001 39850001 0.175772 Downstream:0 Transcript:XM_015104630.2:protein_coding Gene:GINS4:protein_coding
NC_040277.1 39800001 39850001 0.175772 Exon:2:8:RETAINED Transcript:rna-XM_015104630.2 Gene:exon-XM_015104630.2-1:null
NC_040277.1 39800001 39850001 0.175772 Downstream:0 Transcript:rna-XM_012106001.3 Gene:exon-XM_012106001.3-1:null
NC_040277.1 39800001 39850001 0.175772 Transcript:XM_012106001.3:protein_coding Gene:GOLGA7:protein_coding
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=<RF>)
{
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=<RF>)
{
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)
