phmm.pl

一个用于蛋白质序列分类的profile hmm的Perl代码

print "\$Mstate[$i]{\"P\"}=".$Mstate[$i]{"P"}."\n";
print "\$Mstate[$i]{\"Q\"}=".$Mstate[$i]{"Q"}."\n";
print "\$Mstate[$i]{\"R\"}=".$Mstate[$i]{"R"}."\n";
print "\$Mstate[$i]{\"S\"}=".$Mstate[$i]{"S"}."\n";
print "\$Mstate[$i]{\"T\"}=".$Mstate[$i]{"T"}."\n";
print "\$Mstate[$i]{\"V\"}=".$Mstate[$i]{"V"}."\n";
print "\$Mstate[$i]{\"W\"}=".$Mstate[$i]{"W"}."\n";
print "\$Mstate[$i]{\"X\"}=".$Mstate[$i]{"X"}."\n";
print "\$Mstate[$i]{\"Y\"}=".$Mstate[$i]{"Y"}."\n";
print "\n";
print "\$Istate[$i]{\"A\"}=".$Istate[$i]{"A"}."\n";
print "\$Istate[$i]{\"C\"}=".$Istate[$i]{"C"}."\n";
print "\$Istate[$i]{\"D\"}=".$Istate[$i]{"D"}."\n";
print "\$Istate[$i]{\"E\"}=".$Istate[$i]{"E"}."\n";
print "\$Istate[$i]{\"F\"}=".$Istate[$i]{"F"}."\n";
print "\$Istate[$i]{\"G\"}=".$Istate[$i]{"G"}."\n";
print "\$Istate[$i]{\"H\"}=".$Istate[$i]{"H"}."\n";
print "\$Istate[$i]{\"I\"}=".$Istate[$i]{"I"}."\n";
print "\$Istate[$i]{\"K\"}=".$Istate[$i]{"K"}."\n";
print "\$Istate[$i]{\"L\"}=".$Istate[$i]{"L"}."\n";
print "\$Istate[$i]{\"M\"}=".$Istate[$i]{"M"}."\n";
print "\$Istate[$i]{\"N\"}=".$Istate[$i]{"N"}."\n";
print "\$Istate[$i]{\"P\"}=".$Istate[$i]{"P"}."\n";
print "\$Istate[$i]{\"Q\"}=".$Istate[$i]{"Q"}."\n";
print "\$Istate[$i]{\"R\"}=".$Istate[$i]{"R"}."\n";
print "\$Istate[$i]{\"S\"}=".$Istate[$i]{"S"}."\n";
print "\$Istate[$i]{\"T\"}=".$Istate[$i]{"T"}."\n";
print "\$Istate[$i]{\"V\"}=".$Istate[$i]{"V"}."\n";
print "\$Istate[$i]{\"W\"}=".$Istate[$i]{"W"}."\n";
print "\$Istate[$i]{\"X\"}=".$Istate[$i]{"X"}."\n";
print "\$Istate[$i]{\"Y\"}=".$Istate[$i]{"Y"}."\n";
print "\n";
}
=cut


@AA=qw(A C D E F G H I K L M N P Q R S T V W Y);
$m=0;
for($i=0;$i<scalar(@AA);$i++){
   for($j=0;$j<scalar(@AA);$j++){
      for($k=0;$k<scalar(@AA);$k++){
         for($l=0;$l<scalar(@AA);$l++){
             $randpep[$m]{'name'}=$AA[$i].$AA[$j].$AA[$k].$AA[$l];
             $randpep[$m]{'prob'}=forward($randpep[$m]{'name'},"M8");
             print $randpep[$m]{'name'}."\t".$randpep[$m]{'prob'}."\n";
             $m+=1;
         }
      }
   }
}


open(ngAIP,"./nogapsAIP.txt")||die "$!";
@contentng=<ngAIP>;
close(ngAIP);

