群体遗传研究——变异注释--snpEff软件进行建库
2022-09-24 本文已影响0人
千万别加香菜
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信息即可，则最后出现的文件小一些
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)
