以NBS-LRR基因家族为例
1、下载数据
从核苷酸数据库下载所有的NBS-LRR基因的概要文件
如果fasta数据下载，太慢了。如果网络好可以全部下载下来，然后本地写脚本
2、从概要中提取基因ac号，并分类
#! /usr/bin/perl
#处理从ncbi下载下来的summary文件
use autodie;
$file=`find ncbi*.txt`;#打开ncbi下载的文件,如果有多个文件建议手动输入
open IN1 ,"$ARGV[0]"|| $file;
#第一步，统计每个物种的NBS-LRR基因数量
map{$species{$1}++ if(/\[(.+?)\]/)}<IN1>;
#打印结果 种类=>基因个数
map{printf "%-45s => %-s \n",$_,$species{$_};}sort keys %species;
#第二步，提取各物种的nbs-lrr的蛋白质的ac号
seek(IN1,0,0);#重置到文件开头
while(<IN1>){
$species=$1 if(/\[(.+?)\]/);
$ac{$species} .= "$1\t" if(/(.+) GI:/);
}
#输出结果到屏幕
for(sort{$species{$a}<=>$species{$b}} keys %ac){
print "$_:\n$ac{$_}\n";
}
# summary 文件数据样式
#1. NBS-LRR, partial [Musa AAB Group]
#172 aa protein
#AEQ16455.1 GI:349734021
#
#2. NBS-LRR [Panicum miliaceum]
#712 aa protein
#RLN11370.1 GI:1492243987
image.png
3、利用分类的基因ac号，找到感兴趣的物种，将其ac号复制到另外一个文件中去。整个为一行中间用\t隔开。
图片.png
4、用这个文件上传到ncbi下载序列
网址
https://www.ncbi.nlm.nih.gov/sites/batchentrez
图片.png
图片.png
5、如果你下载了全部序列，则可以自己写个脚本，本地提取更方便。
#! /usr/bin/perl -w
use autodie;
open IN,"$ARGV[0]";#ac号文件
open IN2,"$ARGV[1]";#下载的全部序列文件
open OUT,"> ac_seq.fasta";
map{$ac .= "$_"; }<IN>;
map{
if(/>(\S+)/){
$flag=0;
$count+=$flag=1 if($ac=~/$1/)
}
print OUT "$_"if(1==$flag)
}<IN2>;
print"$count";
#ac号文件样式 各个号码之间以\t隔开
#CAB38788.1 CAB38849.1 CAB79574.1 CAB80047.1 AAG12572.1
#AAC64218.1 AAC26126.1AAF19533.1 AAG52149.1 AAG52150.1
图片.png
