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

Fasta_Dupicate_Header

A simple, self-explanatory matlab code to identify duplicate headers in fasta files. 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 %Simple matlab code to check for the duplicate header in fasta files %Store the size of header -> unique(header) ->size of new header %Copy the Table data in excel and compare the two values clear; clc; tic; Path = 'Drive\Path\FileName' ; % FileList = dir(Path); [rFL, cFL] = size(FileList); for i = 3:rFL %i of 1 & 2 are . & .. respectively     Fas_Fname{i-2,1} = FileList(i).name; %FileList is a structure end [rFas,cFas] = size(Fas_Fname); for i = 1:rFas     clear Header Seq Old_Header Unik_header     OpenFile = cell2mat(strcat(Path,Fas_Fname(i)));     [Header, Seq] = fastaread(OpenFile);[rH,cH] = size(Header);     Old_Header = length(Header);     Unik_header = length(unique(Header));     Table{i,1} = Fas_Fname(i);     Table{i,2} = num2str(Old_Header);     Table{i,3} = num2str(Unik_heade
Read more

Gene Copy Number Matrix

Given a cluster file, one can create a gene copy number matrix (GCN). With this self-explanatory simple matlab file it is easy to create one. Input Format: Cluster file.xlsx: (1) A  B  C  D (2) A  A  A (3) B C (4) D D A List File.xlsx: A B C D Output: (of course, the output will have the file with only numbers printed)      A B C D (1) 1  1  1  1 (2) 3  0  0  0 (3) 0  1  1  0 (4) 1  0  0  2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 %Author = Arun Prasanna %Create a gene copy number matrix from cluster information clear ; clc ; tic [ mat1 , mat ] = xlsread ( 'ClusterFile.xlsx' , 'Sheet1' ); clear mat1 [ mat2 , head ] = xlsread ( 'Organism_list.xlsx' , 'Sheet1' ); clear mat2 %species/gene name new_head = head (:, 1 ) ' ; %transpose to make it as header [ rmat , cmat ] = size ( mat ); [ rhead , chead ] = size ( new_head ); out = zeros ( rmat , chead );
Read more

Fasta_Header_Rename

A simple matlab code to rename the headers in fasta file. Self-explanatory variable names. 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 %Author = Arun Prasanna %Rename the headers in fasta file to desired choice %For example: The input fasta file used here had header in >num_name format %strtok is used to strip and extract the required format clear; clc; tic; Path = 'Drive\Path\ToReadFile' ; % FileList = dir(Path); [rFL, cFL] = size(FileList); for i = 3:rFL %i of 1 & 2 are . & .. respectively     Fas_Fname{i-2,1} = FileList(i).name; %FileList is a structure end [rFas,cFas] = size(Fas_Fname); for i = 1:rFas     clear Header Seq ProtID Sp new_Header     OpenFile = cell2mat(strcat(Path,Fas_Fname(i)));     [Header, Seq] = fastaread(OpenFile);[rH,cH] = size(Header);     for j = 1:cH         [ProtID, Sp] = strtok(Header(1,j), '_' ); %First Sp = _name         [Sp, rm] = strtok(Sp, '_' ); %Se
Read more