格式如下: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问题所遇到的程序开发问题。
如果觉得内存溢出网站内容还不错,欢迎将内存溢出网站推荐给程序员好友。
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)