FANDOM


#!/usr/bin/perl
use strict;
use Parallel::ChildManager;
use Digest::MD5;
use Getopt::Long;

my $scratch = "/tmp/pancoregenome.$$";
my $FORMATDB = '/usr/biotools/blast/bin/formatdb';
my $TIGRCUT = "/usr/biotools/indirect/tigrcut"; # puts 50% ALR and 50% cutoff identity on an XML output of blastall
my $BLASTALL = "/usr/biotools/blast/bin/blastall";
my $cache_source = "blastp-core+matrix";


my ($keep,$cpu,$clean,$cex);
&GetOptions (
# main options
 "keep:s"   , \$keep,       # keep temporary files in this dir ( TRY AND CREATED UNLESS EXIST )
 "cpu:s"   ,  \$cpu,        # number of CPUs - 5 is default
 "clean"  =>  \$clean ,      # wipe all existing results for this run
 "cex:s"     =>  \$cex  ,     # label scaling
);

$cex = 0.7 unless defined $cex;
$cpu = 1 unless defined $cpu;

mkdir $scratch unless -d $scratch;

# create config reading from stdin
my @config = parse_config ();

@config = prepare_fasta(@config);

# make blast databases
&make_blastdbs(@config);

# make blast reports
my $cm = new ChildManager($cpu);
make_blastreports(@config);

&group(@config);

@config = core_genome ( @config ) ;
@config = pan_genome ( @config ) ;
@config = new_genes ( @config ) ;
@config = new_families ( @config ) ;

make_table ( @config ) ;

plot(@config);

if ( defined ( $keep ) ) {
    system "mv $keep $keep.old_from_$$" if -e "$keep.old_from_$$";
    system "mv $scratch $keep";
    print STDERR "# kept temporary files in $keep\n";
}

cleanup ();

sub cleanup {
    system "rm -rf $scratch";
}

