通过GeneId在NCBI上批量搜寻序列

来源：互联网发布：香橙派和树莓派知乎编辑：程序博客网时间：2024/04/28 19:16

当你有大批量的GeneId序列的时候，手动一个个在NCBI里面比对肯定是不现实的，本来有Batch Entrez方便批量比对，可是总容易不允许进入。所以写了以下一种流程可以批量比对。

首先如果是比对蛋白质（protein）或者(nucleotide)数据库，那很容易用perl以下代码即可。

#!/usr/bin/perl -wuse Bio::SeqIO;use Bio::DB::GenBank;my $myIdFile="example.txt";#GeneId所在文件;my $outputDir="result/";#输出目录open(myFile,$myIdFile)||die();while(my $myLine=<myFile>){chomp($myLine);$out = Bio::SeqIO->new(-file => ">".$outputDir.$myLine."_protein.fasta" , '-format' => 'fasta');my $query_string = $myLine;my $query = Bio::DB::Query::GenBank->new(          -query=>$query_string,         -db      => 'protein',);# get a genbank database handle$gb = new Bio::DB::GenBank; eval{my $stream = $gb->get_Stream_by_query($query);while (my $seq = $stream->next_seq) {$out->write_seq($seq);}};}

但是如果要是想比对GENE数据库时候，就麻烦了。因为不提供直接远程比对GENE数据库服务库的。Bio::DB::GenBank;模块有对此说明的。

通过对下载链接http://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?tool=portal&sendto=on&log$=seqview&db=nuccore&dopt=fasta&sort=&val=566192888&from=19085984&to=19088706&maxplex=1的分析我们知道需要一个giID，和起始结束的位置。

那我们的思路就是

1.首先进入到http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=gene&id=POPTR_0010s21980g&retmode=xml中间找到giID<Seq-id_gi>568815590</Seq-id_gi>

2.然后再进入到NCBI中GENE库中搜索geneID找寻到起始结束位置比如http://www.ncbi.nlm.nih.gov/gene/?term=POPTR_0010s21980g。

3.因为这两个页面所需求内容都在网页源代码里面，所以都可以直接用LWP::Simple 模块得到相应内容。然后依据 giID 和起始结束位置就能批量产生下载地址。

4.这也是一个问题就是下载地址必须在地址栏打开后通过迅雷嗅针自动打开迅雷下载才行，如果直接复制到迅雷批量下载中是不可以的。所以此处可以写一个javascript脚本+html页面模拟实现即可。（注意在浏览器设置中优先选择迅雷下载，否则不能批量下载）。

好，重点来了。贴上代码：

1.通过genID集合文件找寻giID结果文件输出到seqId.txt文件中。

#!/usr/bin/perl -wuse strict;use LWP::Simple qw(get);my $myIdFile="ee.txt";#geneId号所在文件。my $seq_id=0;my $html="";my @array;open(myFile,$myIdFile)||die();while(my $myUrl=<myFile>){chomp($myUrl);$html = get( "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=gene&id=$myUrl&retmode=xml" );@array=split/\n/,$html;foreach my $myLine(@array){if($myLine=~/<Seq-id_gi>[\d0-9]+<\/Seq-id_gi>/){if ($myLine =~ /<Seq-id_gi>(\S+)<\/Seq-id_gi>/){ $seq_id=$1; }}}open(OutFile,">>seqId.txt")||die();print OutFile $myUrl."&".$seq_id."\n";close(OutFile);}close(myFile);

2.通过genID集合文件找寻起始结束位置。

#!/usr/bin/perl -wuse strict;use LWP::Simple qw(get);my $myIdFile="ee.txt";#GeneId号所在位置。my $begin=0;my $end=0;my $html="";my @array;open(myFile,$myIdFile)||die();while(my $myUrl=<myFile>){chomp($myUrl);$html = get( "http://www.ncbi.nlm.nih.gov/gene/?term=$myUrl" );@array=split/\n/,$html;foreach my $myLine(@array){if($myLine=~/\([\d0-9]+\.\.[\d0-9]/){if ($myLine =~ /\((\S+)\.\.(\S[\d0-9]+)/){ $begin=$1; $end=$2; }}}open(OutFile,">>postion.txt")||die();print OutFile $myUrl."&".$begin."&".$end."\n";close(OutFile);}close(myFile);

3.将两个结果生成url地址。

