哈希表解决提取reads问题

哈希表解决提取reads问题,第1张

概述问题描述,有两个文件,一个是基因和其相对应所包含的readsID以及是正反向测序是否都要的标识。 格式如下:Tea_CCG055929.1&FCC55J8ACXX:8:1101:6302:2021#TAGCTTAT/1&A 基因名&reads名&标识。如果A则正反都要,如果S则只要一条。 另一个则是reads文件标准fq格式。 依据关键字文件找出对应的reads序列存入相应基因名字建立的文件中。 问题描述,有两个文件,一个是基因和其相对应所包含的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 -wuse 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=;          if($Readsname=~/(\S+)\:(\S+):(\S+):(\S+):(\S+)\#/){         $mybeginAdd=*$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=*$myAmplify;         if($mybeginAdd>$MAXSIZE){             next;         }            while(!($array[$mybeginAdd] eq "END")){            if($array[$mybeginAdd]=~/(\S+)\&(\S+)\&(\S+)/){              $Genename=;              $Readsname=;              $myFlag=;              if($myFlag eq "A"){                 if($Readsname=~/(\S+)\/(\S+)/){                    $Readsname=;                 }              }                               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);}
总结

以上是内存溢出为你收集整理的哈希表解决提取reads问题全部内容,希望文章能够帮你解决哈希表解决提取reads问题所遇到的程序开发问题。

如果觉得内存溢出网站内容还不错,欢迎将内存溢出网站推荐给程序员好友。

欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/langs/1273825.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2022-06-08
下一篇 2022-06-08

发表评论

登录后才能评论

评论列表(0条)

保存