K-mer Index问题
|
问题：给定一个DNA序列，这个系列只含有4个字母ATCG，如 S =“CTGTACTGTAT”。给定一个整数值k，从S的第一个位置开始，取一连续k个字母的短串，称之为k-mer（如= 5，则此短串为CTGTA），然后从的第二个位置，取另一（如TGTAC），这样直至的末端，就得一个集合，包含全部。 如对序列S来说，所有5-mer为｛GTACTTACTGACTGTTGTAT｝通常这些k-mer需一种数据索引方法，可被后面的操作快速访问。例如，对5-mer来说，当查询CTGTA，通过这种数据索引方法，可返回其在DNA序列S中的位置为｛1，6｝。 现在以文件形式给定 多个 DNA序列，要求对给定k， 给出并实现一种数据索引方法，可返回任意一个k-mer所在的DNA序列编号和相应序列中出现的位置。 思路： 1、首先采用命A=0,C=1,G=2,T=3. 就相当于4进制数字，然后采用karp-Rabin算法转换成唯一十进制数字。由于用此算法的哈希函数为：hash(value)=value*(4^(k-q-1)); value是该字符对应的值，k是kmer长度,q是此字符在字符串的位置范围在[0-(q-1)]。然后把一个kmer里面所有字符的hash值求和就行了。 2、那么很容易看出来，对于连续的下一个Kmer，就有推理公式了?hashNew=addValue+(hashOld-deleteValue*(4^(k-1)))*4; hashNew就是往右平移一个字符的kmer hash值，hashOld就是平移之前的值，addValue就是平移后右边多的一个字符，deleteValue就是平移后左边少的一个字符。这样整个hash表建立的时间复杂度约为O（m+k）,m是整个文本长度。 3、由于kmer长度如果过长，其hash值过大，会造成内存不够溢出的现象，所以kmer内部定死为10 。那么问题就来了，如何应对不同的kmer值。分三种情况。 第一种：q>10? 这种可以将kmer以10为单位，将hash表中对应值取出，然后对结果进行分析，这边分析方法为建立两个数组一个二维数组unionName储存位置关系，一个一维数组unionScore,计数用。 思路就是首先第一轮初始化unionName[Name][Pos]全部赋值Pos?并初始化unionScore，然后再第二轮匹配如果[Name][Pos-cycle]=Pos-1则将其赋值为当前Pos，cycle为当前循环次数。并将当前循环数存入unionScore[NAME]中。最后当unionScore[NAME]值也就是循环数为k-1，即我们需要的交集了。 第二种：q=10 直接求出hash值，取出相应的值即可。 第三种：q<10 可以用前缀种子+后缀种子交集产生。 前缀种子：在字符串后面补字符直到长度等于K，这个很容易看出来 最小是全补A，最大是全补T，然后将最小值到最大值之间的hash值即为所求。 后缀种子：后缀种子和前缀种子不同就是在字符串左边补齐字符。所以此时需要进行变换。只要对前置种子产生的值变化下就行了。(preValue-minValue)*(4^(K-q))+hash(p) 。其中preValue就是对应的前置种子的hash值，minValue就是前置种子中最小值也就是全补A的情况，hash（p）就是字符串长度为p时候的hash值。? 交集就是先求后缀种子所有的值，再加上 前缀种子中起始位置在[0-(k-1)]中的值。 代码如下：
#usr/bin/perl -w
use strict;
my $myFileName="test.txt";#需要找寻的文件
my $searchKey="GGGA";#找寻的字符
my $myKmerInput=length $searchKey;#K-mer值
my $myKmer=2;
my @name={};
my $tempLength=0;
my $i=0;
my $j=0;
my @Add={};
my @temp={};
my $newHashValue=0;
my $myHashOld=0;
my $flag=0;
my $nameCycle=0;
my $finaNameFlag=0;
if($myKmerInput<=$myKmer){
$flag=1;
}
open(myFile,$myFileName)||die("cannot open $myFileName");
while(my $myLine=<myFile>){
chomp($myLine);
if($myLine=~/>/){
if($myLine=~/>(S+) /){
$name[$nameCycle]=$1;
}
}
elsif($myLine=~/[A-Z]/){
@temp=split//,$myLine;
$tempLength=scalar @temp;
if($tempLength-$myKmer+1 <= 0){
die("Kemer Length is too long");
}
$myHashOld=&myHash($myKmer,@temp);
$Add[$myHashOld].="@1&$nameCycle";
for($i=1;$i<$tempLength-$myKmer+1;$i++){
$j=$i+1;
$newHashValue=&rehash($myHashOld,$temp[$i-1],$temp[$i+$myKmer-1]);
$Add[$newHashValue].="@$j&$nameCycle";
$myHashOld=$newHashValue;
}
$nameCycle++;
}
else{
next;
}
}
close(myFile);
@temp=split//,$searchKey;
if($flag==0){
my @addTemp={};
$tempLength=scalar @temp;
if($tempLength-$myKmer+1 <= 0){
die("Kemer Length is too long");
}
$myHashOld=&myHash($myKmer,@temp);
$addTemp[0]=$Add[$myHashOld];
for($i=1;$i<$tempLength-$myKmer+1;$i++){
$j=$i+1;
$newHashValue=&rehash($myHashOld,$temp[$i+$myKmer-1]);
if($Add[$newHashValue]=~/[0-9]/){
$addTemp[$i]=$Add[$newHashValue];
}
$myHashOld=$newHashValue;
}
if((scalar @addTemp)!=($tempLength-$myKmer+1)){
die("cannot findn");
}
@addTemp=&unionValue(@addTemp);
for($i=0;$i<$finaNameFlag;$i++)
{
if($addTemp[$i]=~/(S+)&;(S+)/){
my $posEnd=$2+$myKmerInput-1;
print "myDNAName=$name[$1]tposBegin=$2tposEnd=$posEndn";
}
}
}
elsif($flag==1){
my $resultFalg=0;
my $minValue=&myHash($myKmerInput,@temp);
my $nextMinValue=$minValue;
if(($myKmer-$myKmerInput)>=1){
for($i=0;$i<$myKmer-$myKmerInput;$i++){
$temp[($i+$myKmerInput)]="A";
}
$minValue=&myHash($myKmer,@temp);
for($i=0;$i<$myKmer-$myKmerInput;$i++){
$temp[($i+$myKmerInput)]="T";
}
}
my $maxValue=&myHash($myKmer,@temp);
my $myNowValue=$minValue;
for($myNowValue;$myNowValue<=$maxValue;$myNowValue++){
my @addTemp={};
@addTemp=split/@/,$Add[$myNowValue];
my $addTempLength=scalar @addTemp;
if($addTempLength==0){
next;
}
for($i=0;$i<$addTempLength;$i++){
if($addTemp[$i]=~/(S+)&;(S+)/){
if($1>$myKmer&&($myKmer!=$myKmerInput)){
next;
}
else{
my $posEnd=$1+$myKmerInput-1;
print "myDNAName=$name[$2]tposBegin=$1tposEnd=$posEndn";
$resultFalg=1;
}
}
}
}
if(($myKmer-$myKmerInput)>=1){
$myNowValue=$minValue;
my $nextMaxValue=($maxValue-$minValue)*(4**($myKmer-$myKmerInput))+$nextMinValue;
for($myNowValue;$myNowValue<=$maxValue;$myNowValue++){
my @addTemp={};
my $nextNowValue=($myNowValue-$minValue)*(4**$myKmerInput)+$nextMinValue;
@addTemp=split/@/,$Add[$nextNowValue];
my $addTempLength=scalar @addTemp;
if($addTempLength==0){
next;
}
for($i=0;$i<$addTempLength;$i++){
if($addTemp[$i]=~/(S+)&;(S+)/){
if($1==1){
next;
}
else{
my $nextPosBegin=$1+$myKmer-$myKmerInput;
my $nextPosEnd=$nextPosBegin+$myKmerInput-1;
print "myDNAName=$name[$2]tposBegin=$nextPosBegintposEnd=$nextPosEndn";
$resultFalg=1;
}
}
}
}
}
if(!$resultFalg){
die("cannot findn");
}
}
# for($i=0;$i<scalar @Add;$i++){
# if(!($Add[$i] eq "")){
# print $i."t".$Add[$i]."n";
# }
# }
sub unionValue{
my(@unionTemp)=@_;
my $unionTempLength=scalar @unionTemp;
my $unionCycle=0;
my $unionNameBefore=0;
my $unionNameNext=0;
my $unionPosBefore=0;
my $unionPosNext=0;
my @unionTempBefore={};
my @unionTempName=();
my @unionScoreName={};
my @finalName={};
my $myDNAName="";
my $unionTempCycle=0;
@unionTempBefore=split/@/,$unionTemp[0];
my $unionTempBeforeLength=scalar @unionTempBefore;
for($unionTempCycle=1;$unionTempCycle<$unionTempBeforeLength;$unionTempCycle++){
if($unionTempBefore[$unionTempCycle]=~/(S+)&;(S+)/){
$unionNameBefore=$2;
$unionPosBefore=$1;
$unionTempName[$unionNameBefore][$unionPosBefore]= $unionPosBefore;
$unionScoreName[$unionNameBefore]=0;
}
}
for($unionCycle=1;$unionCycle<$unionTempLength;$unionCycle++){
@unionTempBefore=split/@/,$unionTemp[$unionCycle];
$unionTempBeforeLength=scalar @unionTempBefore;
for($unionTempCycle=1;$unionTempCycle<$unionTempBeforeLength;$unionTempCycle++){
if($unionTempBefore[$unionTempCycle]=~/(S+)&;(S+)/){
$unionNameBefore=$2;
$unionPosBefore=$1;
}
if($unionPosBefore<$unionCycle){
next;
}
if($unionTempName[$unionNameBefore][$unionPosBefore-$unionCycle]==($unionPosBefore-1)&&(($unionPosBefore-1)!=0)){
$unionTempName[$unionNameBefore][$unionPosBefore-$unionCycle]=$unionPosBefore;
$unionScoreName[$unionNameBefore]=$unionCycle;
if($unionScoreName[$unionNameBefore]==$unionTempLength-1){
$finalName[$finaNameFlag++]=$unionNameBefore."&".($unionPosBefore-$unionCycle);
}
}
}
}
if($finaNameFlag==0){
die("cannot findn");
}
return @finalName;
}
sub rehash{
my ($myOldValue,$deleteChar,$addChar) = @_;
my $myNewValue=0;
$myNewValue=&testChar($addChar)+($myOldValue-&testChar($deleteChar)*(4**($myKmer-1)))*4;
return $myNewValue;
}
sub myHash{
my ($myLength,@list) = @_;
my $result=0;
for(my $k=0;$k<$myLength;$k++){
$result+= &testChar($list[$k])*(4**($myLength-$k-1));
}
return $result;
}
sub testChar{
my ($myChar) = @_;
if($myChar eq 'A')
{
return 0;
}
elsif($myChar eq 'C'){
return 1;
}
elsif($myChar eq 'G'){
return 2;
}
elsif($myChar eq 'T'){
return 3;
}
else{
die("wrong DNA");
}
}
（编辑：李大同） 【声明】本站内容均来自网络，其相关言论仅代表作者个人观点，不代表本站立场。若无意侵犯到您的权利，请及时与联系站长删除相关内容! |