BLAST_TABULAR_OUTPUT Generator


Few MPI versions do not support the very useful -m 8 tabular output format. The following code can be used to parse the raw blast output obtained from mpi_blast to tabular output file. The creditials of original creator and modifications are given in comments section of the code. 
#=========================================================================================
# Modified BLAST Parser  -  Script for parsing the BLAST output in tabular format
# Version 1.1.6 (November 26, 2015)
# Modified by Arun Prasanna (arunprasanna83@gmail.com) from orginal script by Kirill Kryukov (http://kirill-kryukov.com/kirr/). He also has many other useful utilities..Check out !!
#  
# 
# Purpose: Few versions of MPI-BLAST do not support tabular output or often the user faces 
# error while opting for tabular output. Hence the original code is modified to produce 
# tabular Output verbatim -m 8 option. The output file is compliant with other clustering 
# programs like Silix
# 
# Changes made are:
#  1. Added one more regex in $qnaming section to account (24,300 letters) instances 
#  2. Added $per_id, $rounded variable to calculate percentage identity
#  2. Modified print statement to convert in tabular format
#
#=========================================================================================
# BLAST Parser  -  Script for parsing the BLAST output
#
# Version 1.1.5 (June 27, 2012)
#
# Copyright (c) 2011-2012 Kirill Kryukov
#
# This software is provided 'as-is', without any express or implied
# warranty. In no event will the authors be held liable for any damages
# arising from the use of this software.
#
# Permission is granted to anyone to use this software for any purpose,
# including commercial applications, and to alter it and redistribute it
# freely, subject to the following restrictions:
#
# 1. The origin of this software must not be misrepresented; you must not
#    claim that you wrote the original software. If you use this software
#    in a product, an acknowledgment in the product documentation would be
#     appreciated but is not required.
# 2. Altered source versions must be plainly marked as such, and must not be
#    misrepresented as being the original software.
# 3. This notice may not be removed or altered from any source distribution.
#
# This script is a BLAST parser - it will process the BLAST output and convert
# it into a compact form where each hit is represented by a single line.
#
# Usage:
#   perl blast_parser.pl <blastout.txt >parsed.txt
#
# You can also use it to parse the BLAST output on the fly as it is generated,
# saving a lot of disk space (although losing all alignments).
#========================================================================================
use strict;
my ($qname,$sname) = ('','');
my ($qchanged,$schanged,$hchanged) = (0,0,0);
my ($qnaming,$snaming,$qlen,$slen) = (0,0,0,0);
my ($qstart,$qend,$sstart,$send) = (-1,-1,-1,-1);
my ($bits,$expect,$length,$ident,$positives,$gaps,$frame,$per_id) = (-1,'',-1,-1,-1,-1,'',''); #added per_id with '' 
my ($n_queries,$n_hits) = (0,0);

