[Bioperl-l] script to calculate percentage identity for BLAT psl

gopu_36 gopu_36 at yahoo.com
Fri Mar 13 06:55:14 EDT 2009


Hi 
As given in that page I tried as below as one of the script. But I always
get 100% identity. Please let me know the problem in my following code.
Since I don't use any options and am also mapping DNA sequences, I have made
following changes in the code: Please let me know the problem.

1
 my $sizeMul = 1; and commented:
#if ($option{p}) {
         #        $sizeMul = 3;
         #} else {
         #        $sizeMul = 1;
         #}
2
if ($sizeDiff < 0) {
         	$sizeDiff = 0;
                 #if ($option{m}) {
                         #$sizeDiff = 0;
                 #} else {
                         #$sizeDiff = -($sizeDiff);
                 #}
         }

=================================================== 
And perl script is as bleow: Thanks.
 

#!/usr/bin/perl

while(<>) {
        chomp $_;
        my @v = split(/\t/,$_);
        get_pid($_);

}
sub get_pid {
         my @line = @_;
         my $pid = (100.0 - (&pslCalcMilliBad(@line) * 0.1));
	 print "The percentage: $pid\n";
         #return $pid;
}

sub pslCalcMilliBad {
         my @cols = @_;

         # sizeNul depens of dna/Prot
         my $sizeMul = 1;
         #if ($option{p}) {
         #        $sizeMul = 3;
         #} else {
         #        $sizeMul = 1;
         #}

         # cols[0]  matches
         # cols[1]  misMatches
         # cols[2]  repMaches
         # cols[4]  qNumInsert
         # cols[6]  tNumInsert
         # cols[11] qStart
         # cols[12] qEnd
         # cols[15] tStart
         # cols[16] tEnd

         my $qAliSize = $sizeMul * ($cols[12] - $cols[11]);
         my $tAliSize = $cols[16] - $cols[15];

         # I want the minimum of qAliSize and tAliSize
         my $aliSize;
         $qAliSize < $tAliSize ? $aliSize = $qAliSize : $aliSize =  
$tAliSize;

         # return 0 is AliSize == 0
         return 0 if ($aliSize <= 0);

         # size diff
         my $sizeDiff = $qAliSize - $tAliSize;
         if ($sizeDiff < 0) {
         	$sizeDiff = 0;
                 #if ($option{m}) {
                         #$sizeDiff = 0;
                 #} else {
                         #$sizeDiff = -($sizeDiff);
                 #}
         }

         # insert Factor
         my $insertFactor = $cols[4];
         $insertFactor += $cols[6] unless ($option{m});
         my $milliBad = (1000 * ($cols[1]*$sizeMul + $insertFactor +  
&round(3*log( 1 + $sizeDiff)))) / ($sizeMul * ($cols[0] + $cols[2]
+ $cols[1]));
	        return $milliBad;

}

sub round {
         my $number = shift;
         return int($number + .5);
}



gopu_36 wrote:
> 
> Hi,
> 
> I did go through FAQ from BLAT on how to calculate the precentage identity
> from http://genome.ucsc.edu/FAQ/FAQblat#blat4
> As a new comer, I don;t usederstand on how to implement this. Please let
> me know how to plugin the script for my output.psl file. Please let me
> know. It would be of great help.
> 
> Thanks and Regards.
> 

-- 
View this message in context: http://www.nabble.com/script-to-calculate-percentage-identity-for-BLAT-psl-tp22476983p22494163.html
Sent from the Perl - Bioperl-L mailing list archive at Nabble.com.



More information about the Bioperl-l mailing list