cgi-bin/phenome/qtl_analysis.pl

   1 use strict;
   2 use warnings;
   3
   4 my $qtl_analysis_detail_page =
   5   CXGN::Phenome::QtlAnalysisDetailPage->new();
   6
   7 package CXGN::Phenome::QtlAnalysisDetailPage;
   8
   9
  10
  11 use CXGN::Page;
  12 use CXGN::Page::FormattingHelpers qw /info_section_html
  13   page_title_html
  14   columnar_table_html
  15   html_optional_show
  16   info_table_html
  17   tooltipped_text
  18   html_alternate_show
  19   /;
  20
  21 use CXGN::Phenome::Population;
  22 use CXGN::Phenome::Qtl;
  23 use CXGN::Phenome::PopulationDbxref;
  24 use CXGN::Tools::WebImageCache;
  25 use CXGN::People::PageComment;
  26 use CXGN::People::Person;
  27 use CXGN::Chado::Publication;
  28 use CXGN::Chado::Pubauthor;
  29
  30 use GD;
  31 use GD::Graph::bars;
  32 use GD::Graph::lines;
  33 use GD::Graph::points;
  34 #BEGIN { local $SIG{__WARN__} = sub {}; require GD::Graph::Map }
  35 use Statistics::Descriptive;
  36 use Math::Round::Var;
  37 use Number::Format;
  38 use File::Temp qw / tempfile tempdir /;
  39 use File::Copy;
  40 use File::Spec;
  41 use File::Path qw / mkpath /;
  42 use File::Basename;
  43 use File::Spec::Functions qw / catfile catdir/;
  44 use File::stat;
  45 use File::Slurp qw / read_file /;
  46 use Cache::File;
  47 use Path::Class;
  48 use Try::Tiny;
  49 use CXGN::Contact;
  50 use String::Random;
  51 #use Storable qw / store retrieve /;
  52
  53 use base qw / CXGN::Page::Form::SimpleFormPage CXGN::Phenome::Main/;
  54
  55 use CatalystX::GlobalContext qw( $c );
  56
  57 sub new
  58 {
  59     my $class = shift;
  60     my $self  = $class->SUPER::new(@_);
  61     $self->set_script_name("qtl_analysis.pl");
  62
  63     return $self;
  64 }
  65
  66 sub define_object
  67 {
  68     my $self = shift;
  69
  70     $self->set_dbh( CXGN::DB::Connection->new() );
  71     my %args          = $self->get_args();
  72     my $population_id = $args{population_id};
  73      my $stock_id = $args{stock_id};
  74     my $object;
  75     #########################
  76     # this page needs to be re-written with CXGN::Chado::Stock object
  77     # and without SimpleFormPage, since edits should be done on the parent page only
  78     #########################
  79
  80     unless ( ($population_id and $population_id =~ /^\d+$/) || ($stock_id and $stock_id =~ /^\d+$/) )
  81     {
  82         $c->throw_404("A proper <strong>population id or stock id</strong> argument is missing");
  83     }
  84
  85     if ($stock_id)
  86     {
  87         $object = CXGN::Phenome::Population->new_with_stock_id($self->get_dbh, $stock_id);
  88         $population_id = $object->get_population_id;
  89     } else
  90     {
  91         $object = CXGN::Phenome::Population->new($self->get_dbh, $population_id) ;
  92     }
  93
  94     $self->set_object_id($population_id);
  95     $self->set_object(
  96                        CXGN::Phenome::Population->new(
  97                                         $self->get_dbh(), $self->get_object_id()
  98                        )
  99                      );
 100
 101     $self->set_primary_key("population_id");
 102     $self->set_owners( $self->get_object()->get_owners() );
 103
 104     my $trait_id = $args{cvterm_id};
 105     $trait_id = undef if $trait_id =~ m/\D+/;
 106
 107     if ($trait_id)
 108     {
 109         $self->set_trait_id($trait_id);
 110     } else {
 111         $c->throw_404("A proper <strong>cvterm id</strong> argument is missing");
 112     }
 113
 114
 115 }
 116
 117 sub generate_form
 118 {
 119     my $self = shift;
 120
 121     $self->init_form();
 122
 123     my %args = $self->get_args();
 124
 125     my $population    = $self->get_object();
 126     my $population_id = $self->get_object_id();
 127     my $type_id       = $args{type_id};
 128     my $type          = $args{type};
 129     my $pop_name      = $population->get_name();
 130     my $pop_link =
 131 qq |<a href="/qtl/view/$population_id">$pop_name</a> |;
 132
 133     my $sp_person_id = $population->get_sp_person_id();
 134     my $submitter    = CXGN::People::Person->new( $self->get_dbh(),
 135                                               $population->get_sp_person_id() );
 136     my $submitter_name =
 137       $submitter->get_first_name() . " " . $submitter->get_last_name();
 138     my $submitter_link =
 139 qq |<a href="/solpeople/personal-info.pl?sp_person_id=$sp_person_id">$submitter_name </a> |;
 140
 141     my $login_user    = $self->get_user();
 142     my $login_user_id = $login_user->get_sp_person_id();
 143     my $form          = undef;
 144     if (
 145          $self->get_action() =~ /edit|store/
 146          && (    $login_user_id = $submitter
 147               || $self->get_user()->get_user_type() eq 'curator' )
 148        )
 149     {
 150         $form = CXGN::Page::Form::Editable->new();
 151     }
 152     else
 153     {
 154         $form = CXGN::Page::Form::Static->new();
 155     }
 156
 157     $form->add_label(
 158                       display_name => "Name:",
 159                       field_name   => "name",
 160                       contents     => $pop_link,
 161                     );
 162
 163     $form->add_textarea(
 164                          display_name => "Description: ",
 165                          field_name   => "description",
 166                          object       => $population,
 167                          getter       => "get_description",
 168                          setter       => "set_description",
 169                          columns      => 40,
 170                          rows         => 4,
 171                        );
 172
 173     $form->add_label(
 174                       display_name => "Uploaded by: ",
 175                       field_name   => "submitter",
 176                       contents     => $submitter_link,
 177                     );
 178     $form->add_hidden( field_name => "population_id",
 179                        contents   => $args{population_id} );
 180
 181     $form->add_hidden(
 182                        field_name => "sp_person_id",
 183                        contents   => $self->get_user()->get_sp_person_id(),
 184                        object     => $population,
 185                        setter     => "set_sp_person_id",
 186                      );
 187
 188     $form->add_hidden( field_name => "action", contents => "store" );
 189
 190     $self->set_form($form);
 191
 192     if ( $self->get_action =~ /view|edit/ )
 193     {
 194         $self->get_form->from_database();
 195
 196     }
 197     elsif ( $self->get_action =~ /store/ )
 198     {
 199         $self->get_form->from_request( $self->get_args() );
 200
 201     }
 202
 203 }
 204
 205 sub display_page
 206 {
 207     my $self = shift;
 208     $self->get_page->jsan_use("thickbox");
 209
 210
 211     $self->get_page->add_style( text => <<EOS);
 212 a.abstract_optional_show {
 213   color: blue;
 214   cursor: pointer;
 215   white-space: nowrap;
 216 }
 217 div.abstract_optional_show {
 218   background: #f0f0ff;
 219   border: 1px solid #9F9FC7;
 220   margin: 0.2em 1em 0.2em 1em;
 221   padding: 0.2em 0.5em 0.2em 1em;
 222 }
 223 EOS
 224
 225     my %args            = $self->get_args();
 226     my $dbh             = $self->get_dbh();
 227     my $population      = $self->get_object();
 228     my $population_id   = $self->get_object_id();
 229     my $population_name = $population->get_name();
 230     my $term_id         = $self->get_trait_id();
 231     my $term_name       = $self->get_trait_name();
 232
 233     #used to show certain elements to only the proper users
 234     my $login_user      = $self->get_user();
 235     my $login_user_id   = $login_user->get_sp_person_id();
 236     my $login_user_type = $login_user->get_user_type();
 237
 238     $self->get_page()
 239       ->header(" SGN: $term_name values in population $population_name");
 240     print page_title_html(
 241                     "SGN: $term_name values in population $population_name \n");
 242
 243     my $population_html = $self->get_edit_link_html() . qq |<a href="qtl_form.pl">[New QTL Population]</a><br/>|;
 244
 245     #print all editable form  fields
 246     $population_html .= $self->get_form()->as_table_string();
 247     my $population_obj = $self->get_object();
 248
 249
 250     my $page =
 251 "../phenome/qtl_analysis.pl?population_id=$population_id&amp;cvterm_id=$term_id";
 252     $args{calling_page} = $page;
 253
 254     my $pubmed;
 255     my $url_pubmed = qq | http://www.ncbi.nlm.nih.gov/pubmed/|;
 256
 257     my @publications = $population->get_population_publications();
 258     my $abstract_view;
 259     my $abstract_count = 0;
 260
 261     foreach my $pub (@publications)
 262     {
 263         my (
 264              $title,    $abstract, $authors, $journal,
 265              $pyear,    $volume,   $issue,   $pages,
 266              $obsolete, $pub_id,   $accession
 267            );
 268         $abstract_count++;
 269
 270         my @dbxref_objs = $pub->get_dbxrefs();
 271         my $dbxref_obj  = shift(@dbxref_objs);
 272         $obsolete =
 273           $population_obj->get_population_dbxref($dbxref_obj)->get_obsolete();
 274
 275         if ( $obsolete eq 'f' )
 276         {
 277             $pub_id = $pub->get_pub_id();
 278
 279             $title    = $pub->get_title();
 280             $abstract = $pub->get_abstract();
 281             $pyear    = $pub->get_pyear();
 282             $volume   = $pub->get_volume();
 283             $journal  = $pub->get_series_name();
 284             $pages    = $pub->get_pages();
 285             $issue    = $pub->get_issue();
 286
 287             $accession = $dbxref_obj->get_accession();
 288             my $pub_info =
 289 qq|<a href="/publication/$pub_id/view" >PMID:$accession</a> |;
 290             my @authors;
 291             my $authors;
 292
 293             if ($pub_id)
 294             {
 295
 296                 my @pubauthors_ids = $pub->get_pubauthors_ids($pub_id);
 297
 298                 foreach my $pubauthor_id (@pubauthors_ids)
 299                 {
 300                     my $pubauthor_obj =
 301                       CXGN::Chado::Pubauthor->new( $self->get_dbh,
 302                                                    $pubauthor_id );
 303                     my $last_name   = $pubauthor_obj->get_surname();
 304                     my $first_names = $pubauthor_obj->get_givennames();
 305                     my @first_names = split( /,/, $first_names );
 306                     $first_names = shift(@first_names);
 307                     push @authors, ( "$first_names" . "  " . "$last_name" );
 308                     $authors = join( ", ", @authors );
 309                 }
 310             }
 311
 312             $abstract_view = html_optional_show(
 313                 "abstracts$abstract_count",
 314                 'Abstract',
 315 qq|$abstract <b> <i>$authors.</i> $journal. $pyear. $volume($issue). $pages.</b>|,
 316                 0,                         #< do not show by default
 317                 'abstract_optional_show'
 318                 ,    #< don't use the default button-like style
 319                                                );
 320
 321             $pubmed .=
 322 qq|<div><a href="$url_pubmed$accession" target="blank">$pub_info</a> $title $abstract_view </div> |;
 323         }
 324     }
 325
 326     print info_section_html( title    => 'Population details',
 327                              contents => $population_html, );
 328
 329     my $is_public = $population->get_privacy_status();
 330     my ( $submitter_obj, $submitter_link ) = $self->submitter();
 331
 332     if (    $is_public
 333          || $login_user_type eq 'curator'
 334          || $login_user_id == $population->get_sp_person_id() )
 335     {
 336
 337         my $phenotype = "";
 338         my @phenotype;
 339
 340         my ( $indl_id, $indl_name, $indl_value ) =
 341             $population->get_all_indls_cvterm($term_id);
 342
 343         my ( $min, $max, $avg, $std, $count ) =
 344             $population->get_pop_data_summary($term_id);
 345
 346         for ( my $i = 0 ; $i < @$indl_name ; $i++ )
 347         {
 348
 349             push @phenotype,
 350             [
 351              map { $_ } (
 352                  qq | <a href="/phenome/individual.pl?individual_id=$indl_id->[$i]">$indl_name->[$i]</a>|,
 353                 $indl_value->[$i]
 354              )
 355             ];
 356         }
 357
 358         my ( $phenotype_data, $data_view, $data_download );
 359         my $all_indls_count = scalar(@$indl_name);
 360
 361         if (@phenotype)
 362         {
 363             $phenotype_data = columnar_table_html(
 364                 headings => [
 365                     'Plant accession',
 366                     'Value',
 367
 368                 ],
 369                 data         => \@phenotype,
 370                 __alt_freq   => 2,
 371                 __alt_width  => 1,
 372                 __alt_offset => 3,
 373                 __align      => 'l',
 374                                              );
 375
 376             $data_view = html_optional_show(
 377                       "phenotype",
 378                       'Trait data',
 379                       qq |$phenotype_data|,
 380                       0,                                #<  don't show data by default
 381                                        );
 382             $data_download .=
 383 qq { Download population: <span><a href="/qtl/download/phenotype/$population_id"><b>Phenotype data</b></a> | <a href="/qtl/download/genotype/$population_id"><b>Genotype data</b></a></span> };
 384     }
 385
 386
 387         my (
 388              $image_pheno, $title_pheno, $image_map_pheno,
 389              $plot_html
 390            );
 391         ( $image_pheno, $title_pheno, $image_map_pheno ) = $self->
 392           population_distribution();
 393
 394       $plot_html .= qq | <table  cellpadding = 5><tr><td> |;
 395         $plot_html .= $image_pheno . $image_map_pheno;
 396         $plot_html .= qq | </td><td> |;
 397         $plot_html .= $title_pheno . qq | <br/> |;
 398
 399
 400         my @phe_summ =  ( [ 'No. of obs units', $all_indls_count ],
 401                           [ 'Minimum', $min ],
 402                           [ 'Maximum', $max ],
 403                           [ 'Mean', $avg ],
 404                           [ 'Standard deviation', $std ]
 405                         );
 406
 407         my @summ;
 408         foreach my $phe_summ ( @phe_summ )
 409         {
 410             push @summ, [ map { $_ } ( $phe_summ->[0], $phe_summ->[1] ) ];
 411         }
 412
 413         my $summ_data  = columnar_table_html(
 414                                               headings   => [ '',  ''],
 415                                               data       => \@summ,
 416                                               __alt_freq   => 2,
 417                                               __alt_width  => 1,
 418                                               __alt_offset => 3,
 419                                               __align      => 'l',
 420                                             );
 421
 422
 423         $plot_html .= $summ_data;
 424         $plot_html .= qq | </td></tr></table> |;
 425
 426
 427
 428         my ( $qtl_image, $legend);
 429
 430         #using standard deviation of 0.01 as an arbitrary cut off to run
 431         #qtl analysis. Probably, need to think of better solution.
 432         if ( $std >= 0.01 ) {
 433             $qtl_image           = $self->qtl_plot();
 434             $legend = $self->legend();
 435         }
 436         else {
 437             $qtl_image = 'There is no statistically significant phenotypic
 438                           variation for this trait to run
 439                           QTL analysis.';
 440         }
 441
 442
 443         my $qtl_html = qq | <table><tr><td width=70%>$qtl_image</td><td width=30%>$legend</td></tr></table> |;
 444         my $qtl_effects_ref = $self->qtl_effects();
 445         my $explained_variation = $self->explained_variation();
 446         my ($qtl_effects_data, $explained_variation_data);
 447
 448         if ($qtl_effects_ref)
 449         {
 450             $qtl_effects_data  = columnar_table_html(
 451                                               data         => $qtl_effects_ref,
 452                                               __alt_freq   => 2,
 453                                               __alt_width  => 1,
 454                                               __alt_offset => 3,
 455                                               __align      => 'l',
 456                                             );
 457
 458
 459         } else
 460         {
 461             $qtl_effects_data = "No QTL effects estimates were found for QTL(s) of  this trait.";
 462         }
 463
 464         if ($explained_variation) {
 465          $explained_variation_data  = columnar_table_html(
 466                                               data         => $explained_variation,
 467                                               __alt_freq   => 2,
 468                                               __alt_width  => 1,
 469                                               __alt_offset => 3,
 470                                               __align      => 'l',
 471                                             );
 472
 473         } else  {
 474             $explained_variation_data = "No explained variation estimates were found for QTL(s) of this trait.";
 475         }
 476
 477         print info_section_html(
 478                                 title    => 'QTL(s)',
 479                                 contents => $qtl_html,
 480                                 );
 481
 482
 483         print info_section_html( title    => 'Variation explained by QTL(s) ( Interacting QTLs model )',
 484                                  contents => $explained_variation_data
 485                                  );
 486
 487         print info_section_html( title    => 'QTL effects',
 488                                  contents => $qtl_effects_data
 489                                  );
 490
 491         print info_section_html(
 492                                 title    => 'Phenotype frequency distribution',
 493                                 contents => $plot_html,
 494                                );
 495
 496         print info_section_html(
 497                                 title    => 'Download data',
 498                                 contents => $data_view . " " . $data_download,
 499                                 collapsible => 1,
 500                                 collapsed   => 1,
 501                                );
 502
 503     }
 504     else
 505     {
 506         my $message =
 507           "The QTL data for this trait in this population is not public yet.
 508                        If you would like to know more about this data,
 509                        please contact the owner of the data: <b>$submitter_link</b>
 510                        or email to SGN:
 511                        <a href=mailto:sgn-feedback\@sgn.cornell.edu>
 512                        sgn-feedback\@sgn.cornell.edu</a>.\n";
 513
 514         print info_section_html( title    => 'QTL(s)',
 515                                  contents => $message );
 516     }
 517
 518     print info_section_html(
 519                            title       => 'Publication(s)',
 520                            contents    => $pubmed,
 521                            collapsible => 1,
 522                            collapsed   => 1,
 523                            );
 524
 525     if ($population_name)
 526     {
 527         my $page_comment_obj =
 528             CXGN::People::PageComment->new( $self->get_dbh(), "population",
 529                                             $population_id,
 530                                             $self->get_page()->{request}->uri()."?".$self->get_page()->{request}->args()
 531             );
 532         print $page_comment_obj->get_html();
 533     }
 534
 535
 536     $self->get_page()->footer();
 537     $self->remove_permu_file();
 538
 539     exit();
 540 }
 541
 542 # override store to check if a locus with the submitted symbol/name already exists in the database
 543
 544  sub store
 545  {
 546      my $self          = shift;
 547      my $population    = $self->get_object();
 548      my $population_id = $self->get_object_id();
 549      my %args          = $self->get_args();
 550
 551      $self->SUPER::store(0);
 552  }
 553
 554 sub population_distribution
 555 {
 556     my $self      = shift;
 557     my $pop_id    = $self->get_object_id();
 558     my $pop       = $self->get_object();
 559     my $pop_name  = $pop->get_name();
 560     my $term_name = $self->get_trait_name();
 561     my $term_id   = $self->get_trait_id();
 562     my $dbh       = $self->get_dbh();
 563
 564
 565     my $basepath     = $c->get_conf("basepath");
 566     my $tempfile_dir = $c->get_conf("tempfiles_subdir");
 567
 568     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
 569     $self->cache_temp_path();
 570
 571     my $cache = CXGN::Tools::WebImageCache->new();
 572     $cache->set_basedir($basepath);
 573     $cache->set_temp_dir( $tempfile_dir . '/temp_images');
 574     $cache->set_expiration_time(259200);
 575     $cache->set_key( "popluation_distribution" . $pop_id . $term_id );
 576     $cache->set_map_name("popmap$pop_id$term_id");
 577
 578     my ( @value,    @indl_id, @indl_name );
 579
 580     $cache->set_force(0);
 581     if ( !$cache->is_valid() )
 582     {
 583
 584         my ( $indl_id, $indl_name, $value ) = $pop->plot_cvterm($term_id);
 585         @value = @{$value};
 586
 587         my $stat = Statistics::Descriptive::Full->new();
 588         $stat->add_data(@value);
 589         my %f = $stat->frequency_distribution(10);
 590
 591         my ( @keys, @counts );
 592
 593         for ( sort { $a <=> $b } keys %f )
 594         {
 595             my $round = Math::Round::Var->new(0.001);
 596             my $key = $round->round($_);
 597             push @keys,   $key;
 598             push @counts, $f{$_};
 599         }
 600
 601         my $min = $stat->min();
 602         if ( $min != 0 )
 603         {
 604             $min = $min - 0.01;
 605         }
 606
 607         my @keys_range = $min . '-' . $keys[0];
 608
 609         my $range;
 610         my $previous_k   = $keys[0];
 611         my $keys_shifted = shift(@keys);
 612
 613         foreach my $k (@keys)
 614         {
 615             $range = $previous_k . '-' . $k;
 616             push @keys_range, $range;
 617             $previous_k = $k;
 618         }
 619
 620         my $max = $counts[0];
 621         foreach my $i ( @counts[ 1 .. $#counts ] )
 622         {
 623             if ( $i > $max ) { $max = $i; }
 624         }
 625         $max = int( $max + ( $max * 0.1 ) );
 626
 627         my $c_html;
 628         my @c_html;
 629         my ( $lower, $upper );
 630
 631         foreach my $k (@keys_range)
 632         {
 633             ( $lower, $upper ) = split( /-/, $k );
 634             $c_html =
 635 qq | /phenome/indls_range_cvterm.pl?cvterm_id=$term_id&amp;lower=$lower&amp;upper=$upper&amp;population_id=$pop_id |;
 636             push @c_html, $c_html;
 637
 638         }
 639
 640         my @bar_clr = ("orange");
 641         my @data    = ( [@keys_range], [@counts] );
 642         my $graph   = new GD::Graph::bars();
 643
 644         $graph->set_title_font('gdTinyFont');
 645         $graph->set(
 646                      title             => " ",
 647                      x_label           => "$term_name values",
 648                      y_label           => "No. of observation units",
 649                      y_max_value       => $max,
 650                      x_all_ticks       => 1,
 651                      y_all_ticks       => 2,
 652                      y_label_skip      => 5,
 653                      y_plot_values     => 0,
 654                      x_label_skip      => 1,
 655                      width             => 400,
 656                      height            => 400,
 657                      bar_width         => 30,
 658                      x_labels_vertical => 1,
 659                      show_values       => 1,
 660                      textclr           => "black",
 661                      dclrs             => \@bar_clr,
 662                    );
 663
 664         $cache->set_image_data( $graph->plot( \@data )->png );
 665
 666         my $map = new GD::Graph::Map(
 667                                       $graph,
 668                                       hrefs       => [ \@c_html ],
 669                                       noImgMarkup => 1,
 670                                       mapName     => "popmap$pop_id$term_id",
 671                                       info        => "%x: %y lines",
 672                                     );
 673         $cache->set_image_map_data(
 674                       $map->imagemap( "popimage$pop_id$term_id.png", \@data ) );
 675
 676     }
 677
 678     my $image_map = $cache->get_image_map_data();
 679     my $image     = $cache->get_image_tag();
 680     my $title =
 681 qq | Phenotype frequency distribution of experimental lines evaluated for $term_name. Bars represent the number of experimental lines with $term_name values greater than the lowest value but less or equal to the highest value of the range. |;
 682
 683     return $image, $title, $image_map;
 684 }
 685
 686 sub qtl_plot
 687 {
 688     my $self           = shift;
 689     my $dbh            = $self->get_dbh();
 690     my $pop_id         = $self->get_object_id();
 691     my $population     = $self->get_object();
 692     my $pop_name       = $population->get_name();
 693     my $mapversion     = $population->mapversion_id();
 694     my @linkage_groups = $population->linkage_groups();
 695     my $term_name      = $self->get_trait_name();
 696     my $term_id        = $self->get_trait_id();
 697     my $ac             = $population->cvterm_acronym($term_name);
 698     my $basepath       = $c->get_conf("basepath");
 699     my $tempfile_dir   = $c->get_conf("tempfiles_subdir");
 700
 701     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) = $self->cache_temp_path();
 702     my $cache_tempimages = Cache::File->new( cache_root => $tempimages_path );
 703     $cache_tempimages->purge();
 704
 705     my ( @marker,  @chr,  @pos,   @lod, @chr_qtl, @peak_markers );
 706     my ( $qtl_image, $image, $image_t, $image_url, $image_html, $image_t_url,
 707          $thickbox, $title, $l_m, $p_m, $r_m );
 708
 709     my $round       = Number::Format->new();
 710     $qtl_image      = $self->qtl_images_exist();
 711     my $permu_data  = $self->permu_file();
 712     my $stat_option = $self->qtl_stat_option();
 713
 714     unless ( $qtl_image && -s $permu_data > 1 && $stat_option=~/default/)
 715     {
 716         my ( $qtl_summary, $peak_markers_file ) = $self->run_r();
 717
 718         open my $qtl_fh, "<", $qtl_summary or die "can't open $qtl_summary: $!\n";
 719
 720         my $header = <$qtl_fh>;
 721         while ( my $row = <$qtl_fh> )
 722         {
 723             my ( $marker, $chr, $pos, $lod ) = split( /\t/, $row );
 724             push @marker, $marker;
 725             push @chr,    $chr;
 726             push @pos, $round->round($pos, 1);
 727             push @lod, $round->round($lod, 2);
 728         }
 729
 730         my @o_lod = sort(@lod);
 731         my $max   = $o_lod[-1];
 732         $max = $max + (0.5);
 733
 734
 735         open my $markers_fh, "<", $peak_markers_file
 736            or die "can't open $peak_markers_file: !$\n";
 737
 738         $header = <$markers_fh>;
 739         while ( my $row = <$markers_fh> )
 740
 741         {
 742             chomp($row);
 743             my ($trash, $chr_qtl, $peak_marker ) = split( /\t/, $row );
 744             push @chr_qtl, $chr_qtl;
 745             push @peak_markers, $peak_marker;
 746         }
 747
 748         my (@h_markers, @chromosomes, @lk_groups);
 749         my $h_marker;
 750
 751         @lk_groups  = @linkage_groups;
 752         @lk_groups  = sort ( { $a <=> $b } @lk_groups );
 753         my $random  = String::Random->new();
 754        # my $user_id = $c->user->get_object->get_sp_person_id() if $c->user;
 755        # my $qtl_obj = CXGN::Phenome::Qtl->new($user_id);
 756        # my ($user_qtl_dir, $user_dir) = $qtl_obj->get_user_qtl_dir($c);
 757
 758         for ( my $i = 0 ; $i < @linkage_groups ; $i++ )
 759         {
 760             my $lg = shift(@lk_groups);
 761             my $key_h_marker;
 762             #my %keys_to_urls=();
 763
 764             if ($self->qtl_stat_option() eq 'user params')
 765             {
 766                 $key_h_marker = $random->randpattern("CCccCCnnn") . '_'. $lg;
 767                # %keys_to_urls = ("user_params_${lg}" => $key_h_marker);
 768                # store(\%keys_to_urls, "$user_dir/image_urls");
 769
 770             }
 771             elsif ( $self->qtl_stat_option() eq 'default')
 772
 773             {
 774                 $key_h_marker = "$ac" . "_pop_" . "$pop_id" . "_chr_" . $lg;
 775             }
 776
 777             $h_marker = $cache_tempimages->get($key_h_marker)
 778                       ? $self->qtl_stat_option eq 'default'
 779                       : undef
 780                       ;
 781
 782             unless ($h_marker)
 783             {
 784                 push @chromosomes, $lg;
 785                 $p_m = $peak_markers[$i];
 786                 $p_m =~ s/\s//;
 787
 788                 my $permu_threshold_ref = $self->permu_values();
 789                 my %permu_threshold     = %$permu_threshold_ref;
 790                 my @p_keys;
 791                 foreach my $key ( keys %permu_threshold )
 792                 {
 793                     if ( $key =~ m/^\d./ )
 794                     {
 795                         push @p_keys, $key;
 796                     }
 797
 798                 }
 799                 my $lod1 = $permu_threshold{$p_keys[0]};
 800
 801                 $h_marker =
 802                     qq |/phenome/qtl.pl?population_id=$pop_id&amp;term_id=$term_id&amp;chr=$lg&amp;peak_marker=$p_m&amp;lod=$lod1|;
 803
 804                 $cache_tempimages->set( $key_h_marker, $h_marker, '30 days' );
 805                 push @h_markers, $h_marker;
 806             }
 807         }
 808
 809         my $count       = 0;
 810         my $old_chr_chr = 1;
 811         my (
 812              $chr_chr, $image,     $image_t,
 813              $thickbox,  $max_chr, $chr_chr_e, $marker_chr_e,
 814              $pos_chr_e, $lod_chr_e
 815            );
 816         my $chrs = ( scalar(@chromosomes) ) + 1;
 817
 818         for ( my $i = 1 ; $i < $chrs ; $i++ )
 819         {
 820             my ( @marker_chr, @chr_chr, @pos_chr, @lod_chr, @data, @m_html ) =
 821               ();
 822             my ( $marker_chr, $pos_chr, $lod_chr, $max_chr );
 823
 824             $h_marker = shift(@h_markers);
 825
 826             if ( ( $i == $old_chr_chr ) && ( $i != 12 ) )
 827             {
 828                 push @marker_chr, $marker_chr_e;
 829                 push @chr_chr,    $chr_chr_e;
 830                 push @pos_chr, $round->round($pos_chr_e, 1);
 831                 push @lod_chr, $round->round($lod_chr_e, 2);
 832             }
 833
 834             my $cache_qtl_plot = CXGN::Tools::WebImageCache->new();
 835            $cache_qtl_plot->set_basedir($basepath);
 836            $cache_qtl_plot->set_temp_dir( $tempfile_dir . "/temp_images" );
 837 #           $cache_qtl_plot->set_temp_dir( $tempimages_path);
 838           $cache_qtl_plot->set_expiration_time(259200);
 839
 840
 841            if ($self->qtl_stat_option() eq 'user params')
 842             {
 843                 $cache_qtl_plot->set_key($random->randpattern("CCccCCnnn"));
 844                 $cache_qtl_plot->set_force(1);
 845            }
 846            elsif ( $self->qtl_stat_option() eq 'default')
 847            {
 848                 $cache_qtl_plot->set_key("qtlplot" . $i . "small" . $pop_id . $term_id);
 849                 $cache_qtl_plot->set_force(0);
 850            }
 851
 852             if ( !$cache_qtl_plot->is_valid() )
 853             {
 854
 855                 for ( my $j = 0 ; $j < @marker ; $j++ )
 856                 {
 857
 858                     $chr_chr = $chr[$j];
 859
 860                     if ( $i == $chr_chr )
 861                     {
 862                         $marker_chr = $marker[$j];
 863
 864                         $pos_chr = $pos[$j];
 865                         $lod_chr = $lod[$j];
 866
 867                         push @marker_chr, $marker_chr;
 868                         push @chr_chr,    $chr_chr;
 869                         push @pos_chr, $round->round($pos_chr, 1);
 870                         push @lod_chr, $round->round($lod_chr, 2);
 871
 872                         ( $chr_chr_e, $marker_chr_e, $pos_chr_e, $lod_chr_e ) =
 873                           ();
 874                     }
 875
 876                     elsif ( $i != $chr_chr )
 877                     {
 878
 879                         $chr_chr_e    = $chr[$j];
 880                         $marker_chr_e = $marker[$j];
 881                         $pos_chr_e    = $pos[$j];
 882                         $lod_chr_e    = $lod[$j];
 883                     }
 884
 885                 }
 886
 887                 @data = ( [ (@pos_chr) ], [@lod_chr] );
 888                 my $graph = new GD::Graph::lines( 110, 110 );
 889                 $graph->set_title_font('gdTinyFont');
 890                 $graph->set(
 891                              title             => " ",
 892                              x_label           => "Chr $i (cM)",
 893                              y_label           => "LOD",
 894                              y_max_value       => 10,
 895                              x_all_ticks       => 5,
 896                              y_all_ticks       => 1,
 897                              y_label_skip      => 1,
 898                              y_plot_values     => 1,
 899                              x_label_skip      => 5,
 900                              x_plot_values     => 1,
 901                              x_labels_vertical => 1,
 902                              textclr           => "black"
 903                            );
 904
 905                 $cache_qtl_plot->set_image_data( $graph->plot( \@data )->png );
 906
 907             }
 908
 909             $image      = $cache_qtl_plot->get_image_tag();
 910             $image_url  = $cache_qtl_plot->get_image_url();
 911
 912
 913 ###########thickbox
 914             my $cache_qtl_plot_t = CXGN::Tools::WebImageCache->new();
 915             $cache_qtl_plot_t->set_basedir($basepath);
 916             $cache_qtl_plot_t->set_temp_dir( $tempfile_dir . "/temp_images" );
 917            # $cache_qtl_plot_t->set_temp_dir( $tempimages_path);
 918             $cache_qtl_plot_t->set_expiration_time(259200);
 919
 920            if ($self->qtl_stat_option() eq 'user params')
 921            {
 922                 $cache_qtl_plot_t->set_key($random->randpattern("CCccccnnn"));
 923                 $cache_qtl_plot_t->set_force(1);
 924            }
 925
 926             elsif ( $self->qtl_stat_option() eq 'default')
 927
 928            {
 929                 $cache_qtl_plot_t->set_key("qtlplot_" . $i . "_thickbox_" . $pop_id . $term_id);
 930                 $cache_qtl_plot_t->set_force(0);
 931            }
 932
 933
 934             if ( !$cache_qtl_plot_t->is_valid() )
 935             {
 936                 my @o_lod_chr = sort { $a <=> $b } @lod_chr;
 937                 $max_chr = pop(@o_lod_chr);
 938                 $max_chr = $max_chr + (0.5);
 939
 940                 my $graph_t = new GD::Graph::lines( 420, 420 );
 941                 $graph_t->set_title_font('gdTinyFont');
 942                 $graph_t->set(
 943                                title             => " ",
 944                                x_label           => "Chromosome $i (cM)",
 945                                y_label           => "LOD",
 946                                y_max_value       => $max_chr,
 947                                x_all_ticks       => 5,
 948                                y_all_ticks       => 1,
 949                                y_label_skip      => 1,
 950                                y_plot_values     => 1,
 951                                x_label_skip      => 5,
 952                                x_plot_values     => 1,
 953                                x_labels_vertical => 1,
 954                                textclr           => "black"
 955                              );
 956
 957                 $cache_qtl_plot_t->set_image_data(
 958                                                 $graph_t->plot( \@data )->png );
 959
 960             }
 961
 962             $image_t     = $cache_qtl_plot_t->get_image_tag();
 963             $image_t_url = $cache_qtl_plot_t->get_image_url();
 964
 965             $thickbox =
 966 qq | <a href="$image_t_url" title="<a href=$h_marker&amp;qtl=$image_t_url><font color=#f87431><b>>>>>Go to the QTL page>>>> </b></font></a>" class="thickbox" rel="gallary-qtl"> <img src="$image_url" alt="Chromosome $i $image_t_url $image_url" /> </a> |;
 967
 968             $qtl_image .= $thickbox;
 969             $title       = "  ";
 970             $old_chr_chr = $chr_chr;
 971         }
 972     }
 973
 974     return $qtl_image;
 975 }
 976
 977 =head2 infile_list
 978  Usage: my $file_in = $self->infile_list();
 979  Desc: returns an R input tempfile containing a tempfile
 980        holding the cvterm acronym, pop id, a filepath to the phenotype dataset file,
 981         a filepath to genotype dataset file, a filepath to the permuation file.
 982  Ret: an R input tempfile name (with abosulte path)
 983  Args:
 984  Side Effects:
 985  Example:
 986
 987
 988 =cut
 989
 990 sub infile_list
 991 {
 992     my $self       = shift;
 993     my $dbh        = $self->get_dbh();
 994     my $pop_id     = $self->get_object_id();
 995     my $population = $self->get_object();
 996     my $term_name  = $self->get_trait_name();
 997     my $term_id    = $self->get_trait_id();
 998
 999     my $ac = $population->cvterm_acronym($term_name);
1000
1001     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1002       $self->cache_temp_path();
1003
1004     my $prod_permu_file  = $self->permu_file();
1005     my $gen_dataset_file = $population->genotype_file($c);
1006     my $phe_dataset_file = $population->phenotype_file($c);
1007     my $crosstype_file   = $self->crosstype_file();
1008     my $stat_file = $self->stat_files();
1009
1010     my $input_file_list_temp =
1011       File::Temp->new(
1012                        TEMPLATE => "infile_list_${ac}_$pop_id-XXXXXX",
1013                        DIR      => $prod_temp_path,
1014                        UNLINK   => 0,
1015                      );
1016     my $file_in = $input_file_list_temp->filename();
1017
1018     my $file_cvin = File::Temp->new(
1019                                      TEMPLATE => 'cvterm_input-XXXXXX',
1020                                      DIR      => $prod_temp_path,
1021                                      UNLINK   => 0,
1022                                    );
1023     my $file_cv_in = $file_cvin->filename();
1024
1025     open my $cv_fh, ">", $file_cv_in or die "can't open $file_cv_in: $!\n";
1026     $cv_fh->print($ac);
1027
1028     my $popid_temp = File::Temp->new(
1029                                      TEMPLATE => 'popid-XXXXXX',
1030                                      DIR      => $prod_temp_path,
1031                                      UNLINK   => 0,
1032                                    );
1033     my $file_popid = $popid_temp->filename();
1034
1035     open my $popid_fh, ">", $file_popid or die "can't open $file_popid: $!\n";
1036     $popid_fh->print($pop_id);
1037
1038     my $file_in_list = join( "\t",
1039                              $file_cv_in,       $file_popid,
1040                              $gen_dataset_file, $phe_dataset_file,
1041                              $prod_permu_file, $crosstype_file, $stat_file);
1042
1043     open my $fi_fh, ">", $file_in or die "can't open $file_in: $!\n";
1044     $fi_fh->print ($file_in_list);
1045
1046     return $file_in;
1047
1048 }
1049
1050 =head2 outfile_list
1051
1052  Usage: my ($file_out, $qtl_summary, $peak_markers) = $self->outfile_list();
1053  Desc: returns an R output tempfile containing a tempfile supposed to hold the qtl
1054        mapping output and another tempfile for the qtl peak markers
1055        and the qtl mapping output and qtl peak makers files separately
1056        (convenient for reading their data when plotting the qtl)
1057  Ret: R output file names (with abosulte path)
1058  Args:
1059  Side Effects:
1060  Example:
1061
1062 =cut
1063
1064 sub outfile_list
1065 {
1066     my $self       = shift;
1067     my $pop_id     = $self->get_object_id();
1068     my $population = $self->get_object();
1069     my $term_id    = $self->get_trait_id();
1070     my $term_name  = $self->get_trait_name();
1071     my $dbh        = $self->get_dbh();
1072
1073     my $ac = $population->cvterm_acronym($term_name);
1074
1075     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1076       $self->cache_temp_path();
1077
1078     my $output_file_list_temp =
1079       File::Temp->new(
1080                        TEMPLATE => "outfile_list_${ac}_$pop_id-XXXXXX",
1081                        DIR      => $prod_temp_path,
1082                        UNLINK   => 0,
1083                      );
1084     my $file_out = $output_file_list_temp->filename();
1085
1086     my $qtl_temp = File::Temp->new(
1087                                  TEMPLATE => "qtl_summary_${ac}_$pop_id-XXXXXX",
1088                                  DIR      => $prod_cache_path,
1089                                  UNLINK   => 0
1090     );
1091     my $qtl_summary = $qtl_temp->filename;
1092
1093     my $marker_temp = File::Temp->new(
1094                             TEMPLATE => "peak_markers_${ac}_$pop_id-XXXXXX",
1095                             DIR      => $prod_cache_path,
1096                             UNLINK   => 0
1097     );
1098
1099     my $peak_markers = $marker_temp->filename;
1100
1101     my $ci_lod = $population->ci_lod_file($c, $ac);
1102     my $qtl_effects = $population->qtl_effects_file($c, $ac);
1103     my $explained_variation = $population->explained_variation_file($c, $ac);
1104     my $file_out_list = join ( "\t"
1105         ,$qtl_summary
1106         ,$peak_markers
1107         ,$ci_lod
1108         ,$qtl_effects
1109         ,$explained_variation
1110         );
1111
1112     open my $fo_fh, ">", $file_out or die "can't open $file_out: $!\n";
1113     $fo_fh->print ($file_out_list);
1114
1115     return $file_out, $qtl_summary, $peak_markers;
1116 }
1117
1118 =head2 cache_temp_path
1119
1120  Usage: my ($solqtl_cache, $solqtl_tempfiles, $solqtl_temp_images) = $self->cache_temp_path();
1121  Desc: creates the 'solqtl/cache', 'solqtl/tempfiles', and 'solqtl/temp_images', subdirs in the /export/prod/tmp,
1122  Ret: returns the dirs above
1123  Args: none
1124  Side Effects:
1125  Example:
1126
1127 =cut
1128
1129 sub cache_temp_path {
1130     my $geno_version = $c->config->{default_genotyping_protocol};
1131     $geno_version    = 'analysis-data' if ($geno_version =~ /undefined/) || !$geno_version;
1132     $geno_version    =~ s/\s+//g;
1133     my $tmp_dir      = $c->site_cluster_shared_dir;
1134     $tmp_dir         = catdir($tmp_dir, $geno_version);
1135
1136     my $solqtl_dir         = catdir($tmp_dir, 'solqtl');
1137     my $solqtl_cache       = catdir($tmp_dir, 'solqtl', 'cache');
1138     my $solqtl_tempfiles   = catdir($tmp_dir, 'solqtl', 'tempfiles');
1139
1140     my $basepath = $c->config->{basepath};
1141     my $temp_dir  = $c->config->{tempfiles_subdir};
1142
1143     my $tempimages   = catdir($basepath, $temp_dir, "temp_images" );
1144
1145     mkpath ([$solqtl_cache, $solqtl_tempfiles, $tempimages], 0, 0755);
1146     return $solqtl_cache, $solqtl_tempfiles, $tempimages;
1147
1148
1149 }
1150
1151
1152
1153 =head2 crosstype_file
1154
1155  Usage: my $cross_file = $self->crosstype_file();
1156  Desc: creates the crosstype file in the /data/prod/tmp/r_qtl/tempfiles,
1157
1158  Ret: crosstype filename (with absolute path)
1159  Args: none
1160  Side Effects:
1161  Example:
1162
1163 =cut
1164
1165 sub crosstype_file {
1166     my $self       = shift;
1167     my $pop_id     = $self->get_object_id();
1168     my $population = $self->get_object();
1169
1170     my $type_id = $population->get_cross_type_id
1171         or die "population '$pop_id' has no cross_type, does not seem to be the product of a cross!";
1172     my ($cross_type) = $self->get_dbh->selectrow_array(<<'', undef, $type_id);
1173                           select cross_type
1174                           from cross_type
1175                           where cross_type_id = ?
1176
1177     my $rqtl_cross_type = { 'Back cross'        => 'bc',
1178                             'F2'                => 'f2',
1179                             'RIL (selfing)' => 'rilself',
1180                             'RIL (sibling mating)'  => 'rilsib'
1181                           }->{$cross_type}
1182                           or die "unknown cross_type '$cross_type' for population '$pop_id'";
1183
1184     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1185         $self->cache_temp_path();
1186
1187     my $cross_temp = File::Temp->new(
1188         TEMPLATE => "cross_type_${pop_id}-XXXXXX",
1189         DIR      => $prod_temp_path,
1190         UNLINK   => 0,
1191        );
1192
1193     $cross_temp->print( $rqtl_cross_type );
1194     return $cross_temp->filename;
1195 }
1196
1197 =head2 run_r
1198
1199  Usage: my ($qtl_summary, $peak_markers) = $self->run_r();
1200  Desc: run R in the cluster; returns the R output files (with abosulate filepath) with qtl mapping data
1201        and peak markers
1202  Ret:
1203  Args:  none
1204  Side Effects:
1205  Example:
1206
1207 =cut
1208
1209 sub run_r {
1210     my $self = shift;
1211
1212     my ($solqtl_cache, $solqtl_temp, $solqtl_tempimages) = $self->cache_temp_path();
1213     my $prod_permu_file = $self->permu_file();
1214     my $input_file      = $self->infile_list();
1215     my ($output_file, $qtl_summary, $peak_markers) = $self->outfile_list();
1216
1217     my $pop_id = $self->get_object_id();
1218     my $trait_id = $self->get_trait_id();
1219
1220     $c->stash->{analysis_tempfiles_dir} = $solqtl_temp;
1221     $c->stash->{r_temp_file} = "solqtl-${pop_id}-${trait_id}";
1222     $c->stash->{input_files}  = $input_file;
1223     $c->stash->{output_files} = $output_file;
1224     $c->stash->{r_script}    = 'R/solGS/qtl_analysis.r';
1225
1226     $c->controller('solGS::AsyncJob')->run_r_script($c);
1227
1228     return $qtl_summary, $peak_markers;
1229
1230 }
1231
1232 =head2 permu_file
1233
1234  Usage: my $permu_file = $self->permu_file();
1235  Desc: creates the permutation file in the /data/prod/tmp/r_qtl/cache,
1236        if it does not exist yet and caches it for R.
1237  Ret: permutation filename (with abosolute path)
1238  Args: none
1239  Side Effects:
1240  Example:
1241
1242 =cut
1243
1244 sub permu_file
1245 {
1246     my $self       = shift;
1247     my $dbh        = $self->get_dbh();
1248     my $pop_id     = $self->get_object_id();
1249     my $population = $self->get_object();
1250     my $pop_name   = $population->get_name();
1251     my $term_name  = $self->get_trait_name();
1252     my $term_id    = $self->get_trait_id();
1253
1254     my $ac = $population->cvterm_acronym($term_name);
1255
1256     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1257       $self->cache_temp_path();
1258
1259     my $file_cache = Cache::File->new( cache_root => $prod_cache_path );
1260
1261     my $key_permu = "$ac" . "_" . $pop_id . "_permu";
1262     my $filename  = "permu_" . $ac . "_" . $pop_id;
1263
1264     my  $permu_file = $file_cache->get($key_permu);
1265
1266     unless ($permu_file)
1267     {
1268
1269         my $permu = undef;
1270
1271         my $permu_file = File::Spec->catfile( $prod_cache_path, $filename );
1272
1273         open my $permu_fh, ">", $permu_file or die "can't open $permu_file: !$\n";
1274         $permu_fh->print($permu);
1275
1276
1277         $file_cache->set( $key_permu, $permu_file, '30 days' );
1278         $permu_file = $file_cache->get($key_permu);
1279     }
1280
1281     return $permu_file;
1282
1283 }
1284
1285 =head2 permu_values
1286
1287  Usage: my $permu_values = $self->permu_values();
1288  Desc: reads the permutation output from R,
1289        creates a hash with the probality level as key and LOD threshold as the value,
1290
1291  Ret: a hash ref of the permutation values
1292  Args: none
1293  Side Effects:
1294  Example:
1295
1296 =cut
1297
1298 sub permu_values
1299 {
1300     my $self            = shift;
1301     my $prod_permu_file = $self->permu_file();
1302
1303     my %permu_threshold;
1304
1305     my $permu_file = fileparse($prod_permu_file);
1306     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1307       $self->cache_temp_path();
1308
1309     $permu_file = File::Spec->catfile( $prod_cache_path, $permu_file );
1310
1311     my $round1 = Math::Round::Var->new(0.1);
1312
1313     open my $permu_fh, "<", $permu_file
1314       or die "can't open $permu_file: !$\n";
1315
1316     my $header = <$permu_fh>;
1317
1318     while ( my $row = <$permu_fh> )
1319     {
1320         my ( $significance, $lod_threshold ) = split( /\t/, $row );
1321         $lod_threshold = $round1->round($lod_threshold);
1322         $permu_threshold{$significance} = $lod_threshold;
1323     }
1324
1325     return \%permu_threshold;
1326
1327 }
1328
1329
1330
1331 =head2 qtl_images_exist
1332
1333  Usage: my $qtl_images_ref = $self->qtl_images_exist();
1334  Desc: checks and returns a scalar ref if the qtl plots (with thickbox and their links to the comparative viewer) exist in the cache
1335  Ret: scalar ref to the images or undef
1336  Args: none
1337  Side Effects:
1338  Example:
1339
1340 =cut
1341
1342 sub qtl_images_exist
1343 {
1344     my $self       = shift;
1345     my $pop_id     = $self->get_object_id();
1346     my $population = $self->get_object();
1347     my $pop_name   = $population->get_name();
1348     my $term_name  = $self->get_trait_name();
1349     my $term_id    = $self->get_trait_id();
1350     my $dbh        = $self->get_dbh();
1351     my $qtl_image  = undef;
1352
1353     my @linkage_groups = $population->linkage_groups();
1354     @linkage_groups = sort ( { $a <=> $b } @linkage_groups );
1355
1356     my $ac = $population->cvterm_acronym($term_name);
1357
1358     my $basepath     = $c->get_conf("basepath");
1359     my $tempfile_dir = $c->get_conf("tempfiles_subdir");
1360
1361     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1362             $self->cache_temp_path();
1363
1364     my $cache_tempimages = Cache::File->new( cache_root => $tempimages_path );
1365     $cache_tempimages->purge();
1366
1367     my $cache_qtl_plot = CXGN::Tools::WebImageCache->new();
1368     $cache_qtl_plot->set_basedir($basepath);
1369     $cache_qtl_plot->set_temp_dir( $tempfile_dir . "/temp_images" );
1370
1371     if ($self->qtl_stat_option eq 'default')
1372     {
1373         my ( $image, $image_t, $image_url, $image_html, $image_t_url,
1374              $thickbox, $title );
1375
1376
1377       IMAGES: foreach my $lg (@linkage_groups)
1378       {
1379           my $key_h_marker = "$ac" . "_pop_" . "$pop_id" . "_chr_" . $lg;
1380           my $h_marker     = $cache_tempimages->get($key_h_marker);
1381
1382           my $key = "qtlplot" . $lg . "small" . $pop_id . $term_id;
1383           $cache_qtl_plot->set_key($key);
1384
1385           if ( $cache_qtl_plot->is_valid )
1386           {
1387               $image      = $cache_qtl_plot->get_image_tag();
1388               $image_url  = $cache_qtl_plot->get_image_url();
1389           }
1390
1391           my $key_t = "qtlplot_" . $lg . "_thickbox_" . $pop_id . $term_id;
1392           $cache_qtl_plot->set_key($key_t);
1393
1394           if ( $cache_qtl_plot->is_valid )
1395           {
1396               $image_t     = $cache_qtl_plot->get_image_tag();
1397               $image_t_url = $cache_qtl_plot->get_image_url();
1398
1399               $thickbox =
1400 qq | <a href="$image_t_url" title= "<a href=$h_marker&amp;qtl=$image_t_url><font color=#f87431><b>>>>Go to the QTL page>>>> </b></font></a>"  class="thickbox" rel="gallary-qtl"> <img src="$image_url" alt="Chromosome $lg $image_t_url $image_url" /> </a> |;
1401
1402               $qtl_image .= $thickbox;
1403               $title = "  ";
1404           }
1405           else
1406           {
1407               $qtl_image = undef;
1408               last IMAGES;
1409
1410           }
1411       }
1412     }
1413     else
1414     {
1415
1416         foreach my $lg (@linkage_groups)
1417         {
1418             my $key_h_marker = $ac . "_pop_" . $pop_id . "_chr_" . $lg;
1419             $cache_tempimages->remove($key_h_marker);
1420
1421             my $key = "qtlplot" . $lg . "small" . $pop_id . $term_id;
1422             $cache_qtl_plot->set_key($key);
1423
1424             if ($cache_qtl_plot->is_valid)
1425             {
1426             $cache_qtl_plot->destroy();
1427             }
1428             my $key_t = "qtlplot_" . $lg . "_thickbox_" . $pop_id . $term_id;
1429             $cache_qtl_plot->set_key($key_t);
1430
1431             if ($cache_qtl_plot->is_valid)
1432             {
1433             $cache_qtl_plot->destroy();
1434             }
1435             $qtl_image = undef;
1436         }
1437     }
1438
1439     return $qtl_image;
1440
1441 }
1442
1443
1444 =head2 stat_files
1445
1446  Usage: my $stat_param_files = $self->stat_files();
1447  Desc:  creates a master file containing individual files
1448         in /data/prod/tmp/r_qtl for each statistical parameter
1449         which are feed to  R.
1450  Ret:  an absolute path to the statistical parameter's
1451        master file (and individual files)
1452  Args: None
1453  Side Effects:
1454  Example:
1455
1456 =cut
1457
1458 sub stat_files
1459 {
1460     my $self   = shift;
1461     my $pop_id = $self->get_object_id();
1462     my $pop    = $self->get_object();
1463     my $user_id;
1464     if ($c->user)
1465     {
1466         $user_id = $c->user->get_object->get_sp_person_id;
1467     } else
1468     {
1469         $user_id = $pop->get_sp_person_id();
1470     }
1471
1472     my $qtl            = CXGN::Phenome::Qtl->new($user_id);
1473     my $user_stat_file = $qtl->get_stat_file($c, $pop_id);
1474
1475     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1476       $self->cache_temp_path();
1477
1478     open my $user_stat_fh, "<", $user_stat_file or die "can't open file: !$\n";
1479
1480     my $stat_files;
1481
1482     while (<$user_stat_fh>)
1483     {
1484         my ( $parameter, $value ) = split( /\t/, $_ );
1485
1486         my $stat_temp = File::Temp->new(
1487                                       TEMPLATE => "${parameter}_$pop_id-XXXXXX",
1488                                       DIR      => $prod_temp_path,
1489                                       UNLINK   => 0
1490         );
1491         my $stat_file = $stat_temp->filename;
1492
1493         open my $sf_fh, ">", $stat_file or die "can't open file: !$\n";
1494         $sf_fh->print($value);
1495
1496
1497         $stat_files .= $stat_file . "\t";
1498
1499     }
1500
1501     my $stat_param_files =
1502       $prod_temp_path . "/" . "stat_temp_files_pop_id_${pop_id}";
1503
1504     open my $stat_param_fh, ">", $stat_param_files or die "can't open file: !$\n";
1505     $stat_param_fh->print ($stat_files);
1506
1507
1508     return $stat_param_files;
1509
1510 }
1511
1512 =head2 stat_param_hash
1513
1514  Usage: my $stat_param = $self->stat_param_hash();
1515  Desc: creates a hash table (with the statistical parameters (as key) and
1516        their corresponding values) out of a tab delimited
1517        statistical parameters file.
1518  Ret: a hashref for statistical parameter key and value pairs table
1519  Args: None
1520  Side Effects:
1521  Example:
1522
1523 =cut
1524
1525 sub stat_param_hash
1526 {
1527     my $self   = shift;
1528     my $pop_id = $self->get_object_id();
1529     my $pop    = $self->get_object();
1530     my $user_id;
1531     if ($c->user) {
1532         $user_id = $c->user->get_object->get_sp_person_id;
1533     } else {
1534         $user_id = $pop->get_sp_person_id();
1535     }
1536     my $qtl            = CXGN::Phenome::Qtl->new($user_id);
1537     my $user_stat_file = $qtl->get_stat_file($c, $pop_id);
1538
1539     open my $user_stat_fh, "<", $user_stat_file or die "can't open file: !$\n";
1540
1541     my %stat_param;
1542
1543     while (<$user_stat_fh>)
1544     {
1545         my ( $parameter, $value ) = split( /\t/, $_ );
1546
1547         $stat_param{$parameter} = $value;
1548
1549     }
1550
1551     return \%stat_param;
1552 }
1553
1554 sub submitter
1555 {
1556     my $self         = shift;
1557     my $population   = $self->get_object();
1558     my $sp_person_id = $population->get_sp_person_id();
1559     my $submitter    = CXGN::People::Person->new( $self->get_dbh(),
1560                                               $population->get_sp_person_id() );
1561     my $submitter_name =
1562       $submitter->get_first_name() . " " . $submitter->get_last_name();
1563     my $submitter_link =
1564 qq |<a href="/solpeople/personal-info.pl?sp_person_id=$sp_person_id">$submitter_name</a> |;
1565
1566     return $submitter, $submitter_link;
1567
1568 }
1569
1570 #move to qtl or population object
1571 sub legend {
1572     my $self = shift;
1573     my $pop  = $self->get_object();
1574     my $user_id;
1575     if ($c->user) {
1576         $user_id = $c->user->get_object->get_sp_person_id;
1577     } else {
1578         $user_id = $pop->get_sp_person_id();
1579     }
1580
1581     my $qtl       = CXGN::Phenome::Qtl->new($user_id);
1582     my $stat_file = $qtl->get_stat_file($c, $pop->get_population_id());
1583     my @stat;
1584     my $ci=1;
1585
1586     open my $sf, "<", $stat_file or die "$! reading $stat_file\n";
1587     while (my $row = <$sf>)
1588     {
1589         chomp($row);
1590         my ( $parameter, $value ) = split( /\t/, $row );
1591         if ($parameter =~/qtl_method/) {$parameter = 'Mapping method';}
1592         if ($parameter =~/qtl_model/) {$parameter = 'Mapping model';}
1593         if ($parameter =~/prob_method/) {$parameter = 'QTL genotype probability method';}
1594         if ($parameter =~/step_size/) {$parameter = 'Genome scan size (cM)';}
1595         if ($parameter =~/permu_level/) {$parameter = 'Permutation significance level';}
1596         if ($parameter =~/permu_test/) {$parameter = 'No. of permutations';}
1597         if ($parameter =~/prob_level/) {$parameter = 'QTL genotype signifance level';}
1598         if ($parameter =~/stat_no_draws/) {$parameter = 'No. of imputations';}
1599
1600         if ( $value eq 'zero' || $value eq 'Marker Regression' )
1601         {
1602             $ci = 0;
1603         }
1604
1605         unless (($parameter =~/No. of imputations/ && !$value ) ||
1606                 ($parameter =~/QTL genotype probability/ && !$value ) ||
1607                 ($parameter =~/Permutation significance level/ && !$value)
1608                )
1609
1610         {
1611             if ($parameter =~/Genome scan/ && $value eq 'zero' || !$value)
1612             {
1613                 $value = '0.0'
1614             }
1615
1616             push @stat, [map{$_} ($parameter, $value)];
1617
1618         }
1619
1620     }
1621
1622     my $sm;
1623     foreach my $st (@stat) {
1624         foreach my $i (@$st) {
1625             if ($i =~/zero/)  {
1626                 foreach my $s (@stat) {
1627                     foreach my $j (@$s) {
1628                         $j =~ s/Maximum Likelihood/Marker Regression/;
1629                         $ci = 0;
1630                     }
1631                 }
1632             }
1633         }
1634     }
1635
1636
1637     my $permu_threshold_ref = $self->permu_values();
1638     my %permu_threshold     = %$permu_threshold_ref;
1639
1640     my @keys;
1641
1642     foreach my $key ( keys %permu_threshold )
1643     {
1644         if ( $key =~ m/^\d./ )
1645         {
1646             push @keys, $key;
1647         }
1648
1649     }
1650     my $lod1 = $permu_threshold{ $keys[0] };
1651
1652     if  (!$lod1)
1653     {
1654         $lod1 = qq |<i>Not calculated</i>|;
1655     }
1656
1657     push @stat,
1658     [
1659      map {$_} ('LOD threshold', $lod1)
1660     ];
1661
1662
1663     if ($ci) {
1664         push @stat,
1665         [
1666          map {$_} ('Confidence interval', 'Based on 95% Bayesian Credible Interval')
1667         ];
1668     }
1669
1670     push @stat,
1671     [
1672      map {$_} ('QTL software', "<a href=http://www.rqtl.org>R/QTL</a>")
1673     ];
1674     push @stat,
1675     [
1676      map {$_} ('Citation', "<a href=http://www.biomedcentral.com/1471-2105/11/525>solQTL</a>")
1677     ];
1678
1679     my $legend_data = columnar_table_html (
1680                                            headings    => [
1681                                                           ' ',
1682                                                           ' ',
1683                                                           ],
1684                                            data        => \@stat,
1685                                           __alt_freq   => 2,
1686                                           __alt_width  => 1,
1687                                           __alt_offset => 3,
1688                                           __align      => 'l',
1689
1690                                           );
1691
1692
1693     return $legend_data;
1694
1695 }
1696
1697 =head2 set_trait_id, get_trait_id
1698
1699  Usage:
1700  Desc: the 'cvterm id' here is not necessarily a cvterm id,
1701        but may also be user submitted trait id...
1702  Ret:
1703  Args:
1704  Side Effects:
1705  Example:
1706
1707 =cut
1708
1709
1710
1711 sub get_trait_id {
1712     my $self = shift;
1713     return $self->{cvterm_id};
1714 }
1715 sub set_trait_id {
1716     my $self = shift;
1717     return $self->{cvterm_id} = shift;
1718 }
1719
1720 =head2 get_trait_name
1721  Usage: my $term_name = $self->get_trait_name()
1722  Desc: retrieves the name of the trait whether
1723        it is stored in the user_trait or cvterm table
1724  Return: a trait name
1725  Args: None
1726  Side Effects:
1727  Example:
1728
1729 =cut
1730 sub get_trait_name {
1731     my $self        = shift;
1732     my $population  = $self->get_object();
1733     my $term_id     = $self->get_trait_id();
1734     my $dbh         = $c->dbc()->dbh();
1735
1736     my ($term_obj, $term_name);
1737     if ( $population->get_web_uploaded() )
1738     {
1739         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $term_id );
1740         $term_name = $term_obj->get_name();
1741     }
1742     else
1743     {
1744         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $term_id );
1745         $term_name = $term_obj->get_cvterm_name();
1746     }
1747
1748     return $term_name;
1749
1750
1751 }
1752
1753 sub qtl_effects {
1754     my $self       = shift;
1755     my $trait_name = $self->get_trait_name();
1756     my $pop        = $self->get_object();
1757     $trait_name   = $pop->cvterm_acronym($trait_name);
1758
1759     my $file = $pop->qtl_effects_file($c, $trait_name);
1760
1761     if ( -s $file > 1 )
1762     {
1763
1764         my @effects =  map  { [ split( /\t/, $_) ]}  read_file( $file );
1765         my $trash   = shift(@effects);
1766
1767         push @effects, map { [ $_ ] } (" ", "QTL effects interpretation example: 2\@100 means
1768                                         QTL at linkage group 2 and position 100cM.",
1769                                         "2\@100:3\@100 means interaction between QTL at linkage
1770                                         group 2 position 100 cM and QTL at linkage group 3
1771                                         position 100 cM. 'a' and 'd' stand for additive and domininace
1772                                         effects, respectively."
1773                                       );
1774         return \@effects;
1775     } else
1776     {
1777         return undef;
1778     }
1779
1780 }
1781
1782 sub explained_variation {
1783     my $self       = shift;
1784     my $trait_name = $self->get_trait_name();
1785     my $pop        = $self->get_object();
1786     $trait_name    = $pop->cvterm_acronym($trait_name);
1787
1788     my $file = $pop->explained_variation_file($c, $trait_name);
1789
1790     if ( -s $file > 1 )
1791     {
1792         my @anova =  map  { [ split( /\t/, $_) ]}  read_file( $file );
1793         $anova[0][0] = "Source";
1794
1795         if ( $anova[1][0] eq 'Model')
1796         {
1797             push @anova, map { [ $_ ] } ("  ", "The ANOVA model is based on a single QTL
1798                                          significant source of variation"
1799                                         );
1800         }
1801         else
1802         {
1803             push @anova, map { [ $_ ] } ( "  ",  "Variance source interpretation example: 2\@100 means
1804                                           QTL at linkage group 2 and position 100cM.",
1805                                           "2\@100:3\@100 means interaction between QTL at linkage
1806                                           group 2 position
1807                                           100 cM and QTL at linkage group 3 position 100 cM."
1808                                        );
1809         }
1810
1811
1812         return \@anova;
1813     } else
1814     {
1815         return undef;
1816     }
1817
1818 }
1819
1820 sub qtl_stat_option {
1821     my $self    = shift;
1822     my $pop_id  = $self->get_object_id();
1823     my $user_id = $c->user->get_object->get_sp_person_id if $c->user;
1824     my $qtl_obj = CXGN::Phenome::Qtl->new($user_id);
1825
1826     my ($user_qtl_dir, $user_dir) = $qtl_obj->get_user_qtl_dir($c);
1827     my $stat_options_file         = "$user_dir/stat_options_pop_${pop_id}.txt";
1828
1829     my $stat_option = -e $stat_options_file && read_file($stat_options_file) =~ /Yes/ ? 'default'
1830                     : !-e $stat_options_file                                          ? 'default'
1831                     :                                                                   'user params'
1832                     ;
1833
1834     return $stat_option;
1835 }
1836
1837 sub remove_permu_file {
1838     my $self = shift;
1839     my $population = $self->get_object();
1840
1841     if ($self->qtl_stat_option eq 'user params')
1842     {
1843         my $ac = $population->cvterm_acronym($self->get_trait_name());
1844         my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1845             $self->cache_temp_path();
1846
1847         my $file_cache = Cache::File->new( cache_root => $prod_cache_path );
1848         my $key_permu = $ac . "_" . $population->get_population_id() . "_permu";
1849
1850         unlink($self->permu_file());
1851         $file_cache->remove($key_permu);
1852
1853     }
1854 }