[Bio] / FigMetagenomeTools / 16sblasts.pl Repository:
ViewVC logotype

View of /FigMetagenomeTools/16sblasts.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (download) (as text) (annotate)
Tue Apr 10 16:37:39 2007 UTC (12 years, 10 months ago) by olson
Branch: MAIN
CVS Tags: mgrast_dev_08112011, mgrast_dev_08022011, mgrast_dev_05262011, mgrast_dev_04082011, mgrast_version_3_2, mgrast_dev_12152011, mgrast_dev_06072011, mgrast_rel_2008_0806, mgrast_dev_10262011, mgrast_dev_02212011, mgrast_rel_2008_0923, mgrast_release_3_0, mgrast_dev_03252011, mgrast_rel_2008_0924, mgrast_rel_2008_1110_v2, 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, mgrast_dev_04012011, mgrast_rel_2008_0919, mgrast_rel_2008_1110, myrast_33, mgrast_rel_2008_0917, mgrast_dev_04052011, mgrast_dev_02222011, HEAD
new stuff that didn't make it in before

#!/usr/bin/perl -w

# generate a summary of the 16S blasts 

use strict;
my $usage=<<EOF;

$0 <options> <list of blast files>

-s 16S correspondance file

-c cutoff (1e-10)
-n number of hits to keep for each query (default 5)

EOF

my ($corrf, $cutoff, $num)=('',1e-10,5);
my @blastfiles;

while (@ARGV) {
 my $test=shift(@ARGV);
 if ($test eq "-s") {$corrf=shift @ARGV}
 elsif ($test eq "-c") {$cutoff=shift @ARGV}
 elsif ($test eq "-n") {$num=shift @ARGV}
 elsif (-e $test) {push @blastfiles, $test}
 else {print "What is $test\n"}
}


die $usage unless ($blastfiles[0] && $corrf);

open(IN, $corrf) || die "Can't open correspondance file $corrf";
my %corr;
while (<IN>) {
 chomp;
 my @a=split /\t/;
 $corr{$a[0]}=$a[1];
}
close IN;

my %hits;
foreach my $blastf (@blastfiles) {
 open(IN, $blastf) || die "Can't open $blastf";
 while (<IN>) {
  chomp;
  my @a=split /\t/;
  next if (defined $cutoff && $a[10] > $cutoff);
  next if (defined $hits{$a[0]} &&  scalar @{$hits{$a[0]}} > $num);
  my $res;
  if ($corr{$a[1]}) {$res = $corr{$a[1]} . "(" . $a[1] . ")"}
  else {print STDERR "No correspondance for $a[1]\n"; $res=$a[1]}
  push @{$hits{$a[0]}}, $res;
 }
}

foreach my $hit (keys %hits) {
 print join "\t", ($hit, @{$hits{$hit}}), "\n";
}


MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3