lib/SGN/Controller/AJAX/Heritability.pm

   1 use strict;
   2
   3 package SGN::Controller::AJAX::Heritability;
   4
   5 use Moose;
   6 use Data::Dumper;
   7 use File::Temp qw | tempfile |;
   8 use File::Slurp;
   9 use File::Spec qw | catfile|;
  10 use File::Basename qw | basename |;
  11 use File::Copy;
  12 use CXGN::Dataset;
  13 use CXGN::Dataset::File;
  14 use CXGN::Tools::Run;
  15 use CXGN::Page::UserPrefs;
  16 use CXGN::Tools::List qw/distinct evens/;
  17 use CXGN::Blast::Parse;
  18 use CXGN::Blast::SeqQuery;
  19 use Cwd qw(cwd);
  20
  21 BEGIN { extends 'Catalyst::Controller::REST' }
  22
  23 __PACKAGE__->config(
  24     default   => 'application/json',
  25     stash_key => 'rest',
  26     map       => { 'application/json' => 'JSON' },
  27     );
  28
  29
  30 sub shared_phenotypes: Path('/ajax/heritability/shared_phenotypes') : {
  31     my $self = shift;
  32     my $c = shift;
  33     my $dataset_id = $c->req->param('dataset_id');
  34
  35     my $exclude_outliers = $c->req->param('dataset_trait_outliers');
  36
  37     my $sp_person_id = $c->user() ? $c->user->get_object()->get_sp_person_id() : undef;
  38     my $people_schema = $c->dbic_schema("CXGN::People::Schema", undef, $sp_person_id);
  39     my $schema = $c->dbic_schema("Bio::Chado::Schema", "sgn_chado", $sp_person_id);
  40     my $ds = CXGN::Dataset->new(people_schema => $people_schema, schema => $schema, sp_dataset_id => $dataset_id);
  41     my $traits = $ds->retrieve_traits();
  42     my @trait_info;
  43     foreach my $t (@$traits) {
  44           my $tobj = CXGN::Cvterm->new({ schema=>$schema, cvterm_id => $t->[0] });
  45         push @trait_info, [ $tobj->cvterm_id(), $tobj->name()];
  46     }
  47
  48     $c->tempfiles_subdir("heritability_files");
  49     my ($fh, $tempfile) = $c->tempfile(TEMPLATE=>"heritability_files/trait_XXXXX");
  50
  51     my $temppath = $c->config->{basepath}."/".$tempfile;
  52     print STDERR "***** temppath = $temppath\n";
  53     my $ds2 = CXGN::Dataset::File->new(people_schema => $people_schema, schema => $schema, sp_dataset_id => $dataset_id, exclude_dataset_outliers => $exclude_outliers, file_name => $temppath, quotes => 0);
  54     my $phenotype_data_ref = $ds2->retrieve_phenotypes();
  55
  56     print STDERR Dumper(@trait_info);
  57     $c->stash->{rest} = {
  58         options => \@trait_info,
  59         tempfile => $tempfile."_phenotype.txt",
  60     };
  61 }
  62
  63
  64
  65 sub extract_trait_data :Path('/ajax/heritability/getdata') Args(0) {
  66     my $self = shift;
  67     my $c = shift;
  68
  69     my $file = $c->req->param("file");
  70     my $trait = $c->req->param("trait");
  71
  72     $file = basename($file);
  73
  74     my $temppath = File::Spec->catfile($c->config->{basepath}, "static/documents/tempfiles/heritability_files/".$file);
  75     print STDERR Dumper($temppath);
  76
  77     my $F;
  78     if (! open($F, "<", $temppath)) {
  79         $c->stash->{rest} = { error => "Can't find data." };
  80         return;
  81     }
  82
  83     my $header = <$F>;
  84     chomp($header);
  85     print STDERR Dumper($header);
  86     my @keys = split("\t", $header);
  87     print STDERR Dumper($keys[1]);
  88     for(my $n=0; $n <@keys; $n++) {
  89         if ($keys[$n] =~ /\|CO\_/) {
  90         $keys[$n] =~ s/\|CO\_.*//;
  91         }
  92     }
  93     my @data = ();
  94
  95     while (<$F>) {
  96         chomp;
  97
  98         my @fields = split "\t";
  99         my %line;
 100         for(my $n=0; $n <@keys; $n++) {
 101             if (exists($fields[$n]) && defined($fields[$n])) {
 102                 $line{$keys[$n]}=$fields[$n];
 103             }
 104         }
 105         #print STDERR Dumper(\%line);
 106         push @data, \%line;
 107     }
 108
 109     $c->stash->{rest} = { data => \@data, trait => $trait};
 110 }
 111
 112 sub generate_results: Path('/ajax/heritability/generate_results') : {
 113     my $self = shift;
 114     my $c = shift;
 115     my $dataset_id = $c->req->param('dataset_id');
 116     my $trait_id = $c->req->param('trait_id');
 117
 118     print STDERR "The dataset is $dataset_id\n";
 119     print STDERR $dataset_id;
 120     print STDERR $trait_id;
 121
 122
 123     my $exclude_outliers = $c->req->param('dataset_trait_outliers');
 124
 125     $c->tempfiles_subdir("heritability_files");
 126     my $heritability_tmp_output = $c->config->{cluster_shared_tempdir}."/heritability_files";
 127     mkdir $heritability_tmp_output if ! -d $heritability_tmp_output;
 128     print STDERR "heritability_files subdir = $heritability_tmp_output\n";
 129     my ($tmp_fh, $tempfile) = tempfile(
 130       "h2_download_XXXXX",
 131       DIR=> $heritability_tmp_output,
 132     );
 133
 134     print STDERR "TEMPFILE NOW = $tempfile\n";
 135
 136     my $pheno_filepath = $tempfile . "_phenotype.txt";
 137
 138     my $sp_person_id = $c->user() ? $c->user->get_object()->get_sp_person_id() : undef;
 139     my $people_schema = $c->dbic_schema("CXGN::People::Schema", undef, $sp_person_id);
 140     my $schema = $c->dbic_schema("Bio::Chado::Schema", "sgn_chado", $sp_person_id);
 141
 142     #my $temppath = $heritability_tmp_output . "/" . $tempfile;
 143     my $temppath = $tempfile;
 144
 145     my $ds = CXGN::Dataset::File->new(people_schema => $people_schema, schema => $schema, sp_dataset_id => $dataset_id, exclude_dataset_outliers => $exclude_outliers, file_name => $temppath, quotes => 0);
 146
 147     my $phenotype_data_ref = $ds->retrieve_phenotypes($pheno_filepath);
 148
 149
 150
 151     my $h2File = $tempfile . "_" . "h2File.json";
 152     my $h2CsvFile = $tempfile . "_" . "h2CsvFile.csv";
 153     my $errorFile = $tempfile . "_" . "error.txt";
 154
 155
 156     $trait_id =~ tr/ /./;
 157     $trait_id =~ tr/\//./;
 158
 159     my $cmd = CXGN::Tools::Run->new({
 160             backend => $c->config->{backend},
 161             submit_host=>$c->config->{cluster_host},
 162             temp_base => $c->config->{cluster_shared_tempdir} . "/heritability_files",
 163             queue => $c->config->{'web_cluster_queue'},
 164             do_cleanup => 0,
 165             # don't block and wait if the cluster looks full
 166             max_cluster_jobs => 1_000_000_000,
 167         });
 168
 169         print STDERR Dumper $pheno_filepath;
 170
 171     # my $job;
 172     $cmd->run_cluster(
 173             "Rscript ",
 174             $c->config->{basepath} . "/R/heritability/h2_blup_rscript.R",
 175             $pheno_filepath,
 176             $trait_id,
 177             $h2File,
 178             $h2CsvFile,
 179             $errorFile
 180     );
 181     $cmd->alive;
 182     $cmd->is_cluster(1);
 183     $cmd->wait;
 184
 185
 186     my $figure_path = $c->{basepath} . "./documents/tempfiles/heritability_files/";
 187     copy($h2File, $figure_path);
 188     copy($h2CsvFile, $figure_path);
 189
 190
 191     my $h2Filebasename = basename($h2File);
 192     my $h2File_response = "/documents/tempfiles/heritability_files/" . $h2Filebasename;
 193
 194     my $h2CsvFilebasename = basename($h2CsvFile);
 195     my $h2CsvFile_response = "/documents/tempfiles/heritability_files/" . $h2CsvFilebasename;
 196
 197     my $errors;
 198     if ( -e $errorFile ) {
 199         open my $fh, '<', $errorFile or die "Can't open error file $!";
 200         $errors = do { local $/; <$fh> };
 201     }
 202
 203     $c->stash->{rest} = {
 204         h2Table => $h2File_response,
 205         dummy_response => $dataset_id,
 206         error => $errors,
 207         h2CsvTable => $h2CsvFile_response
 208     };
 209 }
 210
 211 1
 212