Blat percent identity and score calculations

最新推荐文章于 2019-10-21 16:22:57 发布

keying0520

最新推荐文章于 2019-10-21 16:22:57 发布

阅读量1.3k

点赞数

文章标签： function struct licensing alignment browser file

本文链接：https://blog.csdn.net/keying0520/article/details/6476453

版权

Replicating web-based Blat percent identity and score calculations

Question:
"Using my own command-line Blat server, how can I replicate the percent identity and score calculations produced by web-based Blat?"

Response:
There isn't an option to command-line Blat that gives you the percent ID and the score. Instead, you will have to write your own program to produce the calculations, incorporating some of the functions from the Genome Browser source code.

To calculate the percent ID, incorporate the following code and function into a program that processes your Blat PSL output. The parameter isMrna should be set to TRUE, regardless of whether the input sequence is mRNA or protein.

The percent identity score is calculated like this:

100.0 - pslCalcMilliBad(psl, TRUE) * 0.1

Here is the source for pslCalcMilliBad:

int pslCalcMilliBad(struct psl *psl, boolean isMrna) /* Calculate badness in parts per thousand. */ { int sizeMul = pslIsProtein(psl) ? 3 : 1; int qAliSize, tAliSize, aliSize; int milliBad = 0; int sizeDif; int insertFactor; int total; qAliSize = sizeMul * (psl->qEnd - psl->qStart); tAliSize = psl->tEnd - psl->tStart; aliSize = min(qAliSize, tAliSize); if (aliSize <= 0) return 0; sizeDif = qAliSize - tAliSize; if (sizeDif < 0) { if (isMrna) sizeDif = 0; else sizeDif = -sizeDif; } insertFactor = psl->qNumInsert; if (!isMrna) insertFactor += psl->tNumInsert; total = (sizeMul * (psl->match + psl->repMatch + psl->misMatch)); if (total != 0) milliBad = (1000 * (psl->misMatch*sizeMul + insertFactor + round(3*log(1+sizeDif)))) / total; return milliBad; }

The complexity in milliBad arises primarily from how it handles inserts. Ignoring the inserts, the calculation is simply mismatches expressed as parts per thousand. However, the algorithm factors in insertion penalties as well, which are relatively weak compared to say blasts but still present. When huge inserts are allowed (which is necessary to accommodate introns), it is typically necessary to resort to logarithms like this calculation does.

The pslIsProtein function called by pslCalcMilliBad is:

boolean pslIsProtein(const struct psl *psl) /* is psl a protein psl (are it's blockSizes and scores in protein space) */ { int lastBlock = psl->blockCount - 1; return (((psl->strand[1] == '+' ) && (psl->tEnd == psl->tStarts[lastBlock] + 3*psl->blockSizes[lastBlock])) || ((psl->strand[1] == '-') && (psl->tStart == (psl->tSize-(psl->tStarts[lastBlock] + 3*psl->blockSizes[lastBlock]))))); }

This function automatically determines whether or not the PSL output file contains alignment information for a protein query. Alternatively, you could write the program such that the user specifies if the query is a protein or not.

The score calculation is generated by the following function:

int pslScore(const struct psl *psl) /* Return score for psl. */ { int sizeMul = pslIsProtein(psl) ? 3 : 1; return sizeMul * (psl->match + ( psl->repMatch>>1)) - sizeMul * psl->misMatch - psl->qNumInsert - psl->tNumInsert; }

For help with creating a C program to do perform these calculations, you may want to use the libraries from the Genome Browser source code. See our FAQ on source code licensing and downloads for information on obtaining the source. The file kent/src/lib/psl.c contains the pslCalcMilliBad, pslIsProtein and pslScore functions and also a useful function called pslLoadAll that loads the psl file into a linked list structure. The definition of the psl struct can be found in kent/src/inc/psl.h.

perl:

#!/usr/bin/perl
use strict;
open(IN,"try.psl") or die "Cannot open file gi";
open(OUT, "> out") or die "Cannot create file!";

while(<IN>) {
        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 = @_;

         # sizeMul depens of dna/Prot
         my $sizeMul = 1;                  #核酸等于1，蛋白质等于3

    # 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;
         }

         # insert Factor
         my $insertFactor = $cols[4];

# $insertFactor += $cols[6]; # 若为蛋白质

         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);