[Bio] / FigKernelScripts / blastn_chunks.pl Repository:
ViewVC logotype

View of /FigKernelScripts/blastn_chunks.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (download) (as text) (annotate)
Mon Feb 5 03:12:11 2007 UTC (12 years, 9 months ago) by overbeek
Branch: MAIN
CVS Tags: mgrast_dev_08112011, rast_rel_2009_05_18, mgrast_dev_08022011, rast_rel_2014_0912, rast_rel_2008_06_18, myrast_rel40, rast_rel_2008_06_16, mgrast_dev_05262011, rast_rel_2008_12_18, mgrast_dev_04082011, rast_rel_2008_07_21, rast_rel_2010_0928, rast_2008_0924, mgrast_version_3_2, mgrast_dev_12152011, rast_rel_2008_04_23, mgrast_dev_06072011, rast_rel_2008_09_30, rast_rel_2009_0925, rast_rel_2010_0526, rast_rel_2014_0729, mgrast_dev_02212011, rast_rel_2010_1206, mgrast_release_3_0, mgrast_dev_03252011, rast_rel_2010_0118, mgrast_rel_2008_0924, mgrast_rel_2008_1110_v2, rast_rel_2009_02_05, rast_rel_2011_0119, mgrast_rel_2008_0625, mgrast_release_3_0_4, mgrast_release_3_0_2, mgrast_release_3_0_3, mgrast_release_3_0_1, mgrast_dev_03312011, mgrast_release_3_1_2, mgrast_release_3_1_1, mgrast_release_3_1_0, mgrast_dev_04132011, rast_rel_2008_10_09, mgrast_dev_04012011, rast_release_2008_09_29, mgrast_rel_2008_0806, mgrast_rel_2008_0923, mgrast_rel_2008_0919, rast_rel_2009_07_09, rast_rel_2010_0827, mgrast_rel_2008_1110, myrast_33, rast_rel_2011_0928, rast_rel_2008_09_29, mgrast_rel_2008_0917, rast_rel_2008_10_29, mgrast_dev_04052011, mgrast_dev_02222011, rast_rel_2009_03_26, mgrast_dev_10262011, rast_rel_2008_11_24, rast_rel_2008_08_07, HEAD
Changes since 1.2: +9 -3 lines
initial, undebugged version of edit_feature.cgi

use FIG;

$usage = "usage: blastn_chunks File DB ChunkSz MaxSc MinIden MinLen [BlastParms]";

(
 ($file        = shift @ARGV) &&
 ($db          = shift @ARGV) &&
 ($chunksz     = shift @ARGV) &&
 ($max_sc      = shift @ARGV) &&
 ($min_iden    = shift @ARGV) &&
 ($min_len     = shift @ARGV)
)
    || die $usage;

my $blast_parms = join(" ",@ARGV);

$dir = "$FIG_Config::temp/blast.$$";
&FIG::verify_dir($dir);
$trans1 = &make_chunked_file($file,"$dir/query_file","id1",$chunksz);
$trans2 = &make_chunked_file($db,"$dir/db_file","id2",$chunksz);
&FIG::run("$FIG_Config::ext_bin/blastall $blast_parms -m 8 -i $dir/query_file -d $dir/db_file -FF -p blastn -r 1 -q -1 > $dir/blastout");
open(IN,"<$dir/blastout") || die "could not open $dir/blastout";
open(SORT,"| sort +3n -6 > $dir/sorted.blast") || die "could not sort blast output";
while (defined($_ = <IN>))
{
    chop;
    my($id1,$id2,$iden,undef,undef,undef,$b1,$e1,$b2,$e2,$psc,$bsc) = split(/\t/,$_);

    my($x,$y);
    if (($x = $trans1->{$id1}) && ($y = $trans2->{$id2}))
    {
	my($c1,$off1) = @$x;
	my($c2,$off2) = @$y;
	$b1a          = $b1 + $off1;
	$e1a          = $e1 + $off1;
	$b2a          = $b2 + $off2;
	$e2a          = $e2 + $off2;
	$len_match    = abs($e1-$b1)+1;
	if (($c1 ne $c2) || ($b1a != $b2a) || ($e1a != $e2a))
	{
	    if (($psc < $max_sc) && ($len_match >= $min_len) && ($min_iden <= $iden))
	    {
		print SORT join("\t",($c1,$c2,$iden,$b1a,$e1a,$b2a,$e2a,$psc,$bsc,$len_match)),"\n";
	    }
	}
    }
    else
    {
	print &Dumper($_); die "aborted";
    }
}
close(SORT);

open(IN,"<$dir/sorted.blast") || die "could not open sorted.blast";
while (defined($_ = <IN>)) { print $_ }
close(IN);

system "/bin/rm -rf $dir";

sub make_chunked_file {
    my($in,$out,$pre,$chunksz) = @_;

    my $tran = {};
    open(IN,"<$in")   || die "could not open $in";
    open(OUT,">$out") || die "could not open $out";

    my $n = 1;
    $/ = "\n>";
    while (defined($_ = <IN>))
    {
	chomp;
	if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
	{
	    my $id = $1;
	    my $seq = $2;
	    $seq =~ s/\s//g;
	    my $ln = length($seq);
	    for (my $i=0; ($i < $ln); $i += $chunksz)
	    {
		my $sz = (($i + (2 * $chunksz)) < $ln) ? (2 * $chunksz) : $ln - $i;
		my $chunk = substr($seq,$i,$sz);
		print OUT ">$pre.$n\n$chunk\n";
		$tran->{"$pre.$n"} = [$id,$i];
		$n++;
	    }
	}
    }
    $/ = "\n";
    close(IN);
    close(OUT);
    &FIG::run("$FIG_Config::ext_bin/formatdb -i $out -p F");
    return $tran;
}

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3