哈希表解决提取reads问题
发布时间:2020-12-15 23:45:58 所属栏目:大数据 来源:网络整理
导读:问题描述,有两个文件,一个是基因和其相对应所包含的readsID以及是正反向测序是否都要的标识。 格式如下:Tea_CCG055929.1FCC55J8ACXX:8:1101:6302:2021#TAGCTTAT/1A 基因名reads名标识。如果A则正反都要,如果S则只要一条。 另一个则是reads文件标准fq格式
问题描述,有两个文件,一个是基因和其相对应所包含的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); } (编辑:李大同) 【声明】本站内容均来自网络,其相关言论仅代表作者个人观点,不代表本站立场。若无意侵犯到您的权利,请及时与联系站长删除相关内容! |