print forward(chomp($contentng[$i]),"M8");
for($i=0;$i<scalar(@contentng);$i++){
    chomp($contentng[$i]);
    print $contentng[$i]."\t\t\t\t\t\t".forward($contentng[$i],"M8")."\n";
}



#print forward("RIQRGPGRAFVTIGK","M8");

# the forware algorithm of profile hidden markov model
# the parameters of this function are observed sequece string and the t+1 state
# ( it is a socalled previous state,actually it is the later state)
sub forward{
    local($tmpchr,$sum,$i,@emitstate,$j,$nextx);
    $j=0;
    $sum=0;
    $tmpchr=substr($_[0],length($_[0])-1,1);
    if (substr($_[1],0,1) eq "M"){
        $nextx="next1";
    }
    elsif(substr($_[1],0,1) eq "D"){
        $nextx="next2";
    }
    else{
        $nextx="next3";
    }

    for ($i=0;$i<=$mstate;$i++){
         if ($Mstate[$i]{$tmpchr}!=0){
             if (substr($_[1],0,1) eq "M"){
                if (substr($_[1],1,length($_[1])-1)-substr($Mstate[$i]{'name'},1,length($Mstate[$i]{'name'})-1)>0){
                    $emitstate[$j]=$Mstate[$i];
                    $j+=1;
                }
             }
             else{
                if (substr($_[1],1,length($_[1])-1)-substr($Mstate[$i]{'name'},1,length($Mstate[$i]{'name'})-1)>=0){
                    $emitstate[$j]=$Mstate[$i];
                    $j+=1;
                }
             }
         }
         if ($Istate[$i]{$tmpchr}!=0){
             if (substr($_[1],0,1) eq "M"){
                if (substr($_[1],1,length($_[1])-1)-substr($Istate[$i]{'name'},1,length($Istate[$i]{'name'})-1)>0){
                    $emitstate[$j]=$Istate[$i];
                    $j+=1;
                }
             }
             else{
                if (substr($_[1],1,length($_[1])-1)-substr($Istate[$i]{'name'},1,length($Istate[$i]{'name'})-1)>=0){
                    $emitstate[$j]=$Istate[$i];
                    $j+=1;
                }
             }
         }
    }
    if (length($_[0])==1){
        for($i=0;$i<$j;$i++){
            $sum+=trans_path_prob("M0",$emitstate[$i]{'name'})*$emitstate[$i]{$tmpchr}*trans_path_prob($emitstate[$i]{'name'},$_[1]);
            #print "$emitstate[$i]{'name'}$_[1]\t".trans_path_prob("M0",$emitstate[$i]{'name'})."\t".$emitstate[$i]{$tmpchr}."\t".trans_path_prob($emitstate[$i]{'name'},$_[1])."\t";
            #print trans_path_prob("M0",$emitstate[$i]{'name'})*$emitstate[$i]{$tmpchr}*trans_path_prob($emitstate[$i]{'name'},$_[1])."\n";
         }
         #print "$sum\n";
         return $sum;
    }
    else{
         for($i=0;$i<$j;$i++){
             $sum+=forward(substr($_[0],0,length($_[0])-1),$emitstate[$i]{'name'})*$emitstate[$i]{$tmpchr}*trans_path_prob($emitstate[$i]{'name'},$_[1]);
         }
    }
    return($sum);
}

