哈希表解决提取reads问题

最新推荐文章于 2024-11-05 01:41:28 发布

澜澜学编程

最新推荐文章于 2024-11-05 01:41:28 发布

阅读量579

点赞数

分类专栏：算法生物信息文章标签： perl 算法

本文链接：https://blog.csdn.net/wjlhxm123/article/details/45111059

版权

算法生物信息专栏收录该内容

2 篇文章 0 订阅

订阅专栏

问题描述，有两个文件，一个是基因和其相对应所包含的readsID以及是正反向测序是否都要的标识。

格式如下：Tea_CCG055929.1&FCC55J8ACXX:8:1101:6302:2021#TAGCTTAT/1&A

基因名&reads名&标识。如果A则正反都要，如果S则只要一条。

另一个则是reads文件标准fq格式。

依据关键字文件找出对应的reads序列存入相应基因名字建立的文件中。

reads文件大约13G，关键字文件大约2G。如果用平常循环比对，耗时将非常长。时间复杂度也为O（mn），若采用哈希表存储结构先将关键字文件存储进去，将reads文件一一比对，时间复杂度将降低到O(n)。

哈希表存储结构关键是哈希函数的构造，用perl写了一个统计关键字文件有多少行以及最后一个数字最大为多少（如例子格式最后一个数字为2021）；最终统计出总共有35741312条数据，最大值为100685 。装载因子采用0.7 所以哈希函数为 hash（KEY）=KEY×500; 存储数组为0-50342500。冲突解决方式为开放地址法，最终地址为H=hash（KEY）+di；

代码如下：

#/usr/bin/perl -w
use strict;
my $myReads="read1.fq";#reads文件
my $myFileAdd=".fastq1";
my $outDir="gene/";#输出目录
my $mySmalt="example.txt";#smalt结果文件

my $MAXSIZE=50342500;
my $myAmplify=500;
my @array;
my $i=0;
my $mybeginAdd=0;
my $GeneName="";
my $ReadsName="";
my $myFlag="";
my $outFlag=0;
my $countNum=0;
my $myResult="";
for($i=0;$i<=$MAXSIZE;$i++){
   $array[$i]="END";
}
open(myFile,$mySmalt)||die("connot open file 1");
while(my $myLine=<myFile>){
  if($myLine=~/(\S+)\&(\S+)\&(\S+)/){
          $ReadsName=$2;
          if($ReadsName=~/(\S+)\:(\S+):(\S+):(\S+):(\S+)\#/){
         $mybeginAdd=$5*$myAmplify;

            $mybeginAdd=&myHash($mybeginAdd);

             $array[$mybeginAdd]=$myLine;

           }
           }

}
close(myFile);

open(myFile,$myReads)||die("connot open file 2");
while(my $myLine=<myFile>){
chomp($myLine);
if($myLine=~/@/){
     if($myLine=~/\@(\S+)\:(\S+):(\S+):(\S+):(\S+)\#/){
         $mybeginAdd=$5*$myAmplify;

         if($mybeginAdd>$MAXSIZE){
             next;
         }
            while(!($array[$mybeginAdd] eq "END")){

            if($array[$mybeginAdd]=~/(\S+)\&(\S+)\&(\S+)/){
              $GeneName=$1;
              $ReadsName=$2;
              $myFlag=$3;
              if($myFlag eq "A"){
                 if($ReadsName=~/(\S+)\/(\S+)/){
                    $ReadsName=$1;

                 }

              }
                               if($myLine=~/$ReadsName/){

                               $outFlag=1;
                               $countNum=0;
                               open(OUTFILE,">>".$outDir.$GeneName.$myFileAdd)||die("cannot open outfile");
                               print OUTFILE $myLine."\n";
                               close(OUTFILE);
                               last;
                               }

              }

              $mybeginAdd++;
              if($mybeginAdd>$MAXSIZE){
                         $mybeginAdd=$mybeginAdd%$MAXSIZE;
                      }
            }

           }
     }
     elsif($outFlag==1){
        if($countNum==2){
          $myResult.=$myLine."\n";
          open(OUTFILE,">>".$outDir.$GeneName.$myFileAdd)||die("cannot open outfile");
                               print OUTFILE $myResult;
                               close(OUTFILE);
                               $outFlag=0;
                               $myResult="";
        }
        else
        {
          $myResult.=$myLine."\n";
        }
        $countNum++;

     }
}
close(myFile);

sub myHash{
                 my ($myAdd) =@_;

                 $myAdd=$myAdd+0;

                 while(!($array[$myAdd] eq "END")){
                      $myAdd++;
                      if($myAdd>$MAXSIZE){

                         $myAdd=$myAdd%$MAXSIZE;
                      }

                 }
                   return ($myAdd);
}