[Bio] / DrugTargets / get_homologs.pl Repository:
ViewVC logotype

View of /DrugTargets/get_homologs.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (download) (as text) (annotate)
Thu Aug 25 23:22:54 2005 UTC (14 years, 3 months ago) by fangfang
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +6 -1 lines
Fixed the conservation script bug thanks to Bob.
Added category codes instructions page.

$usage =  "get_homologs [-cons] [-max=set_size] [-detail] FIG_ID \n";

while (($arg = shift @ARGV) =~ /^-/) {
    if ($arg =~ /-max=(\d+)/i) {
	$max = $1;
    } elsif ($arg =~ /^-detail/i) {
	$detail = 1;
    } elsif ($arg =~ /^-cons/i) {
	$for_cons = 1;
    }
}

($arg) || die "Usage: $usage";

use FIG;
use FIG_Config;

my $fig = new FIG; 
my $peg = $arg;
my $org = ($peg =~ /fig\|(\d+\.\d+)/, $1);
#print $org. "\n";

my $list_file;
my @pinned;

if ($for_cons) {
    if ($fig->taxonomy_of($org) =~ /firmicutes/i) {
	$list_file = "cons_orgs_list1.txt";
    } else {
	$list_file = "cons_orgs_list2.txt";
    }
    $list_file = "$FIG_Config::fig/var/DrugTargets/$list_file" unless (-e $list_file);
    @pinned = &get_homologs_for_cons($peg);
} else {
    @pinned = &get_homologs_by_sims($peg, $max);
}
#print STDERR &Dumper(\@pinned);

my $pegs = join(" ", @pinned);
#print STDERR $pegs . "\n";

my $output = `$FIG_Config::bin/align_with_clustal $pegs | $FIG_Config::bin/clustal_to_fasta`;
($?) && die "Alignment failed with error $?: $!\n";
print $output;

sub get_homologs_for_cons {
    my ($peg) = @_;
    my $find_homo = "$FIG_Config::bin/find_homologs_in";

    $find_homo = "perl find_homologs_in.pl" unless (-e $find_homo);
    #print STDERR "$list_file\n";
    $list  = `cut -f1 $list_file`;
    $list  =~ s/\n/ /sg;
    $list  =~ s/\#\S+//g;
    $list  =~ s/$org//g;
    #print $list;
    $homos = `$find_homo '$peg' $list`;
    #print $homos;
    #print "perl find_homologs_in.pl $peg $list";
    my @homologs = ("'$peg'");
    
    for (split(/\n/, $homos)) {
	$_ = (/^(\S+)/, "'$1'");
	push(@homologs, $_);
    }
    #print STDERR "$homos";
    return @homologs;
}

sub get_homologs_by_sims {
    my($peg, $max) = @_; 
    my($maxN,$maxP,$genome1,$sim,$id2,$genome2,%seen);
 
    #$maxN = 1000;
    $maxN = 100;
    $maxP = 1.0e-5; 
    my @sims = $fig->sims($peg, $maxN, $maxP, "fig"); 
 
    my @homologs = (); 
    $seen{&FIG::genome_of($peg)} = 1; 
    foreach $sim (@sims) 
    { 
        $id2     = $sim->id2;
	#print join("\t", $sim->iden, $sim->psc, $id2) . "\n";
        $genome2 = &FIG::genome_of($id2); 
	$taxon2  = $fig->taxonomy_of($genome2);
        my @coup; 
        #if ((! $seen{$genome2}) && (@coup = $fig->coupled_to($id2)) && (@coup > 0) && ($taxon2 =~ /Firmicutes/i)) 
	#if ((! $seen{$genome2}) && ($taxon2 =~ /Firmicutes/i)) 
	if (! $seen{$genome2}) 
        { 
            $seen{$genome2} = 1;
	    if ($detail) {
		print STDERR join("\t", $id2,  $sim->iden, $sim->psc, $fig->genus_species($genome2)) . "\n";	    
	    } else {
		print STDERR join("\t", $id2,  $fig->genus_species($genome2)) . "\n";
	    }
	    
            push(@homologs, "'$id2'");
	    last unless (!$max || (@homologs < $max));
        } 
    } 
    return sort { $b->[1] <=> $a->[1] } @homologs; 
} 

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3