lib/SGN/Controller/solGS/Correlation.pm

   1 package SGN::Controller::solGS::Correlation;
   2
   3 use Moose;
   4 use namespace::autoclean;
   5
   6 use Cache::File;
   7 use CXGN::Tools::Run;
   8 use File::Temp qw / tempfile tempdir /;
   9 use File::Spec::Functions qw / catfile catdir/;
  10 use File::Slurp qw /write_file read_file/;
  11 use File::Path qw / mkpath  /;
  12 use File::Copy;
  13 use File::Basename;
  14 use CXGN::Phenome::Population;
  15 use JSON;
  16 use Try::Tiny;
  17
  18 BEGIN { extends 'Catalyst::Controller' }
  19
  20
  21 sub check_pheno_corr_result :Path('/phenotype/correlation/check/result/') Args(1) {
  22     my ($self, $c, $pop_id) = @_;
  23
  24     $c->stash->{pop_id} = $pop_id;
  25
  26     $self->pheno_correlation_output_files($c);
  27     my $corre_output_file = $c->stash->{corre_coefficients_json_file};
  28
  29     my $ret->{result} ='No';
  30
  31     if (-s $corre_output_file && $pop_id =~ /\d+/)
  32     {
  33         $ret->{result} = 'yes';
  34     }
  35
  36     $ret = to_json($ret);
  37
  38     $c->res->content_type('application/json');
  39     $c->res->body($ret);
  40
  41 }
  42
  43
  44 sub correlation_phenotype_data :Path('/correlation/phenotype/data/') Args(0) {
  45     my ($self, $c) = @_;
  46
  47     my $pop_id = $c->req->param('population_id');
  48     $c->stash->{pop_id} = $pop_id;
  49     my $referer = $c->req->referer;
  50
  51     my $phenotype_file;
  52
  53     if( $pop_id =~ /uploaded/)
  54     {
  55         my $phenotype_dir = $c->stash->{solgs_prediction_upload_dir};
  56         my $userid        = $c->user->id;
  57         $phenotype_file   = "phenotype_data_${userid}_${pop_id}";
  58         $phenotype_file   = $c->controller('solGS::solGS')->grep_file($phenotype_dir, $phenotype_file);
  59     }
  60     elsif ($referer =~ /qtl/)
  61     {
  62         $self->create_correlation_phenodata_file($c);
  63         $phenotype_file =  $c->stash->{phenotype_file};
  64     }
  65     else
  66     {
  67         my $phenotype_dir = $c->stash->{solgs_cache_dir};
  68         $phenotype_file   = 'phenotype_data_' . $pop_id;
  69         $phenotype_file   = $c->controller('solGS::solGS')->grep_file($phenotype_dir, '\'^' . $phenotype_file . '\'');
  70     }
  71
  72     unless ($phenotype_file)
  73     {
  74         $self->create_correlation_phenodata_file($c);
  75         $phenotype_file =  $c->stash->{phenotype_file};
  76     }
  77
  78     my $ret->{status} = 'failed';
  79
  80     if (-s $phenotype_file)
  81     {
  82         $ret->{status} = 'success';
  83     }
  84     else
  85     {
  86         $ret->{status} = 'This population set has no phenotype data.';
  87     }
  88
  89     $ret = to_json($ret);
  90
  91     $c->res->content_type('application/json');
  92     $c->res->body($ret);
  93
  94 }
  95
  96
  97 sub correlation_genetic_data :Path('/correlation/genetic/data/') Args(0) {
  98     my ($self, $c) = @_;
  99
 100     my $corr_pop_id = $c->req->param('corr_population_id');
 101     my $pop_type    = $c->req->param('type');
 102     my $model_id    = $c->req->param('model_id');
 103
 104     my $index_file  = $c->req->param('index_file');
 105
 106     $c->stash->{model_id} = $model_id;
 107     $c->stash->{pop_id}   = $model_id;
 108
 109     $c->stash->{prediction_pop_id} = $corr_pop_id if $pop_type =~ /selection/;
 110
 111     $c->stash->{selection_index_file} = $index_file;
 112     $self->combine_gebvs_of_traits($c);
 113     my $combined_gebvs_file = $c->stash->{combined_gebvs_file};
 114
 115     my $ret->{status} = 'failed';
 116
 117     if ( -s $combined_gebvs_file )
 118     {
 119         $ret->{status} = 'success';
 120         $ret->{gebvs_file} = $combined_gebvs_file;
 121     }
 122
 123     $ret = to_json($ret);
 124
 125     $c->res->content_type('application/json');
 126     $c->res->body($ret);
 127
 128 }
 129
 130
 131 sub combine_gebvs_of_traits {
 132     my ($self, $c) = @_;
 133
 134     $c->controller("solGS::solGS")->get_gebv_files_of_traits($c);
 135     my $gebvs_files = $c->stash->{gebv_files_of_valid_traits};
 136
 137     if (!-s $gebvs_files)
 138     {
 139         $gebvs_files = $c->stash->{gebv_files_of_traits};
 140     }
 141
 142     my $index_file  = $c->stash->{selection_index_file};
 143
 144     my @files_no = map { split(/\t/) } read_file($gebvs_files);
 145
 146     if (scalar(@files_no) > 1 )
 147     {
 148
 149         if ($index_file)
 150         {
 151             write_file($gebvs_files, {append => 1}, "\t". $index_file )
 152         }
 153
 154         my $pred_pop_id = $c->stash->{prediction_pop_id};
 155         my $model_id    = $c->stash->{model_id};
 156         my $identifier  =  $pred_pop_id ? $model_id . "_" . $pred_pop_id :  $model_id;
 157
 158         my $combined_gebvs_file = $c->controller("solGS::solGS")->create_tempfile($c, "combined_gebvs_${identifier}");
 159
 160         $c->stash->{input_files}  = $gebvs_files;
 161         $c->stash->{output_files} = $combined_gebvs_file;
 162         $c->stash->{r_temp_file}  = "combining-gebvs-${identifier}";
 163         $c->stash->{r_script}     = 'R/combine_gebvs_files.r';
 164
 165         $c->controller("solGS::solGS")->run_r_script($c);
 166         $c->stash->{combined_gebvs_file} = $combined_gebvs_file;
 167     }
 168     else
 169     {
 170         $c->stash->{combined_gebvs_files} = 0;
 171     }
 172 }
 173
 174
 175 sub create_correlation_phenodata_file {
 176     my ($self, $c)  = @_;
 177     my $referer = $c->req->referer;
 178
 179     if ($referer =~ /qtl/)
 180     {
 181         my $pop_id = $c->stash->{pop_id};
 182
 183         my $pheno_exp = "phenodata_${pop_id}";
 184         my $dir       = catdir($c->config->{solqtl}, 'cache');
 185
 186         my $phenotype_file = $c->controller("solGS::solGS")->grep_file($dir, $pheno_exp);
 187
 188         unless ($phenotype_file)
 189         {
 190             my $pop =  CXGN::Phenome::Population->new($c->dbc->dbh, $pop_id);
 191             $phenotype_file =  $pop->phenotype_file($c);
 192         }
 193
 194         my $new_file = catfile($c->stash->{correlation_dir}, "phenotype_data_${pop_id}.csv");
 195
 196         copy($phenotype_file, $new_file)
 197             or die "could not copy $phenotype_file to $new_file";
 198
 199         $c->stash->{phenotype_file} = $new_file;
 200     }
 201     else
 202     {
 203       $c->controller("solGS::solGS")->phenotype_file($c);
 204     }
 205
 206 }
 207
 208
 209 sub create_correlation_dir {
 210     my ($self, $c) = @_;
 211
 212     $c->controller("solGS::solGS")->get_solgs_dirs($c);
 213
 214 }
 215
 216
 217 sub pheno_correlation_output_files {
 218     my ($self, $c) = @_;
 219
 220     my $pop_id = $c->stash->{pop_id};
 221
 222     $self->create_correlation_dir($c);
 223     my $corre_dir = $c->stash->{correlation_dir};
 224
 225     my $file_cache  = Cache::File->new(cache_root => $corre_dir);
 226     $file_cache->purge();
 227
 228     my $key_table = 'corre_coefficients_table_' . $pop_id;
 229     my $key_json  = 'corre_coefficients_json_' . $pop_id;
 230     my $corre_coefficients_file      = $file_cache->get($key_table);
 231     my $corre_coefficients_json_file = $file_cache->get($key_json);
 232
 233     unless ($corre_coefficients_file && $corre_coefficients_json_file )
 234     {
 235         $corre_coefficients_file= catfile($corre_dir, "corre_coefficients_table_${pop_id}");
 236
 237         write_file($corre_coefficients_file);
 238         $file_cache->set($key_table, $corre_coefficients_file, '30 days');
 239
 240         $corre_coefficients_json_file = catfile($corre_dir, "corre_coefficients_json_${pop_id}");
 241
 242         write_file($corre_coefficients_json_file);
 243         $file_cache->set($key_json, $corre_coefficients_json_file, '30 days');
 244     }
 245
 246     $c->stash->{corre_coefficients_table_file} = $corre_coefficients_file;
 247     $c->stash->{corre_coefficients_json_file}  = $corre_coefficients_json_file;
 248 }
 249
 250
 251 sub genetic_correlation_output_files {
 252     my ($self, $c) = @_;
 253
 254     my $corre_pop_id = $c->stash->{corre_pop_id};
 255     my $model_id     = $c->stash->{model_id};
 256     my $type         = $c->stash->{type};
 257
 258     my $pred_pop_id = $c->stash->{prediction_pop_id};
 259     $model_id    = $c->stash->{model_id};
 260     my $identifier  =  $type =~ /selection/ ? $model_id . "_" . $corre_pop_id :  $corre_pop_id;
 261
 262     my $solgs_controller = $c->controller("solGS::solGS");
 263     my $corre_json_file  = $solgs_controller->create_tempfile($c, "genetic_corre_json_${identifier}");
 264     my $corre_table_file = $solgs_controller->create_tempfile($c, "genetic_corre_table_${identifier}");
 265
 266     $c->stash->{genetic_corre_table_file} = $corre_table_file;
 267     $c->stash->{genetic_corre_json_file}  = $corre_json_file;
 268 }
 269
 270
 271 sub pheno_correlation_analysis_output :Path('/phenotypic/correlation/analysis/output') Args(0) {
 272     my ($self, $c) = @_;
 273
 274     my $pop_id = $c->req->param('population_id');
 275     $c->stash->{pop_id} = $pop_id;
 276
 277     $self->pheno_correlation_output_files($c);
 278     my $corre_json_file = $c->stash->{corre_coefficients_json_file};
 279
 280     my $ret->{status} = 'failed';
 281
 282     if (!-s $corre_json_file)
 283     {
 284         $self->run_pheno_correlation_analysis($c);
 285         $corre_json_file = $c->stash->{corre_coefficients_json_file};
 286     }
 287
 288     if (-s $corre_json_file)
 289     {
 290         $ret->{status}   = 'success';
 291         $ret->{data}     = read_file($corre_json_file);
 292     }
 293
 294     $ret = to_json($ret);
 295
 296     $c->res->content_type('application/json');
 297     $c->res->body($ret);
 298
 299 }
 300
 301
 302 sub genetic_correlation_analysis_output :Path('/genetic/correlation/analysis/output') Args(0) {
 303     my ($self, $c) = @_;
 304
 305     $c->stash->{corre_pop_id} = $c->req->param('corr_population_id');
 306     $c->stash->{model_id}     = $c->req->param('model_id');
 307     $c->stash->{type}         = $c->req->param('type');
 308
 309     my $corr_pop_id = $c->req->param('corr_population_id');
 310     my $model_id    = $c->req->param('model_id');
 311     my $type        = $c->req->param('type');
 312
 313     my $gebvs_file = $c->req->param('gebvs_file');
 314     $c->stash->{data_input_file} = $gebvs_file;
 315
 316     $self->genetic_correlation_output_files($c);
 317
 318     if (-s $gebvs_file)
 319     {
 320         $self->run_genetic_correlation_analysis($c);
 321     }
 322
 323     my $ret->{status} = 'failed';
 324     my $corre_json_file = $c->stash->{genetic_corre_json_file};
 325
 326     if (-s $corre_json_file)
 327     {
 328         $ret->{status}   = 'success';
 329         $ret->{data}     = read_file($corre_json_file);
 330     }
 331
 332     $ret = to_json($ret);
 333
 334     $c->res->content_type('application/json');
 335     $c->res->body($ret);
 336
 337 }
 338
 339
 340 sub run_pheno_correlation_analysis {
 341     my ($self, $c) = @_;
 342
 343     my $pop_id = $c->stash->{pop_id};
 344
 345     $self->create_correlation_phenodata_file($c);
 346     $c->stash->{data_input_file} = $c->stash->{phenotype_file};
 347
 348     $self->pheno_correlation_output_files($c);
 349     $c->stash->{corre_table_output_file} = $c->stash->{corre_coefficients_table_file};
 350     $c->stash->{corre_json_output_file}  = $c->stash->{corre_coefficients_json_file};
 351
 352     $c->controller("solGS::solGS")->formatted_phenotype_file($c);
 353
 354     $c->stash->{referer} = $c->req->referer;
 355
 356     $c->stash->{correlation_type} = "pheno_correlation_${pop_id}";
 357     $c->stash->{correlation_script} = "R/phenotypic_correlation.r";
 358
 359     $self->run_correlation_analysis($c);
 360
 361 }
 362
 363
 364 sub run_genetic_correlation_analysis {
 365     my ($self, $c) = @_;
 366
 367     my $pop_id = $c->stash->{corre_pop_id};
 368
 369     $self->genetic_correlation_output_files($c);
 370     $c->stash->{corre_table_output_file} = $c->stash->{genetic_corre_table_file};
 371     $c->stash->{corre_json_output_file}  = $c->stash->{genetic_corre_json_file};
 372
 373     $c->stash->{referer} = $c->req->referer;
 374
 375     $c->stash->{correlation_type} = "genetic_correlation_${pop_id}";
 376     $c->stash->{correlation_script} = "R/genetic_correlation.r";
 377     $self->run_correlation_analysis($c);
 378
 379 }
 380
 381
 382 sub download_phenotypic_correlation : Path('/download/phenotypic/correlation/population') Args(1) {
 383     my ($self, $c, $id) = @_;
 384
 385     $self->create_correlation_dir($c);
 386     my $corr_dir = $c->stash->{correlation_dir};
 387     my $corr_file = catfile($corr_dir,  "corre_coefficients_table_${id}");
 388
 389     unless (!-e $corr_file || -s $corr_file <= 1)
 390     {
 391         my @corr_data;
 392         my $count=1;
 393
 394         foreach my $row ( read_file($corr_file) )
 395         {
 396             if ($count==1) {  $row = 'Traits,' . $row;}
 397             $row =~ s/NA//g;
 398             $row = join(",", split(/\s/, $row));
 399             $row .= "\n";
 400
 401             push @corr_data, [ $row ];
 402             $count++;
 403         }
 404
 405         $c->res->content_type("text/plain");
 406         $c->res->body(join "",  map{ $_->[0] } @corr_data);
 407
 408
 409     }
 410 }
 411
 412
 413 sub run_correlation_analysis {
 414     my ($self, $c) = @_;
 415
 416     my $pop_id = $c->stash->{pop_id};
 417
 418     $self->create_correlation_dir($c);
 419     my $corre_dir = $c->stash->{correlation_dir};
 420
 421     my $data_input_file = $c->stash->{data_input_file};
 422
 423     my $corre_table_file = $c->stash->{corre_table_output_file};
 424     my $corre_json_file  = $c->stash->{corre_json_output_file};
 425
 426     my $formatted_phenotype_file = $c->stash->{formatted_phenotype_file};
 427
 428     my $referer        = $c->stash->{referer};
 429     my $corre_analysis = $c->stash->{correlation_type};
 430     my $corre_script   = $c->stash->{correlation_script};
 431
 432     if (-s $data_input_file)
 433     {
 434         CXGN::Tools::Run->temp_base($corre_dir);
 435
 436         my ( $corre_commands_temp, $corre_output_temp ) =
 437             map
 438         {
 439             my (undef, $filename ) =
 440                 tempfile(
 441                     catfile(
 442                         CXGN::Tools::Run->temp_base(),
 443                         "$corre_analysis-$_-XXXXXX",
 444                          ),
 445                 );
 446             $filename
 447         } qw / in out /;
 448
 449     {
 450         my $corre_commands_file = $c->path_to($corre_script);
 451         copy( $corre_commands_file, $corre_commands_temp )
 452             or die "could not copy '$corre_commands_file' to '$corre_commands_temp'";
 453     }
 454
 455       try
 456       {
 457           print STDERR "\nsubmitting correlation job to the cluster..\n";
 458           my $r_process = CXGN::Tools::Run->run_cluster(
 459               'R', 'CMD', 'BATCH',
 460               '--slave',
 461               "--args $formatted_phenotype_file $referer $corre_table_file $corre_json_file $data_input_file",
 462               $corre_commands_temp,
 463               $corre_output_temp,
 464               {
 465                   working_dir => $corre_dir,
 466                   max_cluster_jobs => 1_000_000_000,
 467               },
 468               );
 469
 470           $r_process->wait;
 471           print STDERR "\ndone with correlation analysis..\n";
 472       }
 473       catch
 474       {
 475
 476             my $err = $_;
 477             $err =~ s/\n at .+//s; #< remove any additional backtrace
 478             #     # try to append the R output
 479
 480             try
 481             {
 482                 $err .= "\n=== R output ===\n".file($corre_output_temp)->slurp."\n=== end R output ===\n"
 483             };
 484
 485             $c->stash->{script_error} = "Correlation analysis failed.";
 486
 487       };
 488     }
 489 }
 490
 491
 492 sub begin : Private {
 493     my ($self, $c) = @_;
 494
 495     $c->controller("solGS::solGS")->get_solgs_dirs($c);
 496
 497 }
 498
 499
 500
 501 ####
 502 1;
 503 ####