cgi-bin/phenome/population_indls.pl

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