通過GeneId在NCBI上批量搜尋序列

當你有大批量的GeneId序列的時候，手動一個個在NCBI裡面比對肯定是不現實的，本來有Batch Entrez友善批量比對，可是總容易不允許進入。是以寫了以下一種流程可以批量比對。

首先如果是比對蛋白質（protein）或者(nucleotide)資料庫，那很容易用perl以下代碼即可。

#!/usr/bin/perl -w
use 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 -w
use 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 -w
use 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 -w
use 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" οnclick="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 -w
use 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 -w
use 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);

結束，總共花了三天時間摸索出來的一套流程。對于沒有伺服器做本地比對的童鞋來說，還是可以友善借鑒的。外加一句，摸索之路真是吐血。

通過GeneId在NCBI上批量搜尋序列

繼續閱讀

AWStats安裝與使用

安裝AWStats

日志分析軟體awstats的安裝使用指南

AWStats使用

linux下安裝awstats6.5

轉載唐鳳的部落格 :: Perl 6 十周年慶

WINDOWS下安裝MRTG全攻略

configure/make/make install的作用

ubuntu下gvim配置檔案.vimrc

SPOJ QTREE4 Query on a tree IV

如何配置Eclipse進行Perl開發

npm install stylus --save失敗

Excel多行轉多列

在Windows上編譯Wireshark源代碼 .

Learning Perl: 1.3. How Can I Get Perl?

Perl與網絡監控