加入收藏 | 设为首页 | 会员中心 | 我要投稿 李大同 （https://www.lidatong.com.cn/）- 科技、建站、经验、云计算、5G、大数据,站长网!
当前位置： 首页 > 大数据 > 正文

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

发布时间：2020-12-15 23:46:07 所属栏目：大数据 来源：网络整理
导读：当你有大批量的GeneId序列的时候，手动一个个在NCBI里面比对肯定是不现实的，本来有Batch Entrez方便批量比对，可是总容易不允许进入。所以写了以下一种流程可以批量比对。 首先如果是比对蛋白质（protein）或者(nucleotide)数据库 ，那很容易 用perl以下代

当你有大批量的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" 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 -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,$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);
结束，总共花了三天时间摸索出来的一套流程。对于没有服务器做本地比对的童鞋来说，还是可以方便借鉴的。外加一句，摸索之路真是吐血。

（编辑：李大同）

【声明】本站内容均来自网络，其相关言论仅代表作者个人观点，不代表本站立场。若无意侵犯到您的权利，请及时与联系站长删除相关内容!
    相关内容
      推荐文章
        站长推荐
        热点阅读

          【免责声明】本站内容转载自互联网，其发布内容言论不代表本站观点，如果其链接、内容的侵犯您的权益，烦请提交相关链接至邮箱bqsm@foxmail.com我们将及时予以处理。

          建议您使用1920×1080分辨率、谷歌浏览器Google Chrome、Microsoft Edge以获得本站的最佳浏览效果

          Copygight © 2008-2022 https://www.lidatong.com.cn/ All Rights Reserved. 李大同