#print trans_path_prob("M7","M7")."\n";
# 计算current state 经过D state到 previous state之间的转移概率
# trans_path_prob(curState,preState)
sub trans_path_prob{
    local($i,$prob);
    $i=substr($_[1],1,length($_[1])-1)-substr($_[0],1,length($_[0])-1);
    return 0 if ($_[1] eq $_[0] and substr($_[0],0,1) eq "M");
    if ($i>1){
            if (substr($_[0],0,1) eq "M"){
               $prob=$Mstate[substr($_[0],1,length($_[0])-1)]{'next2'}[1];
            }
            else{
               $prob=$Istate[substr($_[0],1,length($_[0])-1)]{'next2'}[1];
            }
            while($i>2){
               $prob*=$Dstate[substr($_[1],1,length($_[1])-1)-($i-1)]{'next2'}[1];
               $i-=1;
            }
            if (substr($_[1],0,1) eq "M"){
               $prob*=$Dstate[substr($_[1],1,length($_[1])-1)-1]{'next1'}[1];
            }
            else{
               $prob*=$Dstate[substr($_[1],1,length($_[1])-1)-1]{'next1'}[1];
               $prob*=$Dstate[substr($_[1],1,length($_[1])-1)]{'next3'}[1];
            }
    }
    elsif($i==1){
            if (substr($_[1],0,1) eq "M"){
               if (substr($_[0],0,1) eq "M"){
                  $prob=$Mstate[substr($_[0],1,length($_[0])-1)]{'next1'}[1];
               }
               else{
                  $prob=$Istate[substr($_[0],1,length($_[0])-1)]{'next1'}[1];
               }
            }
            else{
               if (substr($_[0],0,1) eq "M"){
                  $prob=$Mstate[substr($_[0],1,length($_[0])-1)]{'next2'}[1];
                  $prob*=$Dstate[substr($_[0],1,length($_[0])-1)+1]{'next3'}[1];
               }
               else{
                  $prob=$Istate[substr($_[0],1,length($_[0])-1)]{'next2'}[1];
                  $prob*=$Dstate[substr($_[0],1,length($_[0])-1)+1]{'next3'}[1];
               }
            }
    }
    elsif($i==0){ # 如果$i=0则肯定为MiIi或IiIi
            if ($_[1] eq $_[0]){
               $prob=$Istate[substr($_[0],1,length($_[0])-1)]{'next3'}[1];
            }
            else{
               $prob=$Mstate[substr($_[0],1,length($_[0])-1)]{'next3'}[1];
            }
    }
    else{
           return 0;
    }
    return ($prob);
}



# 统计一条状态序列中的状态转移
# campstate($seqstate[$i]) 输入状态序列
sub campstate{
        local($i,$n,$tempchr,$prestate,$j);
        $prestate=\@Mstate;
        $j=0;
        $n=2;
        $i=0;
        while ($i<length($_[0])){

              if ($i+$n<length($_[0])){
                 if (substr($_[0],$i+$n,1) ne "M" and substr($_[0],$i+$n,1) ne "I" and substr($_[0],$i+$n,1) ne "D"){
                    $n+=1;
                 }
              }
              $tempchr=substr($_[0],$i,$n);
              $i+=$n;
                if (substr($tempchr,0,1) eq "M"){
                    $$prestate[$j]{'next1'}[1]+=1;
                }
                elsif(substr($tempchr,0,1) eq "D"){
                    $$prestate[$j]{'next2'}[1]+=1;
                }
                else{
                    $$prestate[$j]{'next3'}[1]+=1;
                }

                if (substr($tempchr,0,1) eq "M"){
                    $prestate=\@Mstate;
                }
                elsif(substr($tempchr,0,1) eq "D"){
                    $prestate=\@Dstate;
                }
                else{
                    $prestate=\@Istate;
                }
                $j=substr($tempchr,1,length($tempchr)-1);
        }
        $$prestate[$j]{'next1'}[1]+=1;
    
}


# 将对齐后的各序列转化成其对应的状态序列的函数
# 算法：如果序列中的某个位置不为"-"则在该位置加入相应的状态；
#　否则看该位置是否为match state，如果是加入D和状态编号，否则不进行任何操作
# str2state（$content[$i]）输入对齐后的序列
sub str2state{
    local($temp);
    local($i);
    local($tmpstate);
    $tmpstate="";
    for($i=0;$i<length($_[0])-1;$i++){
        $temp=substr($_[0],$i,1);
        if ($temp ne "-"){
            $tmpstate.=$state[$i];
        }
        else{
             if (substr($state[$i],0,1) eq "M"){
                  $tmpstate=$tmpstate."D".substr($state[$i],1,length($state[$i])-1);
             }
        }
    }
    return($tmpstate);
}