sub make_blastreports {
    my @config = @_;
    my %tab_h;
    foreach my $a ( 0 .. $#config ) {
        foreach my $b ( 0 .. $#config ) {
            $tab_h{"$a-$b"} = $cm->start();
            unless ( $tab_h{"$a-$b"} ) {
                my $blastreport_scratch = "$scratch/$a-$b.blastout.gz";
                print STDERR "# make blast report $blastreport_scratch\n";
                my $jobid = md5 ( "$scratch/$a.fsa" , "$scratch/$b.fsa" ) ;
                system "/usr/biotools/indirect/cacher --id='$jobid' --source='$cache_source' -action get > $blastreport_scratch";
                if ( $? != 0 or $clean ) {
                    print STDERR "# jobid $jobid not in cache - redoing\n";
                    system "$BLASTALL -F 0 -i $scratch/$a.fsa -p blastp -e 1e-5 -m 7 -d $scratch/$b.fsa | perl $TIGRCUT | gawk '{print \$1\"\\t\"\$2}' | gzip > $blastreport_scratch";
                    system "/usr/biotools/indirect/cacher --id='$jobid' --source='$cache_source' -action put -expire 100 < $blastreport_scratch";
                } else {
                    print STDERR "# fetched jobid $jobid from cache\n";
                }
                my $waiting = ($#config+1)**2;
                foreach my $a ( 0 .. $#config ) {
                    foreach my $b ( 0 .. $#config ) {
                        $waiting-- if -e "$scratch/$a-$b.blastout.gz";
                    }
                }
                print STDERR "# jobid $jobid finished - $waiting left\n";
                exit;
            }
        }
    }
    $cm->wait_all_children;
}

sub make_blastdbs {
    my @config = @_;
    foreach my $id ( 0 .. $#config ) {
        my $target = $config[$id]->{target};
        print STDERR "# building blast database '$target.psq'\n";
        system "$FORMATDB -i $target -p T -t $id";
    }
}

sub prepare_fasta {
    my @config = @_;
    foreach my $id ( 0 .. $#config ) {
        my ($nprot,$skipped,$checksum) = fasta2fasta( $config[$id]->{source} , "$scratch/$id.fsa");
        warn "# prepare_fasta(): nprot=$nprot, skipped=$skipped\n";
        $config[$id]->{target} = "$scratch/$id.fsa";
        $config[$id]->{'total genes'} = $nprot;
        $config[$id]->{skipped} = $skipped;
        $config[$id]->{checksum} = $checksum;
        die "nprot == 0 for $config[$id]->{source}\n" if $nprot == 0;
    }
    return @config;
}

sub fasta2fasta {
    my ($input,$output) = @_;
    my $temp = "$output.temp.raw";
    my $nprot = 0;
    my $skipped = 0;
    open RAW , "| saco_convert -I fasta -O raw | sort > $temp" or die $!;
    open FASTA , $input;
    while (<FASTA>) {
        print RAW;
    }   
    close FASTA;
    close RAW;
    my $checksum = md5 ( $temp ) ;
    open FASTA , "| saco_convert -I tab -O fasta > $output";
    open TAB , $temp or die $!;
    while (<TAB>) {
        if ( ! /^([A-Za-z]+)/ ) {
            $skipped++;
            next;
        }
        $nprot++;
        print FASTA "$checksum.$nprot\t$1\n";
    }
    close FASTA;
    close TAB;
    unlink $temp;
    return ($nprot,$skipped,$checksum);
}

sub md5 {
    my @files = @_;
    my $md5 = Digest::MD5->new;
    foreach my $file (@files) {
        open(FILE, $file) or die "Can't open '$file': $!";
        binmode(FILE);
        while (<FILE>) {
            $md5->add($_);
        }
        close(FILE);
    }
    return $md5->hexdigest;
}

sub parse_config {
    my @config;

    while ( <> ) {
        next if /^#/;
        chomp;
        my ($description,$source) = split /\t/;
        warn "# description=$description, source=$source\n";
        my $id = $#config + 1;
        ( $config[$id]->{description} , $config[$id]->{source} , $config[$id]->{id} ) = ( $description , $source , $id );
    }

    return @config;
}

sub group {
    my @config = @_;
    foreach my $id (0 .. $#config) {
        open GRP , "| perl /usr/biotools/indirect/group > $scratch/group_$id.dat";
        foreach my $x (0 .. $id) {
            foreach my $y (0 .. $id) {
                foreach my $z ($x , $y) {
                    foreach my $n ( 1 .. $config[$z]->{'total genes'} ) {
                        print GRP $config[$z]->{checksum}.".$n\n";
                    }
                }
                open BLASTOUT , "gunzip -c $scratch/$x-$y.blastout.gz |";
                while (<BLASTOUT>) {
                    print GRP $_;
                }
                close BLASTOUT;
            }
        }
        close GRP;
    }
}

sub pan_genome {
    my @config = @_;
    foreach my $id (0 .. $#config) {
        $config[$id]->{'pan genome'} = 0;
        open GRP , "$scratch/group_$id.dat";
        while (<GRP>) {
            $config[$id]->{'pan genome'}++;
        }
        close GRP;
    }
    return @config;
}

sub core_genome {
    my @config = @_;
    foreach my $id (0 .. $#config) {
        $config[$id]->{'core genome'} = 0;
        open GRP , "$scratch/group_$id.dat";
        while (<GRP>) {
            my ($group_id,$count,@members) = split /\t/;
            my %repr;
            foreach (@members) {
                next unless /^(.*)\.(.*)/;
                $repr{$1} = 1;
            }
            $config[$id]->{'core genome'}++ if scalar ( keys %repr ) == $id + 1;
        }
        close GRP;
    }
    return @config;
}

sub new_genes {
    my @config = @_;
foreach my $a ( 0 .. $#config ) {
        my %HITS;
        $config[$a]->{'new genes'} = 0;
        print STDERR "# parsing blast reports (id $a)\n";
        foreach my $b ( 0 .. $a ) {
            open GUNZIP , "gunzip -c $scratch/$a-$b.blastout.gz $scratch/$b-$a.blastout.gz |";
            while (<GUNZIP>) {
                chomp;
                my ($q,$s) = split /\t/;
                $HITS{$q}{$s} = 1;
            }
            close GUNZIP;
        }
        open FSA , "$scratch/$a.fsa" or die $!;
        while (<FSA>) {
            chomp;
            next unless /^>(.*)/;
            my $q = $1;
            my $hits_in_other_samples;
            foreach (keys %{$HITS{$q}}) {
                next unless /^(.*)\./;
                if ( $config[$a]->{checksum} ne $1) {
                    $hits_in_other_samples = 1;
                    last;
                }
            }
            $config[$a]->{'new genes'}++ if ! defined $hits_in_other_samples;
        }
    }
    return @config;
}

sub new_families {
    my @config = @_;
    foreach my $id (0 .. $#config) {
        $config[$id]->{'new families'} = 0;
        open GRP , "$scratch/group_$id.dat";
        while (<GRP>) {
            my ($group_id,$count,@members) = split /\t/;
            my $is_new_family = 1;
            foreach (@members) {
                next unless /^(.*)\.(.*)/;
                $is_new_family = 0 if $1 ne $config[$id]->{checksum};
            }
            $config[$id]->{'new families'} += $is_new_family;
        }
        close GRP;
    }
    return @config;
}

sub make_table {
    my @config = @_;
    foreach my $id (0 .. $#config) {
        print $config[$id]->{'new families'},"\t";
        print $config[$id]->{'new genes'},"\t";
        print $config[$id]->{'core genome'},"\t";
        print $config[$id]->{'pan genome'},"\n";
    }
}

sub plot {
    my @config = @_;
    my $tbl = "$scratch/tbl";
    open TBL , ">$tbl";
    my @keys = ("id","description","total genes","new genes","new families","pan genome","core genome");
    print TBL join ("\t",@keys)."\n";
    foreach my $id ( 0 .. $#config ) {
        foreach my $key (@keys) {
            print TBL $config[$id]->{$key},"\t";
            }
    print TBL "\n";
    }
    close TBL;
    open R , "| /usr/bin/R --vanilla > /dev/null";
    print R "
postscript('$scratch/ps');
layout( matrix(c(1,2), 1, 2, byrow = TRUE),widths=c(0.5,0.5))
data <- read.table('$tbl',skip=1,sep='\t',dec='.',header=FALSE);
tN <- data[,3]
x<-rbind(data[,4],data[,5])
r <- barplot(x,beside=TRUE,cex.names=$cex,names.arg=as.vector(data[,1])+1,ylim=c(0,1.3*max(data[,6])))
lines(r[2,], type='b', data[,7], col='red', lwd=4)
lines(r[2,], type='b', data[,6], col='blue', lwd=4)
legend(1,1.3*max(data[,6]),c('New genes','New gene families','Core genome','Pan genome'),col=c('darkgray','lightgray','red','blue'), lwd=c(4,4,4,4))
plot.new()
plot.window(xlim=c(0,2),ylim=c(-nrow(data),1))
for (i in 1:nrow(data)) {
 text(0,-i*0.6, adj=0,paste(i,': ',data[i,2]),cex=$cex)
}
dev.off();
";
    close R;
    system "cat $scratch/ps";
}

__END__

NAME
    coregenome - derive core-/pan genome as well as count new genes/families in
    a list of genomes/samples


SYNOPSIS
    perl pancoreplot [-cpu N] [-clean] list

DESCRIPTION
    coregenome reads a list from STDIN, defining the proteomes to be compared.
    The order of the list indicate in what order the genomes will appear in the
    plot. Each line must contain a description and data source, separated by
    tab. Example:

    Campylobacter jejuni subsp. jejuni NCTC 11168    saco_extract -t -I genbank -O fasta /home/projects/pfh/genomes/data/AL111168/AL111168.gbk |
    Campylobacter jejuni RM1221    saco_extract -t -I genbank -O fasta /home/projects/pfh/genomes/data/CP000025/CP000025.gbk |
    Campylobacter fetus subsp. fetus 82-40    saco_extract -t -I genbank -O fasta /home/projects/pfh/genomes/data/CP000487/CP000487.gbk |
    Campylobacter jejuni subsp. jejuni 81-176    saco_extract -t -I genbank -O fasta /home/projects/pfh/genomes/data/CP000538/CP000538.gbk |
    Campylobacter jejuni subsp. doylei 269.97    saco_extract -t -I genbank -O fasta /home/projects/pfh/genomes/data/CP000768/CP000768.gbk |
    Campylobacter curvus 525.92    saco_extract -t -I genbank -O fasta /home/projects/pfh/genomes/data/CP000767/CP000767.gbk |
    Campylobacter hominis ATCC BAA-381    saco_extract -t -I genbank -O fasta /home/projects/pfh/genomes/data/CP000776/CP000776.gbk |
    Campylobacter concisus 13826    saco_extract -t -I genbank -O fasta /home/projects/pfh/genomes/data/CP000792/CP000792.gbk |
    Campylobacter jejuni subsp. jejuni 81116    saco_extract -t -I genbank -O fasta /home/projects/pfh/genomes/data/CP000814/CP000814.gbk |

    Note that the above exampe does not specify files as the data source, rather a command to obtain the data (tailing pipe is important)

    Each proteome is blasted against all previous proteoms in the list. All
    blast results are kept in the directory './cache', and the naming is
    derived based in MD5 checksums of the raw proteins sequences. Example:

    GenomeA    A.fsa
    GenomeB    B.fsa
    GenomeC    C.fsa

        Providing the list as above will cause all proteoms to be BLAST againast all.

    The program will for each proteome X calculate the follwing:

    - Number of proteins observed in X
    - Number of new genes observed in genome in X compared to 0 .. (X-1)
    - Number of new gene families observed in genome in X compared to 0 .. (X-1)
    - Size of pan genome at genome X
    - Size of core genome at genome X

    The program assumes two proteins being similar if the aligment identified by
    BLASTP (evalue <= 1e-5) spans 50% or more of the longest of the sequence AND
    that 50% of the residues are conserved within the alignment. Only the reciprocal match
    is considered.

OPTIONS
    -clean
       Ignores existing and cached results and regenerates the blast reports

    -cpu <N>
       Specify how many BLAST searches run in parallel.

    -table <filename>
       Specify the filename to which the result table will be written
 
VERSION
    This version (2.2, November 27, 2008) replaces all older versions of coregenome as older version had known bugs.

EXAMPLE
    mysql -B -N -e "select genbank,concat('saco_extract -t -I genbank -O fasta \
      /home/projects/pfh/genomes/data/',genbank,'/',genbank,'.gbk |') from pfh_public.genbank_complete_seq as s , \
      pfh_public.genbank_complete_prj as p where s.pid = p.pid and p.organism_name like 'campy%' \
      and genbank not like 'genome%' and segment_name like 'chromo%' order by released" | \
      perl coregenome-2.2 -cpu 10 > output.ps


AUTHOR
    Peter Fischer Hallin, pfh@cbs.dtu.dk

Ad blocker interference detected!


Wikia is a free-to-use site that makes money from advertising. We have a modified experience for viewers using ad blockers

Wikia is not accessible if you’ve made further modifications. Remove the custom ad blocker rule(s) and the page will load as expected.