#/usr/bin/perl -wuse strict;my $begin=0;my $end=0;my $seq_id=0;my $seqIdFile="seqId.txt";#giID号my $posFile="postion.txt";#gene所在位置my $keyWord="";my $mySourceUrl="";open(Myfile,$seqIdFile)||die();while(my $Line=<Myfile>){    if($Line=~/(\S+)\&(\S+)/){        $keyWord=$1;        $seq_id=$2;    }   open(myPosFile,$posFile)||die();   while(my $posLine=<myPosFile>){      if($posLine=~/(\S+)\&(\S+)\&(\S+)/){                   my $midnum=$2;                   my $midnum2=$3;           if($keyWord eq $1){                $begin=$midnum;                $end=$midnum2;           }      }   }   close(myPosFile);   $mySourceUrl="http://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?tool=portal&sendto=on&log\$=seqview&db=nuccore&dopt=fasta&sort=&val=$seq_id&from=$begin&to=$end&maxplex=1";              open(OUTFILE,">>url.txt")||die();              print OUTFILE $mySourceUrl."\n";              close(OUTFILE);}close(Myfile);

4.打开url.html将上面生成的所有地址全部放入文本框中（PS：最好用crtl+A选中文本框所有内容后再crtl+V粘贴进去，防止多一行空白行）

<!DOCTYPE html><html xmlns="http://www.w3.org/1999/xhtml"><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>    <title></title></head><body><div >    <textarea id="myText">    </textarea></div><div >    <input type="button" onclick="geturl()" value="提交"/></div></body></html><script type="text/ecmascript">    function geturl() {        var myUrl = document.getElementById("myText").value;        var arr = myUrl.split('\n');        var i=0;        var mytime1=setInterval(openUrl,5000);        function openUrl() {                       window.open(arr[i],'_self');                        if (i == arr.length) {                             clearInterval(mytime1);            }            else {                i++;            }        }    }</script>

5.等待就好，迅雷会自动合并提示是否下载，等到下载地址全部加载完毕，合并任务组下载（PS ：一定点击方便后面处理）如果仍然卡住可以先在NCBI中任意下一条fasta序列，然后再进行第四步骤。不过记得要清空迅雷下载内容，否则出现重复下载也不行。

6.下载完毕，这时候序列都下载好了，但是名字不对对不对？而且fasta文件中>后面内容也不包含GENEID号。这时候需要批量对文件名字处理，要用到之前生成的giID号和起始结束位置文件。然后先用下面的代码把两个文件集合起来生产jihe.txt文件。

#/usr/bin/perl -wuse strict;my $begin=0;my $end=0;my $seq_id=0;my $seqIdFile="seqId.txt";my $posFile="postion.txt";my $keyWord="";my $mySourceUrl="";open(Myfile,$seqIdFile)||die();while(my $Line=<Myfile>){    if($Line=~/(\S+)\&(\S+)/){        $keyWord=$1;        $seq_id=$2;    }   open(myPosFile,$posFile)||die();   while(my $posLine=<myPosFile>){      if($posLine=~/(\S+)\&(\S+)\&(\S+)/){           if($keyWord eq $1){                $begin=$2;                $end=$3;           }      }   }   close(myPosFile);              open(OUTFILE,">>jihe.txt")||die();              print OUTFILE $keyWord."&".$seq_id."&".$begin."&".$end."\n";              close(OUTFILE);}close(Myfile);

7.再把jihe.txt文件和刚才任务组下载的文件夹平级，再调用下面的perl程序。生成正确名字外加正确的序列文件。然后将以前的文件手动删除下即可。注意：find.txt里面表示正确找到底基因序列。如果发现数目少了，可以进行比对。剩下的可以手动比对去找。

#/usr/bin/perl -wuse strict;my $mydir="任务组_20150404_1447/";#下载的任务组文件夹my $mySeqId="jihe.txt";#整合后的文件my $outName="";my $midName="";my $result="";my $myKeyWord="";opendir(DIR,$mydir)||die();while(my $myFile=readdir(DIR)){  if($myFile=~/[\da-z]/){  $result="";  open(myDirFile,$mydir.$myFile)||die();      $outName="";  while(my $myLine=<myDirFile>){  if($myLine=~/>/){         open(mySeqId,$mySeqId)||die();       while(my $mySeqLine=<mySeqId>){          if($mySeqLine=~/(\S+)\&(\S+)\&(\S+)\&(\S+)/){                   $midName=$1;                 $myKeyWord=$2.":".$3."-".$4;              if($myLine=~/$myKeyWord/){              $outName=$midName;              }          }  }  close(mySeqId);  if($myLine=~/(\S+) ?/){my @array=split/$1/,$myLine;     $result.=">".$outName." ".$array[1];}    }  $result.=$myLine;  }  close(myDirFile);  open(OutFile,">".$outName.".fasta")||die();  print OutFile $result;  close(OutFile); open(OutFile,">>find.txt")||die();  print OutFile $outName."\n";  close(OutFile);  }}close(DIR);

结束，总共花了三天时间摸索出来的一套流程。对于没有服务器做本地比对的童鞋来说，还是可以方便借鉴的。外加一句，摸索之路真是吐血。

0 0