Info | Home

BioPHP - Alignment of DNA

Original code submitted by joseba
Code bellow is covered by GNU GPL v2 license.

Description

Last change: 2010/10/18 17:04 | Edit description | Recent Changes | Original description

Alignment of two DNA sequences. 
Smith-Waterman alignment method is used.

Code

Last change: 2010/10/18 17:04 | Edit Code | Recent Changes | Download | Original code

function alignment_of_dna($seqa,$seqb){
        $match = 2;
        $mismatch = -1;
        $gap = -4;

        $a = preg_split('//', $seqa, -1, PREG_SPLIT_NO_EMPTY);
        $b = preg_split('//', $seqb, -1, PREG_SPLIT_NO_EMPTY);
        $maxa=sizeof($a);
        $maxb=sizeof($b);
        $lenn=max ($maxa,$maxb);

        // Creación de la matriz
        // He reducido el código para hacerlo mas simple y rapido, pero tan solo ahorra un 20% del tiempo
        // Con matrices muy grandes, PHP no sabe trabajar muy bien (es poco eficaz).
        $mx=0;
        for ($i=0;$i<$maxa;$i++){
                for ($j=0;$j<$maxb;$j++){
                        if($b[$j]==$a[$i]){
                                $x=$matriz[$j-1][$i-1]+$match;
                        }else{
                                $x=max (0,$matriz[$j-1][$i-1]-1,$matriz[$j][$i-1]-4,$matriz[$j-1][$i]-4);
                        }
                        $matriz[$j][$i]=$x;
                        if ($mx<$x){$mx=$x; $mj=$j; $mi=$i;}
                }
        }
        // Matriz terminada

        $j=$mj;
        $i=$mi;

        $matrizz[$j][$i]=1;
        while ($i>0 or $j>0):
                $aa=$matriz[$j-1][$i-1];
                $ab=$matriz[$j][$i-1];
                $ac=$matriz[$j-1][$i];
                if($aa<>'//' or $aa==0){
                        if($aa>=$ab and $aa>=$ac){
                                $j=$j-1;
                                $i=$i-1;
                        }
                        if($ab>$aa){$i=$i-1;}
                        if($ac>$aa){$j=$j-1;}
                }else{
                        if($ab<>'//'){$i=$i-1;}
                        if($ac<>'//'){$j=$j-1;}
                }
                if($j<0){$j=0;}
                if($i<0){$i=0;}
                $matrizz[$j][$i]=1;
        endwhile;

        $j=$mj;
        $i=$mi;

        while ($i<strlen($seqa)-1 or $j<strlen($seqb)-1):
                $aa=$matriz[$j+1][$i+1];
                $ab=$matriz[$j][$i+1];
                $ac=$matriz[$j+1][$i];
                if($aa<>'//'){
                        if($aa>=$ab and $aa>=$ac){
                                $j=$j+1;
                                $i=$i+1;
                        }
                        if($ab>$aa){$i=$i+1;}
                        if($ac>$aa){$j=$j+1;}
                }else{
                        if($ab<>'//'){$i=$i+1;}
                        if($ac<>'//'){$j=$j+1;}
                }
                if($j>$lenn){$j=$lenn;}
                if($i>$lenn){$i=$lenn;}
                $matrizz[$j][$i]=1;
        endwhile;

        $j=0;
        $i=0;
        $t=1;

        while ($i<strlen($seqa)-2 and $j<strlen($seqb)-2 and $t=1):
                $t=0;
                if($matrizz[$j+1][$i+1]==1){
                        $t=1;
                        $sseqa.=$a[$i];
                  $sseqb.=$b[$j];
                  $i=$i+1;
                  $j=$j+1;
                }
                if($matrizz[$j][$i+1]==1){
                        $t=1;
                        $sseqa.=$a[$i];
                        $sseqb.="-";
                        $i=$i+1;
                }
                if($matrizz[$j+1][$i]==1){
                        $t=1;
                        $sseqa.="-";
                        $sseqb.=$b[$j];
                        $j=$j+1;
                }
        endwhile;

        if($matrizz[$j+1][$i+1]==1){
                $sseqa.=$a[$i];
                $sseqb.=$b[$j];
                $i=$i+1;
                $j=$j+1;
                $t=1;
        }
        if($t==0 and $matrizz[$j][$i+1]==1){
                $sseqa.=$a[$i];
                $sseqb.="-";
                $i=$i+1;
                }
        if($t==0 and $matrizz[$j+1][$i]==1){
                $sseqa.="-";
                $sseqb.=$b[$j];
                $j=$j+1;
        }
        if($i+1==$maxa){
                for ($ii=$j;$ii<$maxb;$ii++){$sseqb.=$b[$ii];}
                $sseqa.=$a[$i];
                for ($ii=$j;$ii<$maxb-1;$ii++){$sseqa.="-";}
        }
        if($j+1==$maxb){
                for ($ii=$i;$ii<$maxa;$ii++){$sseqa.=$a[$ii];}
                $sseqb.=$b[$j];
                for ($ii=$i;$ii<$maxa-1;$ii++){$sseqb.="-";}
        }


        // tengo que quitar la última letra del alineamiento, de lo contrario, se repite la ultima posición
        // por que se da la repetición?
        $results["seqa"]=substr($sseqa,0,strlen($sseqa)-1);
        $results["seqb"]=substr($sseqb,0,strlen($sseqb)-1);
        //$results["seqa"]=$sseqa;
        //$results["seqb"]=$sseqb;


        return $results;
}