while (<>)
{
    chomp;

    if ($qnaming)
    {
        if ($_ =~ /^\s*Length=(\d+)$/ or $_ =~ /^\s*\((\d+) letters\)$/ or $_ =~ /^\s*\((\d+,\d+) letters\)$/) #added third or statement to check for comma
        {
            $qlen = $1;
            $qnaming = 0;
            #print STDERR "Query = $qname\n";
            #print STDERR ">";
        }
        else { $qname .= " $_"; }
        next;
    }
    #print "The value in qlen is: $qlen\n";
    #print "The value in qname is: $qname\n";

    if ($snaming)
    {
        if (/^\s*Length\s*=\s*(\d+)$/)
        {
            $slen = $1;
            $snaming = 0;
            #print STDERR "    Subject = $sname\n";
            #print STDERR ":";
        }
        else { s/^\s+//; $sname .= " $_"; }
        next;
    }
    #print "The value in slen is: $slen\n";
    #print "The value in sname is: $sname\n";

    
    if (/^Query[\:\s]\s(\d+)\s.*\s(\d+)$/)   #Changed if to while to check if it can still look for fragments
    {
        if ($qstart < 0 or $1 < $qstart) { $qstart = $1; }
        if ($qend < 0 or $1 > $qend) { $qend = $1; }
        if ($2 < $qstart) { $qstart = $2; }
        if ($2 > $qend) { $qend = $2; }
        next;
    }
    #print "The value in qstart is: $qstart\n";
    #print "The value in qend is: $qend\n";
    
    if (/^Sbjct[\:\s]\s(\d+)\s.*\s(\d+)$/)    #Changed if to while to check if it can still look for fragments
    {
        if ($sstart < 0 or $1 < $sstart) { $sstart = $1; }
        if ($send < 0 or $1 > $send) { $send = $1; }
        if ($2 < $sstart) { $sstart = $2; }
        if ($2 > $send) { $send = $2; }
        next;
    }
    #print "The value in sstart is: $sstart\n";
    #print "The value in send is: $send\n";
    
    if (/^\sScore/)
    {
        if ($snaming) { print STDERR "Error: can't complete subject name: $sname\n"; next; }
        if ($qnaming) { print STDERR "Error: can't complete query name: $qname\n"; next; }
        #&flush_hit();   #not letting it flush hit
        ($qstart,$qend,$sstart,$send) = (-1,-1,-1,-1);
        ($bits,$expect,$length,$ident,$positives,$gaps,$frame) = (-1,'',-1,-1,-1,-1,'');

        if (/^\sScore\s+=\s+([\d\+e\.]+)\s+bits/) { $bits = int($1); } else { print STDERR "Error: Can't parse score: \"$_\"\n"; }
        if (/\s+Expect\(?\d*\)?\s+=\s+([e\-0-9\.]+)/) { $expect = $1; } else { print STDERR "Error: Can't parse expectation: \"$_\"\n"; }
        $hchanged = 1;
        next;
    }

    if (/^\sIdentities/)
    {
        ($length,$ident,$positives,$gaps,$frame) = (-1,-1,-1,-1,'');
        if (/^\sIdentities\s+=\s+(\d+)\/(\d+)\s+\(\d+\%\)/)
        {
            ($ident,$length) = ($1,$2);
        }
        if (/\s+Positives\s+=\s+(\d+)\/\d+\s+\(\d+\%\)/)
        {
            ($positives) = ($1);
        }
        if (/Gaps\s+=\s+(\d+)\//)
        {
            $gaps = $1;
        }
        else
        {
            $gaps = 0;
        }
        if ($length < 0)
        {
            print STDERR "Error: unexpected BLAST ' Identities = ' line format: $_\n";
        }
        next;
    }

    if (/\s[Ff]rame\s=\s(\S+)$/)
    {
        $frame = $1;
        if ($frame !~ /^(-?\+?\d|-?\+?\d\/-?\+?\d)$/) { print STDERR "Error: Unexpected frame: \"$frame\"\n"; }
        next;
    }

    if (/\sStrand\s=\s(\S+)$/)
    {
        $frame = $1;
        if ($frame !~ /^Plus\/(Plus|Minus)$/) { print STDERR "Error: Unexpected strand: \"$frame\"\n"; }
        next;
    }

    if (/^Query=\s+(.+)$/)
    {
        my $new_qname = $1;
        &flush_hit();
        $qname = $new_qname;
        $qnaming = 1;
        $qchanged = 1;
        $n_queries++;
        next;
    }

    if (/^>(.+)$/)
    {
        my $new_sname = $1;
        &flush_hit();
        $sname = $new_sname;
        $snaming = 1;
        $schanged = 1;
        next;
    }
}

&flush_hit();
print STDERR "OK, processed $n_queries queries, $n_hits hits\n";


sub flush_hit
{
    unless ($hchanged) { return; }
    if ($qstart<0 || $qend<0 || $sstart<0 || $send<0)
    {
        print STDERR "Error: no coordinates!\n";
        $hchanged = 0;
        return;
    }
    if ($bits<0 || $expect eq '' || $length<0 || $ident<0)
    {
        print STDERR "Error: parameters missing!\n";
        $hchanged = 0;
        return;
    }
    if ($frame eq '') { $frame = '-'; }

    #if ($frame eq '')
    #{
    #    print STDERR "Error: Frame/Strand is missing!\n";
    #    $hchanged = 0;
    #    return;
    #}

    if ($positives < 0) { $positives = '-'; }
    if ($gaps < 0) { $gaps = '-'; }

    if ($qchanged)
    {
        $qname =~ s/\s+$//;
        #print "\n$qname\n";  #print "\n$qname ($qlen)\n";
        $qchanged = 0;
        $schanged = 1;
    }

    if ($schanged)
    {
        $sname =~ s/\s{2,}/ /g;
        $sname =~ s/\s+$//;
        $sname =~ s/^\s+//;
        $sname =~ s/^lcl\|//;
        #print "\t$sname\t"; #print "     $sname ($slen)\n";
        $schanged = 0;
    }
    
    $per_id = ($ident/$length)*100;
    my $rounded = sprintf "%.2f", $per_id;

    print "$qname\t";
    print "$sname\t";
    print "$rounded\t";
    print "$length\t";
    print "$positives\t";
    print "$gaps\t";
    print "$qstart\t";
    print "$qend\t";
    print "$sstart\t";
    print "$send\t";
    print "$expect\t"; 
    print "$bits\n";
    $hchanged = 0;
    $frame = '';

    $n_hits++;
}

Comments

Post a Comment

Popular posts from this blog

Pick Matching lines with list of keywords

#Simple code to find the occurrence of list of search terms in single line in a huge file. Search_terms = [ 'A' , 'B' , 'C' , 'D' ] with open ( 'BigFile.txt' , 'r' ) as infile : entries = infile . read () each_line = entries . splitlines () new_list = [] for ix , row in enumerate ( each_line ): element = row . split ( "\t" ) if set ( Search_terms ) == set ( element ): #Use <= if you want no strict option new_list . append ( element ) print ix + 1 #List the matching line number ! thefile = open ( 'Output.txt' , 'w' ) for item in new_list : print >> thefile , item
Read more

Install Parallel versions of Python from source

Few programs require certain versions of python. Especially, when you do not have root permission in your unix machine, you can still install a python version in your /home/usr directory. You can run it in parallel to in-built version. Follow these steps: 1. Unpack with  tar -xvf python-x.x.tar 2. ./configure 3. make altinstall prefix=~ exec-prefix=~ (./lib and ./bin will be created in your home directory(~)) 4. create alias for python-x.x => i.e cd ~/bin/ > ln -s python-x.x python ! If you create this outside bin directory, you will get a warning (ln: symbolic link already exists !! remember it refers to inbuilt python) 5. Complete alias creation by editing: vi ~/.bashrc with alias python='~/bin/python' Now any software you try to install with python setup.py it will access the version you installed in your home ! You can still install softwares using your version by skipping steps 4 & 5. But, everytime you need to keep specifying the export PATH ...
Read more

Map multiple annotations using pandas

A simple pandas solution to map multiple annotations for a protein.  A protein or gene file will have annotations curated by different methods. Most frequently, biologists will encounter more than one annotation for a single protein. It is a task in itself to pick the right annotation.  One of the simple ways is to consolidate them and pick the right ones after enough evidence is known.  Coding i n matlab or other languages may require more number of lines to achieve the same output. Here, a simple 'groupby' of pandas can do produce the same outputin seconds ! # Map multiple annotation for protein ID and join with a delimiter ''' Sample input A_xx Annotation1 A_xx Annotation2 B_xx Annotation1 Sample output A_xx Annotation1, Annotation2 B_xx Annotation1 ''' import pandas as pd data = pd . read_csv( 'input_file.txt' , delimiter = ' \t ' ) dfc = data . groupby([ 'PROTID' ])[ 'ANNOTATION...
Read more