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

View of /FigKernelScripts/ma_pipeline.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.15 - (download) (as text) (annotate)
Wed Sep 15 15:54:52 2010 UTC (9 years, 5 months ago) by overbeek
Branch: MAIN
CVS Tags: mgrast_dev_08112011, mgrast_dev_08022011, rast_rel_2014_0912, myrast_rel40, mgrast_dev_05262011, mgrast_dev_04082011, rast_rel_2010_0928, mgrast_version_3_2, mgrast_dev_12152011, mgrast_dev_06072011, rast_rel_2014_0729, mgrast_dev_02212011, rast_rel_2010_1206, mgrast_release_3_0, mgrast_dev_03252011, rast_rel_2011_0119, 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, myrast_33, rast_rel_2011_0928, mgrast_dev_04052011, mgrast_dev_02222011, mgrast_dev_10262011, HEAD
Changes since 1.14: +4 -4 lines
get rid of RMA Key

use strict;
use Data::Dumper;

my $usage = "usage: ma_pipeline Genome";

my($genome,$orgD);
(
 ($genome = shift @ARGV) && ($orgD = "$FIG_Config::organisms/$genome")
 )
    || die $usage;

(-d $orgD) || die "invalid organism directory";

# (my $rma_key = shift @ARGV) 
#    || die $usage;

my $exp_ed = "$orgD/UserSpace/ExpressionData";
(-d $exp_ed)
    || die "Where is $exp_ed?";
(-s "$exp_ed/probes")
    || die "Where is $exp_ed/probes?";
(-s "$orgD/contigs")
    || die "Where is $orgD/contigs?";
(-s "$orgD/Features/peg/tbl")
    || die "Where is $orgD/Features/peg/tbl?";

my $tbl = (-s "$exp_ed/tbl") ? "$exp_ed/tbl" : "$orgD/Features/peg/tbl";
&run("make_probes_to_genes $exp_ed/probes $orgD/contigs $tbl $exp_ed/peg.probe.table $exp_ed/probe.occ.table 2> $exp_ed/problems");
&run("remove_multiple_occurring_probes $exp_ed/peg.probe.table $exp_ed/probe.occ.table > $exp_ed/peg.probe.table.no.multiple");
chdir $exp_ed;

my $peg_probe_table = "$exp_ed/peg.probe.table.no.multiple";
my $probe_no_match = "$exp_ed/probe.no.match";

&make_missing_probes($peg_probe_table, "$exp_ed/probes", $probe_no_match);

&run("Rscript $FIG_Config::bin/RunRMA.R $peg_probe_table $probe_no_match $exp_ed/Experiments $exp_ed 2> $exp_ed/stderr.RMA");

# &run("Rscript $FIG_Config::bin/RunRMA.R $rma_key $exp_ed $exp_ed/Experiments 2> $exp_ed/stderr.RMA");

&run("$FIG_Config::bin/call_coregulated_clusters_on_chromosome $genome $exp_ed/raw_data.tab > $exp_ed/coregulated.clusters 2> $exp_ed/stderr.coregulated.clusters");
&run("$FIG_Config::bin/make_coreg_conjectures_based_on_subsys  $genome $exp_ed/raw_data.tab > $exp_ed/coregulated.subsys   2> $exp_ed/stderr.coregulated.subsys");
&run("cat $exp_ed/coregulated.clusters $exp_ed/coregulated.subsys | cut -f1 | $FIG_Config::bin/merge_gene_sets | filter_on_known $exp_ed/raw_data.tab > $exp_ed/merged.clusters");

&run("get_ON_probes $genome $exp_ed/probe.occ.table $exp_ed/peg.probe.table $exp_ed/raw_data.tab > $exp_ed/probes.always.on");
&run("cut -f2 $exp_ed/probes.always.on | sort -u > $exp_ed/pegs.always.on");
&run("Rscript $FIG_Config::bin/Pipeline.R $exp_ed/pegs.always.on $exp_ed/merged.clusters > $exp_ed/comments.by.Pipeline.R");

&run("Rscript $FIG_Config::bin/SplitGeneSets.R $exp_ed/merged.clusters 0.7 > $exp_ed/split.clusters");

sub make_missing_probes
{
    my($probe_table, $probes, $output) = @_;
    open(MATCH,"<", $probe_table) or die "Cannot open $probe_table: $!";
    open(PROBES,"<", $probes) or die "Cannot open $probes: $!";
    open(OUTPUT, ">", $output) or die "Cannot open $output: $!";
    my %locations;
    while(<MATCH>)
    {
	chomp;
	my($peg,$loc)=split "\t";
	$locations{$loc} = $peg;
    }
    
    while(<PROBES>)
    {
	chomp;
	my($loc,$seq) = split "\t";
	print OUTPUT $loc, "\n" if ! exists $locations{$loc};
    }
    close(MATCH);
    close(PROBES);
    close(OUTPUT);
}

sub run {
    shift if UNIVERSAL::isa($_[0],__PACKAGE__);
    my($cmd) = @_;

    if ($ENV{FIG_VERBOSE}) {
        my @tmp = `date`;
        chomp @tmp;
        print STDERR "$tmp[0]: running $cmd\n";
    }
    (system($cmd) == 0) || die("FAILED: $cmd");
}


MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3