cgi-bin/phenome/population_indls.pl

   1 #!/usr/bin/perl -w
   2
   3 =head1 DESCRIPTION
   4
   5 Creates a trait/cvterm page with a description of
   6 the population on which the trait/cvterm was evaluated,
   7 displays the frequency distribution of its phenotypic data
   8 and most importantly produces the on-the-fly  QTL analysis
   9 output for the trait and more....
  10
  11 =head1 AUTHOR
  12
  13 Isaak Y Tecle (iyt2@cornell.edu)
  14
  15 =cut
  16
  17 use strict;
  18 use warnings;
  19
  20 our $c;
  21 my $population_indls_detail_page =
  22   CXGN::Phenome::PopulationIndlsDetailPage->new();
  23
  24 package CXGN::Phenome::PopulationIndlsDetailPage;
  25
  26 use CXGN::Page;
  27 use CXGN::Page::FormattingHelpers qw /info_section_html
  28   page_title_html
  29   columnar_table_html
  30   html_optional_show
  31   info_table_html
  32   tooltipped_text
  33   html_alternate_show
  34   /;
  35
  36 use CXGN::Phenome::Population;
  37 use CXGN::Phenome::Qtl;
  38 use CXGN::Phenome::PopulationDbxref;
  39 use CXGN::Tools::WebImageCache;
  40 use SGN::Context;
  41 use CXGN::People::PageComment;
  42 use CXGN::People::Person;
  43 use CXGN::Chado::Publication;
  44 use CXGN::Chado::Pubauthor;
  45
  46 use GD;
  47 use GD::Graph::bars;
  48 use GD::Graph::lines;
  49 use GD::Graph::points;
  50 use GD::Graph::Map;
  51 use Statistics::Descriptive;
  52 use Math::Round::Var;
  53 use File::Temp qw /tempfile tempdir/;
  54 use File::Copy;
  55 use File::Spec;
  56 use File::Basename;
  57 use File::stat;
  58 use Cache::File;
  59 use CXGN::Scrap::AjaxPage;
  60 use CXGN::Contact;
  61 use Storable qw / store /;
  62
  63
  64 use CXGN::Page::UserPrefs;
  65 use base qw / CXGN::Page::Form::SimpleFormPage CXGN::Phenome::Main/;
  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,                                #<  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           = $self->qtl_plot();
 425
 426         my $legend = $self->legend($population);
 427         my $qtl_html = qq | <table><tr><td width=70%>$qtl_image</td><td width=30%>$legend</td></tr></table> |;
 428
 429         print info_section_html(
 430                                 title    => 'QTL(s)',
 431                                 contents => $qtl_html,
 432                                 );
 433
 434         print info_section_html(
 435                                 title    => 'Phenotype Frequency Distribution',
 436                                 contents => $plot_html,
 437                                );
 438
 439         print info_section_html(
 440                                 title    => 'Phenotype Data',
 441                                 contents => $data_view . " " . $data_download,
 442                                );
 443
 444     }
 445     else
 446     {
 447         my $message =
 448           "The QTL data for this trait in this population is not public yet.
 449                        If you would like to know more about this data,
 450                        please contact the owner of the data: <b>$submitter_link</b>
 451                        or email to SGN:
 452                        <a href=mailto:sgn-feedback\@sgn.cornell.edu>
 453                        sgn-feedback\@sgn.cornell.edu</a>.\n";
 454
 455         print info_section_html( title    => 'QTL(s)',
 456                                  contents => $message );
 457     }
 458
 459     print info_section_html(
 460                            title => 'Literature Annotation',
 461                            contents => $pubmed,
 462                            );
 463
 464     if ($population_name)
 465     {
 466         my $page_comment_obj =
 467             CXGN::People::PageComment->new( $self->get_dbh(), "population",
 468                                             $population_id,
 469                                             $self->get_page()->{request}->uri()."?".$self->get_page()->{request}->args()
 470             );
 471         print $page_comment_obj->get_html();
 472     }
 473
 474     $self->get_page()->footer();
 475
 476     exit();
 477 }
 478
 479 # override store to check if a locus with the submitted symbol/name already exists in the database
 480
 481 sub store
 482 {
 483     my $self          = shift;
 484     my $population    = $self->get_object();
 485     my $population_id = $self->get_object_id();
 486     my %args          = $self->get_args();
 487
 488     $self->SUPER::store(0);
 489
 490     exit();
 491 }
 492
 493 sub population_distribution
 494 {
 495     my $self = shift;
 496     my $doc  = CXGN::Scrap::AjaxPage->new();
 497
 498     my ( $pop_id, $cvterm_id ) =
 499       $doc->get_encoded_arguments( "population_id", "cvterm_id" );
 500
 501     my $dbh = CXGN::DB::Connection->new();
 502
 503     my ( $term_obj, $term_name, $term_id );
 504
 505     my $pop = CXGN::Phenome::Population->new( $dbh, $pop_id );
 506
 507     if ( $pop->get_web_uploaded() )
 508     {
 509         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $cvterm_id );
 510         $term_name = $term_obj->get_name();
 511         $term_id   = $term_obj->get_user_trait_id();
 512     }
 513     else
 514     {
 515         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $cvterm_id );
 516         $term_name = $term_obj->get_cvterm_name();
 517         $term_id   = $term_obj->get_cvterm_id();
 518     }
 519
 520     my $vh           = SGN::Context->new();
 521     my $basepath     = $vh->get_conf("basepath");
 522     my $tempfile_dir = $vh->get_conf("tempfiles_subdir");
 523
 524     my $cache = CXGN::Tools::WebImageCache->new();
 525     $cache->set_basedir($basepath);
 526     $cache->set_temp_dir( $tempfile_dir . "/temp_images" );
 527     $cache->set_expiration_time(259200);
 528     $cache->set_key( "popluation_distribution" . $pop_id . $term_id );
 529     $cache->set_map_name("popmap$pop_id$term_id");
 530
 531     my $pop_name;
 532     my ( $variance, $std,     $mean );
 533     my ( @value,    @indl_id, @indl_name );
 534
 535     $cache->set_force(0);
 536     if ( !$cache->is_valid() )
 537     {
 538         my $pop_obj = CXGN::Phenome::Population->new( $dbh, $pop_id );
 539         $pop_name = $pop_obj->get_name();
 540         my ( $indl_id, $indl_name, $value ) = $pop_obj->plot_cvterm($term_id);
 541         my @indl_id   = @{$indl_id};
 542         my @indl_name = @{$indl_name};
 543         @value = @{$value};
 544
 545         my $round = Math::Round::Var->new(0.001);
 546
 547         my $stat = Statistics::Descriptive::Full->new();
 548
 549         $stat->add_data(@value);
 550
 551         my $stat_para = Statistics::Descriptive::Sparse->new();
 552         $stat_para->add_data(@value);
 553         $std  = $stat_para->standard_deviation();
 554         $mean = $stat_para->mean();
 555
 556         my %f = $stat->frequency_distribution(10);
 557
 558         my ( @keys, @counts );
 559
 560         for ( sort { $a <=> $b } keys %f )
 561         {
 562             my $key = $round->round($_);
 563             push @keys,   $key;
 564             push @counts, $f{$_};
 565         }
 566
 567         my $min = $stat->min();
 568         if ( $min != 0 )
 569         {
 570             $min = $min - 0.01;
 571         }
 572
 573         my @keys_range = $min . '-' . $keys[0];
 574
 575         my $range;
 576         my $previous_k   = $keys[0];
 577         my $keys_shifted = shift(@keys);
 578         foreach my $k (@keys)
 579         {
 580             $range = $previous_k . '-' . $k;
 581             push @keys_range, $range;
 582             $previous_k = $k;
 583         }
 584
 585         my $max = $counts[0];
 586         foreach my $i ( @counts[ 1 .. $#counts ] )
 587         {
 588             if ( $i > $max ) { $max = $i; }
 589         }
 590         $max = int( $max + ( $max * 0.1 ) );
 591
 592         my $c_html;
 593         my @c_html;
 594         my ( $lower, $upper );
 595
 596         foreach my $k (@keys_range)
 597         {
 598             ( $lower, $upper ) = split( /-/, $k );
 599             $c_html =
 600 qq | /phenome/indls_range_cvterm.pl?cvterm_id=$term_id&amp;lower=$lower&amp;upper=$upper&amp;population_id=$pop_id |;
 601             push @c_html, $c_html;
 602
 603         }
 604
 605         my @bar_clr = ("orange");
 606         my @data    = ( [@keys_range], [@counts] );
 607         my $graph   = new GD::Graph::bars();
 608
 609         $graph->set_title_font('gdTinyFont');
 610         $graph->set(
 611                      title             => " ",
 612                      x_label           => "Ranges for $term_name",
 613                      y_label           => "Frequency",
 614                      y_max_value       => $max,
 615                      x_all_ticks       => 1,
 616                      y_all_ticks       => 2,
 617                      y_label_skip      => 5,
 618                      y_plot_values     => 0,
 619                      x_label_skip      => 1,
 620                      width             => 400,
 621                      height            => 400,
 622                      bar_width         => 30,
 623                      x_labels_vertical => 1,
 624                      show_values       => 1,
 625                      textclr           => "black",
 626                      dclrs             => \@bar_clr,
 627                    );
 628
 629         $cache->set_image_data( $graph->plot( \@data )->png );
 630
 631         my $map = new GD::Graph::Map(
 632                                       $graph,
 633                                       hrefs       => [ \@c_html ],
 634                                       noImgMarkup => 1,
 635                                       mapName     => "popmap$pop_id$term_id",
 636                                       info        => "%x: %y lines",
 637                                     );
 638         $cache->set_image_map_data(
 639                       $map->imagemap( "popimage$pop_id$term_id.png", \@data ) );
 640
 641     }
 642
 643     my $image_map = $cache->get_image_map_data();
 644     my $image     = $cache->get_image_tag();
 645     my $title =
 646 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. |;
 647
 648     return $image, $title, $image_map;
 649 }
 650
 651 sub qtl_plot
 652 {
 653
 654     my $self = shift;
 655     my $doc  = CXGN::Scrap::AjaxPage->new();
 656
 657     my ( $pop_id, $cvterm_id ) =
 658       $doc->get_encoded_arguments( "population_id", "cvterm_id" );
 659
 660     my $dbh = $self->get_dbh();
 661
 662     my $population     = $self->get_object();
 663     my $pop_name       = $population->get_name();
 664     my $mapversion     = $population->mapversion_id();
 665     my @linkage_groups = $population->linkage_groups();
 666
 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 $vh           = SGN::Context->new();
 686     my $basepath     = $vh->get_conf("basepath");
 687     my $tempfile_dir = $vh->get_conf("tempfiles_subdir");
 688
 689     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
 690       $self->cache_temp_path();
 691     my $cache_tempimages = Cache::File->new( cache_root => $tempimages_path );
 692     $cache_tempimages->purge();
 693
 694     my ( @marker,  @chr,  @pos,   @lod );
 695     my ( @chr_qtl, @left, @right, @peak );
 696     my ( $qtl_image, $image, $image_t, $image_url, $image_html, $image_t_url,
 697          $thickbox, $title, $l_m, $p_m, $r_m );
 698
 699     my $round1 = Math::Round::Var->new(0.1);
 700     my $round2 = Math::Round::Var->new(1);
 701
 702     $qtl_image  = $self->qtl_images_exist();
 703     my $permu_data = $self->permu_values_exist();
 704
 705     unless ( $qtl_image && $permu_data )
 706     {
 707
 708         my ( $qtl_summary, $flanking_markers ) = $self->run_r();
 709
 710         open QTLSUMMARY, "<$qtl_summary" or die "can't open $qtl_summary: $!\n";
 711
 712         my $header = <QTLSUMMARY>;
 713         while ( my $row = <QTLSUMMARY> )
 714         {
 715             my ( $marker, $chr, $pos, $lod ) = split( /\t/, $row );
 716             push @marker, $marker;
 717             push @chr,    $chr;
 718             $pos = $round2->round($pos);
 719             push @pos, $pos;
 720             $lod = $round1->round($lod);
 721             push @lod, $lod;
 722         }
 723
 724         my @o_lod = sort(@lod);
 725         my $max   = $o_lod[-1];
 726         $max = $max + (0.5);
 727
 728         close QTLSUMMARY;
 729
 730         open MARKERS, "<$flanking_markers"
 731           or die "can't open $flanking_markers: !$\n";
 732
 733         $header = <MARKERS>;
 734         while ( my $row = <MARKERS> )
 735
 736         {
 737             chomp($row);
 738             my ($trash, $chr_qtl, $left, $peak, $right, $peakmarker ) = split( /\t/, $row );
 739             push @chr_qtl, $chr_qtl;
 740             push @left,    $left;
 741             push @right,   $right;
 742             push @peak, $peakmarker;
 743         }
 744
 745         close MARKERS;
 746         my (@h_markers, @chromosomes, @lk_groups);
 747         my $h_marker;
 748
 749
 750         @lk_groups = @linkage_groups;
 751         @lk_groups = sort ( { $a <=> $b } @lk_groups );
 752         for ( my $i = 0 ; $i < @left ; $i++ )
 753         {
 754             my $lg           = shift(@lk_groups);
 755             my $key_h_marker = "$ac" . "_pop_" . "$pop_id" . "_chr_" . $lg;
 756             $h_marker = $cache_tempimages->get($key_h_marker);
 757
 758             unless ($h_marker)
 759             {
 760
 761                 push @chromosomes, $lg;
 762                 $l_m = $left[$i];
 763                 $r_m = $right[$i];
 764                 $p_m = $peak[$i];
 765                 s/\s//g for $l_m, $r_m, $p_m;
 766
 767                 my $l_pos =
 768                     $population->get_marker_position( $mapversion, $l_m );
 769                 my $r_pos =
 770                     $population->get_marker_position( $mapversion, $r_m );
 771
 772
 773
 774                 my $permu_threshold_ref = $self->permu_values();
 775                 my %permu_threshold     = %$permu_threshold_ref;
 776                 my @p_keys;
 777                 foreach my $key ( keys %permu_threshold )
 778                 {
 779                     if ( $key =~ m/^\d./ )
 780                     {
 781                         push @p_keys, $key;
 782                     }
 783
 784                 }
 785                 my $lod1 = $permu_threshold{ $p_keys[0] };
 786                 $h_marker =
 787                     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|;
 788
 789
 790                 $cache_tempimages->set( $key_h_marker, $h_marker, '30 days' );
 791             }
 792
 793             push @h_markers, $h_marker;
 794         }
 795
 796         my $count       = 0;
 797         my $old_chr_chr = 1;
 798         my (
 799              $chr_chr, $image,     $image_t,
 800              $thickbox,  $max_chr, $chr_chr_e, $marker_chr_e,
 801              $pos_chr_e, $lod_chr_e
 802            );
 803         my $chrs = ( scalar(@chromosomes) ) + 1;
 804
 805         for ( my $i = 1 ; $i < $chrs ; $i++ )
 806         {
 807             my ( @marker_chr, @chr_chr, @pos_chr, @lod_chr, @data, @m_html ) =
 808               ();
 809             my ( $marker_chr, $pos_chr, $lod_chr, $max_chr );
 810
 811             $h_marker = shift(@h_markers);
 812
 813             if ( ( $i == $old_chr_chr ) && ( $i != 12 ) )
 814             {
 815                 push @marker_chr, $marker_chr_e;
 816                 push @chr_chr,    $chr_chr_e;
 817                 $pos_chr_e = $round2->round($pos_chr_e);
 818                 push @pos_chr, $pos_chr_e;
 819                 $lod_chr = $round1->round($lod_chr_e);
 820                 push @lod_chr, $lod_chr_e;
 821             }
 822
 823             my $cache_qtl_plot = CXGN::Tools::WebImageCache->new();
 824             $cache_qtl_plot->set_basedir($basepath);
 825             $cache_qtl_plot->set_temp_dir( $tempfile_dir . "/temp_images" );
 826             $cache_qtl_plot->set_expiration_time(259200);
 827             $cache_qtl_plot->set_key(
 828                                 "qtlplot" . $i . "small" . $pop_id . $term_id );
 829             $cache_qtl_plot->set_force(0);
 830
 831             if ( !$cache_qtl_plot->is_valid() )
 832             {
 833
 834                 for ( my $j = 0 ; $j < @marker ; $j++ )
 835                 {
 836
 837                     $chr_chr = $chr[$j];
 838
 839                     if ( $i == $chr_chr )
 840                     {
 841                         $marker_chr = $marker[$j];
 842
 843                         $pos_chr = $pos[$j];
 844                         $lod_chr = $lod[$j];
 845
 846                         push @marker_chr, $marker_chr;
 847                         push @chr_chr,    $chr_chr;
 848                         $pos_chr = $round2->round($pos_chr);
 849                         push @pos_chr, $pos_chr;
 850                         $lod_chr = $round1->round($lod_chr);
 851                         push @lod_chr, $lod_chr;
 852
 853                         ( $chr_chr_e, $marker_chr_e, $pos_chr_e, $lod_chr_e ) =
 854                           ();
 855                     }
 856
 857                     elsif ( $i != $chr_chr )
 858                     {
 859
 860                         $chr_chr_e    = $chr[$j];
 861                         $marker_chr_e = $marker[$j];
 862                         $pos_chr_e    = $pos[$j];
 863                         $lod_chr_e    = $lod[$j];
 864                     }
 865
 866                 }
 867
 868                 @data = ( [ (@pos_chr) ], [@lod_chr] );
 869                 my $graph = new GD::Graph::lines( 110, 110 );
 870                 $graph->set_title_font('gdTinyFont');
 871                 $graph->set(
 872                              title             => " ",
 873                              x_label           => "Chr $i (cM)",
 874                              y_label           => "LOD",
 875                              y_max_value       => 10,
 876                              x_all_ticks       => 5,
 877                              y_all_ticks       => 1,
 878                              y_label_skip      => 1,
 879                              y_plot_values     => 1,
 880                              x_label_skip      => 5,
 881                              x_plot_values     => 1,
 882                              x_labels_vertical => 1,
 883                              textclr           => "black"
 884                            );
 885
 886                 $cache_qtl_plot->set_image_data( $graph->plot( \@data )->png );
 887
 888             }
 889
 890             $image      = $cache_qtl_plot->get_image_tag();
 891             $image_url  = $cache_qtl_plot->get_image_url();
 892
 893
 894 ###########thickbox
 895             my $cache_qtl_plot_t = CXGN::Tools::WebImageCache->new();
 896             $cache_qtl_plot_t->set_basedir($basepath);
 897             $cache_qtl_plot_t->set_temp_dir( $tempfile_dir . "/temp_images" );
 898             $cache_qtl_plot_t->set_expiration_time(259200);
 899             $cache_qtl_plot_t->set_key(
 900                           "qtlplot_" . $i . "_thickbox_" . $pop_id . $term_id );
 901             $cache_qtl_plot_t->set_force(0);
 902
 903             if ( !$cache_qtl_plot_t->is_valid() )
 904             {
 905                 my @o_lod_chr = sort { $a <=> $b } @lod_chr;
 906                 $max_chr = pop(@o_lod_chr);
 907                 $max_chr = $max_chr + (0.5);
 908
 909                 my $graph_t = new GD::Graph::lines( 420, 420 );
 910                 $graph_t->set_title_font('gdTinyFont');
 911                 $graph_t->set(
 912                                title             => " ",
 913                                x_label           => "Chromosome $i (cM)",
 914                                y_label           => "LOD",
 915                                y_max_value       => $max_chr,
 916                                x_all_ticks       => 5,
 917                                y_all_ticks       => 1,
 918                                y_label_skip      => 1,
 919                                y_plot_values     => 1,
 920                                x_label_skip      => 5,
 921                                x_plot_values     => 1,
 922                                x_labels_vertical => 1,
 923                                textclr           => "black"
 924                              );
 925
 926                 $cache_qtl_plot_t->set_image_data(
 927                                                 $graph_t->plot( \@data )->png );
 928
 929             }
 930
 931             $image_t     = $cache_qtl_plot_t->get_image_tag();
 932             $image_t_url = $cache_qtl_plot_t->get_image_url();
 933
 934             $thickbox =
 935 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> |;
 936
 937             $qtl_image .= $thickbox;
 938             $title       = "  ";
 939             $old_chr_chr = $chr_chr;
 940         }
 941     }
 942
 943     return $qtl_image;
 944 }
 945
 946 =head2 infile_list
 947  Usage: my $file_in = $self->infile_list();
 948  Desc: returns an R input tempfile containing a tempfile
 949        holding the cvterm acronym, pop id, a filepath to the phenotype dataset file,
 950         a filepath to genotype dataset file, a filepath to the permuation file.
 951  Ret: an R input tempfile name (with abosulte path)
 952  Args:
 953  Side Effects:
 954  Example:
 955
 956
 957 =cut
 958
 959 sub infile_list
 960 {
 961
 962     my $self = shift;
 963     my $doc  = CXGN::Scrap::AjaxPage->new();
 964
 965     my ( $pop_id, $cvterm_id ) =
 966       $doc->get_encoded_arguments( "population_id", "cvterm_id" );
 967
 968     my $dbh = $self->get_dbh();
 969
 970     my ( $term_obj, $term_name, $term_id );
 971     my $population = $self->get_object();
 972
 973     if ( $population->get_web_uploaded() )
 974     {
 975         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $cvterm_id );
 976         $term_name = $term_obj->get_name();
 977         $term_id   = $term_obj->get_user_trait_id();
 978     }
 979     else
 980     {
 981         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $cvterm_id );
 982         $term_name = $term_obj->get_cvterm_name();
 983         $term_id   = $term_obj->get_cvterm_id();
 984     }
 985
 986     my $ac = $population->cvterm_acronym($term_name);
 987
 988     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
 989       $self->cache_temp_path();
 990
 991     my $prod_permu_file  = $self->permu_file();
 992     my $gen_dataset_file = $self->genotype_file();
 993     my $phe_dataset_file = $self->phenotype_file();
 994     my $crosstype_file   = $self->crosstype_file();
 995
 996     my $input_file_list_temp =
 997       File::Temp->new(
 998                        TEMPLATE => "infile_list_${ac}_$pop_id-XXXXXX",
 999                        DIR      => $prod_temp_path,
1000                        UNLINK   => 0,
1001                      );
1002     my $file_in = $input_file_list_temp->filename();
1003
1004     my $file_cvin = File::Temp->new(
1005                                      TEMPLATE => 'cv_input-XXXXXX',
1006                                      DIR      => $prod_temp_path,
1007                                      UNLINK   => 0,
1008                                    );
1009     my $file_cv_in = $file_cvin->filename();
1010
1011     open CV, ">$file_cv_in" or die "can't open $file_cv_in: $!\n";
1012     print CV $ac;
1013     close CV;
1014
1015     my $file_in_list = join( "\t",
1016                              $file_cv_in,       "P$pop_id",
1017                              $gen_dataset_file, $phe_dataset_file,
1018                              $prod_permu_file, $crosstype_file);
1019
1020     open FI, ">$file_in" or die "can't open $file_in: $!\n";
1021     print FI $file_in_list;
1022     close FI;
1023
1024     return $file_in;
1025
1026 }
1027
1028 =head2 outfile_list
1029
1030  Usage: my ($file_out, $qtl_summary, $flanking_markers) = $self->outfile_list();
1031  Desc: returns an R output tempfile containing a tempfile supposed to hold the qtl
1032        mapping output and another tempfile for the qtl flanking markers
1033        and the qtl mapping output and qtl flanking markers files separately
1034        (convenient for reading their data when plotting the qtl)
1035  Ret: R output file names (with abosulte path)
1036  Args:
1037  Side Effects:
1038  Example:
1039
1040 =cut
1041
1042 sub outfile_list
1043 {
1044     my $self = shift;
1045
1046     my $doc = CXGN::Scrap::AjaxPage->new();
1047
1048     my ( $pop_id, $cvterm_id ) =
1049       $doc->get_encoded_arguments( "population_id", "cvterm_id" );
1050
1051     my $dbh = $self->get_dbh();
1052
1053     my ( $term_obj, $term_name, $term_id );
1054     my $population = $self->get_object();
1055
1056     if ( $population->get_web_uploaded() )
1057     {
1058         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $cvterm_id );
1059         $term_name = $term_obj->get_name();
1060         $term_id   = $term_obj->get_user_trait_id();
1061     }
1062     else
1063     {
1064         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $cvterm_id );
1065         $term_name = $term_obj->get_cvterm_name();
1066         $term_id   = $term_obj->get_cvterm_id();
1067     }
1068
1069     my $ac = $population->cvterm_acronym($term_name);
1070
1071     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1072       $self->cache_temp_path();
1073
1074     my $output_file_list_temp =
1075       File::Temp->new(
1076                        TEMPLATE => "outfile_list_${ac}_$pop_id-XXXXXX",
1077                        DIR      => $prod_temp_path,
1078                        UNLINK   => 0,
1079                      );
1080     my $file_out = $output_file_list_temp->filename();
1081
1082     my $qtl_temp = File::Temp->new(
1083                                  TEMPLATE => "qtl_summary_${ac}_$pop_id-XXXXXX",
1084                                  DIR      => $prod_temp_path,
1085                                  UNLINK   => 0
1086     );
1087     my $qtl_summary = $qtl_temp->filename;
1088
1089     my $marker_temp = File::Temp->new(
1090                             TEMPLATE => "flanking_markers_${ac}_$pop_id-XXXXXX",
1091                             DIR      => $prod_temp_path,
1092                             UNLINK   => 0
1093     );
1094
1095     my $flanking_markers = $marker_temp->filename;
1096
1097     my $file_out_list = join(
1098         "\t",
1099         $qtl_summary,
1100         $flanking_markers,
1101                             );
1102     open FO, ">$file_out" or die "can't open $file_out: $!\n";
1103     print FO $file_out_list;
1104     close FO;
1105
1106     return $file_out, $qtl_summary, $flanking_markers;
1107 }
1108
1109 =head2 cache_temp_path
1110
1111  Usage: my ($prod_cache_path, $prod_temp_path, $tempimages_path) = $self->cache_temp_path();
1112  Desc: creates the 'r_qtl' dir in the '/data/prod/tmp/' dir;
1113       'cache' and 'tempfiles' in the /data/prod/tmp/r_qtl/,
1114        and 'temp_images' in the /data/local/cxgn/sgn/documents/tempfiles'
1115  Ret: /data/prod/tmp/r_qtl/cache, /data/prod/tmp/r_qtl/tempfiles,
1116       /data/local/cxgn/sgn/documents/tempfiles/temp_images
1117  Args: none
1118  Side Effects:
1119  Example:
1120
1121 =cut
1122
1123 sub cache_temp_path
1124 {
1125     my $vh           = SGN::Context->new();
1126     my $basepath     = $vh->get_conf("basepath");
1127     my $tempfile_dir = $vh->get_conf("tempfiles_subdir");
1128
1129     my $tempimages_path =
1130       File::Spec->catfile( $basepath, $tempfile_dir, "temp_images" );
1131
1132     my $prod_temp_path = $vh->get_conf('r_qtl_temp_path');
1133     mkdir $prod_temp_path;
1134     my $prod_cache_path = "$prod_temp_path/cache";
1135     mkdir $prod_cache_path;
1136     $prod_temp_path = "$prod_temp_path/tempfiles";
1137     mkdir $prod_temp_path;
1138     -d $prod_temp_path
1139       or die "temp dir '$prod_temp_path' not found, and could not create!";
1140     -r $prod_temp_path or die "temp dir '$prod_temp_path' not readable!";
1141     -w $prod_temp_path or die "temp dir '$prod_temp_path' not writable!";
1142
1143     return $prod_cache_path, $prod_temp_path, $tempimages_path;
1144
1145 }
1146
1147 =head2 genotype_file
1148
1149  Usage: my $gen_file = $self->genotype_file();
1150  Desc: creates the genotype file in the /data/prod/tmp/r_qtl/cache,
1151        if it does not exist yet and caches it for R.
1152  Ret: genotype filename (with abosolute path)
1153  Args: none
1154  Side Effects:
1155  Example:
1156
1157 =cut
1158
1159 sub genotype_file
1160 {
1161     my $self       = shift;
1162     my $pop_id     = $self->get_object_id();
1163     my $population = $self->get_object();
1164
1165     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1166       $self->cache_temp_path();
1167     my $file_cache = Cache::File->new( cache_root => $prod_cache_path );
1168     $file_cache->purge();
1169
1170     my $key_gen          = "popid_" . $pop_id . "_genodata";
1171     my $gen_dataset_file = $file_cache->get($key_gen);
1172
1173     unless ($gen_dataset_file)
1174     {
1175         my $genodata     = $population->genotype_dataset();
1176         my $geno_dataset = ${$genodata};
1177
1178         my $filename = "genodata_" . $pop_id . ".csv";
1179         my $file     = "$prod_cache_path/$filename";
1180
1181         open OUT, ">$file" or die "can't open $file: !$\n";
1182         print OUT $geno_dataset;
1183         close OUT;
1184
1185         $file_cache->set( $key_gen, $file, '30 days' );
1186         $gen_dataset_file = $file_cache->get($key_gen);
1187     }
1188
1189     return $gen_dataset_file;
1190
1191 }
1192
1193 =head2 phenotype_file
1194
1195  Usage: my $gen_file = $self->phenotype_file();
1196  Desc: creates the phenotype file in the /data/prod/tmp/r_qtl/cache,
1197        if it does not exist yet and caches it for R.
1198  Ret: phenotype filename (with abosolute path)
1199  Args: none
1200  Side Effects:
1201  Example:
1202
1203 =cut
1204
1205 sub phenotype_file
1206 {
1207     my $self       = shift;
1208     my $pop_id     = $self->get_object_id();
1209     my $population = $self->get_object();
1210
1211     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1212       $self->cache_temp_path();
1213     my $file_cache = Cache::File->new( cache_root => $prod_cache_path );
1214
1215     my $key_phe          = "popid_" . $pop_id . "_phenodata";
1216     my $phe_dataset_file = $file_cache->get($key_phe);
1217
1218     unless ($phe_dataset_file)
1219     {
1220         my $phenodata     = $population->phenotype_dataset();
1221         my $pheno_dataset = ${$phenodata};
1222         my $filename      = "phenodata_" . $pop_id . ".csv";
1223
1224         my $file = "$prod_cache_path/$filename";
1225
1226         open OUT, ">$file" or die "can't open $file: !$\n";
1227         print OUT $pheno_dataset;
1228         close OUT;
1229
1230         $file_cache->set( $key_phe, $file, '30 days' );
1231         $phe_dataset_file = $file_cache->get($key_phe);
1232     }
1233
1234     return $phe_dataset_file;
1235
1236 }
1237
1238 =head2 crosstype_file
1239
1240  Usage: my $gen_file = $self->crosstype_file();
1241  Desc: creates the crosstype file in the /data/prod/tmp/r_qtl/temp,
1242
1243  Ret: crossotype filename (with abosolute path)
1244  Args: none
1245  Side Effects:
1246  Example:
1247
1248 =cut
1249
1250 sub crosstype_file
1251 {
1252     my $self       = shift;
1253     my $pop_id     = $self->get_object_id();
1254     my $population = $self->get_object();
1255
1256     my $cross_type = 'bc' if ($population->get_cross_type_id() == 2);
1257     $cross_type = 'f2' if ($population->get_cross_type_id() == 1);
1258
1259     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1260         $self->cache_temp_path();
1261
1262     my $cross_temp = File::Temp->new(
1263                                      TEMPLATE => "cross_type_${pop_id}-XXXXXX",
1264                                      DIR      => $prod_temp_path,
1265                                      UNLINK   => 0,
1266                                     );
1267
1268
1269     my $cross_file = $cross_temp->filename;
1270
1271     open CF, ">$cross_file" or die "can't open $cross_file: $!\n";
1272     print CF $cross_type;
1273     close FO;
1274
1275     return $cross_file;
1276
1277 }
1278
1279
1280
1281 =head2 run_r
1282
1283  Usage: my ($qtl_summary, $flanking_markers) = $self->run_r();
1284  Desc: run R in the cluster; copies permutation file from the /data/prod..
1285        to the tempimages dir; returns the R output files (with abosulate filepath) with qtl mapping data
1286        and flanking markers
1287  Ret:
1288  Args:  none
1289  Side Effects:
1290  Example:
1291
1292 =cut
1293
1294 sub run_r
1295 {
1296     my $self = shift;
1297
1298     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1299       $self->cache_temp_path();
1300     my $prod_permu_file = $self->permu_file();
1301     my $file_in         = $self->infile_list();
1302     my ( $file_out, $qtl_summary, $flanking_markers ) = $self->outfile_list();
1303     my $stat_file = $self->stat_files();
1304
1305     CXGN::Tools::Run->temp_base($prod_temp_path);
1306
1307     my ( $r_in_temp, $r_out_temp ) =
1308       map {
1309         my ( undef, $filename ) =
1310           tempfile(
1311                     File::Spec->catfile(
1312                                          CXGN::Tools::Run->temp_base(),
1313                                          "population_indls.pl-$_-XXXXXX",
1314                                        ),
1315                   );
1316         $filename
1317       } qw / in out /;
1318
1319     #copy our R commands into a cluster-accessible tempfile
1320     my $doc = CXGN::Scrap::AjaxPage->new();
1321
1322     {
1323         my $r_cmd_file = $doc->path_to('/cgi-bin/phenome/cvterm_qtl.r');
1324         copy( $r_cmd_file, $r_in_temp )
1325           or die "could not copy '$r_cmd_file' to '$r_in_temp'";
1326     }
1327
1328     # now run the R job on the cluster
1329     my $r_process = CXGN::Tools::Run->run_cluster(
1330         'R', 'CMD', 'BATCH',
1331         '--slave',
1332         "--args $file_in $file_out $stat_file",
1333         $r_in_temp,
1334         $r_out_temp,
1335         {
1336            working_dir => $prod_temp_path,
1337
1338            # don't block and wait if the cluster looks full
1339            max_cluster_jobs => 1_000_000_000,
1340         },
1341     );
1342
1343     sleep 1 while $r_process->alive;    #< wait for R to finish
1344                                         #unlink( $r_in_temp, $r_out_temp );
1345
1346     copy( $prod_permu_file, $tempimages_path )
1347       or die "could not copy '$prod_permu_file' to '$tempimages_path'";
1348
1349     return $qtl_summary, $flanking_markers;
1350
1351 }
1352
1353 =head2 permu_file
1354
1355  Usage: my $permu_file = $self->permu_file();
1356  Desc: creates the permutation file in the /data/prod/tmp/r_qtl/cache,
1357        if it does not exist yet and caches it for R.
1358  Ret: permutation filename (with abosolute path)
1359  Args: none
1360  Side Effects:
1361  Example:
1362
1363 =cut
1364
1365 sub permu_file
1366 {
1367     my $self = shift;
1368     my $doc  = CXGN::Scrap::AjaxPage->new();
1369     my ( $pop_id, $cvterm_id ) =
1370       $doc->get_encoded_arguments( "population_id", "cvterm_id" );
1371
1372     my $dbh = CXGN::DB::Connection->new();
1373
1374     my $population = CXGN::Phenome::Population->new( $dbh, $pop_id );
1375     my $pop_name = $population->get_name();
1376
1377     my ( $term_obj, $term_name, $term_id );
1378
1379     if ( $population->get_web_uploaded() )
1380     {
1381         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $cvterm_id );
1382         $term_name = $term_obj->get_name();
1383         $term_id   = $term_obj->get_user_trait_id();
1384     }
1385     else
1386     {
1387         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $cvterm_id );
1388         $term_name = $term_obj->get_cvterm_name();
1389         $term_id   = $term_obj->get_cvterm_id();
1390     }
1391
1392     my $ac = $population->cvterm_acronym($term_name);
1393
1394     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1395       $self->cache_temp_path();
1396
1397     my $file_cache = Cache::File->new( cache_root => $prod_cache_path );
1398
1399     my $key_permu = "$ac" . "_" . $pop_id . "_permu";
1400     my $filename  = "permu_" . $ac . "_" . $pop_id;
1401
1402     my $permu_file = $file_cache->get($key_permu);
1403
1404     unless ($permu_file)
1405     {
1406
1407         my $permu = undef;
1408
1409         my $permu_file = "$prod_cache_path/$filename";
1410
1411         open OUT, ">$permu_file" or die "can't open $permu_file: !$\n";
1412         print OUT $permu;
1413         close OUT;
1414
1415         $file_cache->set( $key_permu, $permu_file, '30 days' );
1416         $permu_file = $file_cache->get($key_permu);
1417     }
1418
1419     return $permu_file;
1420
1421 }
1422
1423 =head2 permu_values
1424
1425  Usage: my $permu_values = $self->permu_values();
1426  Desc: reads the permutation output from R,
1427        creates a hash with the probality level as key and LOD threshold as the value,
1428
1429  Ret: a hash ref of the permutation values
1430  Args: none
1431  Side Effects:
1432  Example:
1433
1434 =cut
1435
1436 sub permu_values
1437 {
1438     my $self            = shift;
1439     my $prod_permu_file = $self->permu_file();
1440
1441     my %permu_threshold = {};
1442
1443     my $permu_file = fileparse($prod_permu_file);
1444     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1445       $self->cache_temp_path();
1446     $permu_file = File::Spec->catfile( $tempimages_path, $permu_file );
1447
1448     my $round1 = Math::Round::Var->new(0.1);
1449
1450     open PERMUTATION, "<$permu_file"
1451       or die "can't open $permu_file: !$\n";
1452
1453     my $header = <PERMUTATION>;
1454
1455     while ( my $row = <PERMUTATION> )
1456     {
1457         my ( $significance, $lod_threshold ) = split( /\t/, $row );
1458         $lod_threshold = $round1->round($lod_threshold);
1459         $permu_threshold{$significance} = $lod_threshold;
1460     }
1461
1462     close PERMUTATION;
1463
1464     return \%permu_threshold;
1465
1466 }
1467
1468 =head2 permu_values_exist
1469
1470  Usage: my $permu_value = $self->permu_values_exist();
1471  Desc: checks if there is permutation value in the permutation file.
1472  Ret: undef or some value
1473  Args: none
1474  Side Effects:
1475  Example:
1476
1477 =cut
1478
1479 sub permu_values_exist
1480 {
1481     my $self            = shift;
1482     my $prod_permu_file = $self->permu_file();
1483
1484     my ( $size, $permu_file, $permu_data, $tempimages_path, $prod_cache_path,
1485          $prod_temp_path );
1486
1487     if ($prod_permu_file)
1488     {
1489
1490         $permu_file = fileparse($prod_permu_file);
1491         ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1492           $self->cache_temp_path();
1493     }
1494
1495     if ($permu_file)
1496     {
1497
1498         $permu_file = File::Spec->catfile( $tempimages_path, $permu_file );
1499     }
1500
1501     if ( -e $permu_file )
1502     {
1503
1504         open P, "<$permu_file" or die "can't open $permu_file: !$\n";
1505         my $h = <P>;
1506         while ( $permu_data = <P> )
1507         {
1508             last if ($permu_data);
1509
1510             #       #just checking if there is data in there
1511         }
1512         close P;
1513     }
1514
1515     if ($permu_data)
1516     {
1517         return 1;
1518     }
1519     else
1520     {
1521         return 0;
1522
1523     }
1524
1525 }
1526
1527 =head2 qtl_images_exist
1528
1529  Usage: my $qtl_images_ref = $self->qtl_images_exist();
1530  Desc: checks and returns a scalar ref if the qtl plots (with thickbox and their links to the comparative viewer) exist in the cache
1531  Ret: scalar ref to the images or undef
1532  Args: none
1533  Side Effects:
1534  Example:
1535
1536 =cut
1537
1538 sub qtl_images_exist
1539 {
1540     my $self = shift;
1541     my $doc  = CXGN::Scrap::AjaxPage->new();
1542
1543     my ( $pop_id, $cvterm_id ) =
1544       $doc->get_encoded_arguments( "population_id", "cvterm_id" );
1545
1546     my $dbh = $self->get_dbh();
1547
1548     my $population = $self->get_object();
1549     my $pop_name   = $population->get_name();
1550
1551     my @linkage_groups = $population->linkage_groups();
1552     @linkage_groups = sort ( { $a <=> $b } @linkage_groups );
1553
1554     my ( $term_obj, $term_name, $term_id );
1555
1556     if ( $population->get_web_uploaded() )
1557     {
1558         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $cvterm_id );
1559         $term_name = $term_obj->get_name();
1560         $term_id   = $term_obj->get_user_trait_id();
1561     }
1562     else
1563     {
1564         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $cvterm_id );
1565         $term_name = $term_obj->get_cvterm_name();
1566         $term_id   = $term_obj->get_cvterm_id();
1567     }
1568
1569     my $ac = $population->cvterm_acronym($term_name);
1570
1571     my $vh           = SGN::Context->new();
1572     my $basepath     = $vh->get_conf("basepath");
1573     my $tempfile_dir = $vh->get_conf("tempfiles_subdir");
1574
1575     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1576       $self->cache_temp_path();
1577
1578     my $cache_tempimages = Cache::File->new( cache_root => $tempimages_path );
1579     $cache_tempimages->purge();
1580
1581     my ( $qtl_image, $image, $image_t, $image_url, $image_html, $image_t_url,
1582          $thickbox, $title );
1583
1584
1585   IMAGES: foreach my $lg (@linkage_groups)
1586     {
1587         my $cache_qtl_plot = CXGN::Tools::WebImageCache->new();
1588         $cache_qtl_plot->set_basedir($basepath);
1589         $cache_qtl_plot->set_temp_dir( $tempfile_dir . "/temp_images" );
1590
1591         my $key = "qtlplot" . $lg . "small" . $pop_id . $term_id;
1592         $cache_qtl_plot->set_key($key);
1593
1594         my $key_h_marker = "$ac" . "_pop_" . "$pop_id" . "_chr_" . $lg;
1595         my $h_marker     = $cache_tempimages->get($key_h_marker);
1596
1597         if ( $cache_qtl_plot->is_valid )
1598         {
1599             $image      = $cache_qtl_plot->get_image_tag();
1600             $image_url  = $cache_qtl_plot->get_image_url();
1601         }
1602
1603         my $cache_qtl_plot_t = CXGN::Tools::WebImageCache->new();
1604         $cache_qtl_plot_t->set_basedir($basepath);
1605         $cache_qtl_plot_t->set_temp_dir( $tempfile_dir . "/temp_images" );
1606
1607         my $key_t = "qtlplot_" . $lg . "_thickbox_" . $pop_id . $term_id;
1608         $cache_qtl_plot_t->set_key($key_t);
1609
1610         if ( $cache_qtl_plot_t->is_valid )
1611         {
1612
1613             $image_t     = $cache_qtl_plot_t->get_image_tag();
1614             $image_t_url = $cache_qtl_plot_t->get_image_url();
1615
1616             $thickbox =
1617 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> |;
1618
1619             $qtl_image .= $thickbox;
1620             $title = "  ";
1621
1622
1623         }
1624         else
1625         {
1626             $qtl_image = undef;
1627             last IMAGES;
1628
1629         }
1630
1631     }
1632
1633     return $qtl_image;
1634
1635 }
1636
1637 # =head2 user_stat_file
1638
1639 #  Usage:
1640 #  Desc:
1641 #  Ret:
1642 #  Args:
1643 #  Side Effects:
1644 #  Example:
1645
1646 # =cut
1647
1648 # sub user_stat_file {
1649 #     my $self = shift;
1650 #     my $pop = $self->get_object();
1651 #     my $pop_id = $self->get_object_id();
1652 #     my $sp_person_id = $pop->get_sp_person_id();
1653 #     my $qtl  = CXGN::Phenome::Qtl->new($sp_person_id);
1654 #     #$qtl->set_population_id($pop_id);
1655
1656 #     my ($qtl_dir, $user_dir)  = $qtl->get_user_qtl_dir();
1657
1658 #     my $stat_file = "$user_dir/user_stat_pop_$pop_id.txt";
1659 #     print STDERR "stat_file: $stat_file";
1660
1661 #     if (-e $stat_file) {
1662 #       return $stat_file;
1663 #     } else {return 0;}
1664
1665 # }
1666
1667 =head2 stat_files
1668
1669  Usage: my $stat_param_files = $self->stat_files();
1670  Desc:  creates a master file containing individual files
1671         in /data/prod/tmp/r_qtl for each statistical parameter
1672         which are feed to  R.
1673  Ret:  an absolute path to the statistical parameter's
1674        master file (and individual files)
1675  Args: None
1676  Side Effects:
1677  Example:
1678
1679 =cut
1680
1681 sub stat_files
1682 {
1683     my $self           = shift;
1684     my $pop_id         = $self->get_object_id();
1685     my $pop            = $self->get_object();
1686     my $sp_person_id   = $pop->get_sp_person_id();
1687     my $qtl            = CXGN::Phenome::Qtl->new($sp_person_id);
1688     my $c = SGN::Context->new();
1689     my $user_stat_file = $qtl->get_stat_file($c, $pop_id);
1690
1691     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1692       $self->cache_temp_path();
1693
1694     open F, "<$user_stat_file" or die "can't open file: !$\n";
1695
1696     my $stat_files;
1697
1698     while (<F>)
1699     {
1700         my ( $parameter, $value ) = split( /\t/, $_ );
1701
1702         my $stat_temp = File::Temp->new(
1703                                       TEMPLATE => "${parameter}_$pop_id-XXXXXX",
1704                                       DIR      => $prod_temp_path,
1705                                       UNLINK   => 0
1706         );
1707         my $stat_file = $stat_temp->filename;
1708
1709         open SF, ">$stat_file" or die "can't open file: !$\n";
1710         print SF $value;
1711         close SF;
1712
1713         $stat_files .= $stat_file . "\t";
1714
1715     }
1716
1717     close F;
1718
1719     my $stat_param_files =
1720       $prod_temp_path . "/" . "stat_temp_files_pop_id_${pop_id}";
1721
1722     open STAT, ">$stat_param_files" or die "can't open file: !$\n";
1723     print STAT $stat_files;
1724     close STAT;
1725
1726     return $stat_param_files;
1727
1728 }
1729
1730 =head2 stat_param_hash
1731
1732  Usage: my %stat_param = $self->stat_param_hash();
1733  Desc: creates a hash (with the statistical parameters (as key) and
1734        their corresponding values) out of a tab delimited
1735        statistical parameters file.
1736  Ret: a hash statistics file
1737  Args: None
1738  Side Effects:
1739  Example:
1740
1741 =cut
1742
1743 sub stat_param_hash
1744 {
1745     my $self           = shift;
1746     my $pop_id         = $self->get_object_id();
1747     my $pop            = $self->get_object();
1748     my $sp_person_id   = $pop->get_sp_person_id();
1749     my $qtl            = CXGN::Phenome::Qtl->new($sp_person_id);
1750     my $c = SGN::Context->new();
1751     my $user_stat_file = $qtl->get_stat_file($c, $pop_id);
1752
1753     open F, "<$user_stat_file" or die "can't open file: !$\n";
1754
1755     my %stat_param;
1756
1757     while (<F>)
1758     {
1759         my ( $parameter, $value ) = split( /\t/, $_ );
1760
1761         $stat_param{$parameter} = $value;
1762
1763     }
1764
1765     return \%stat_param;
1766 }
1767
1768 sub submitter
1769 {
1770     my $self         = shift;
1771     my $population   = $self->get_object();
1772     my $sp_person_id = $population->get_sp_person_id();
1773     my $submitter    = CXGN::People::Person->new( $self->get_dbh(),
1774                                               $population->get_sp_person_id() );
1775     my $submitter_name =
1776       $submitter->get_first_name() . " " . $submitter->get_last_name();
1777     my $submitter_link =
1778 qq |<a href="/solpeople/personal-info.pl?sp_person_id=$sp_person_id">$submitter_name</a> |;
1779
1780     return $submitter, $submitter_link;
1781
1782 }
1783
1784 #move to qtl or population object
1785 sub legend {
1786     my $self = shift;
1787     my $pop = shift;
1788     my $sp_person_id   = $pop->get_sp_person_id();
1789     my $qtl            = CXGN::Phenome::Qtl->new($sp_person_id);
1790     my $stat_file = $qtl->get_stat_file($c, $pop->get_population_id());
1791     my @stat;
1792     my $ci;
1793
1794     open my $sf, "<", $stat_file or die "$! reading $stat_file\n";
1795     while (my $row = <$sf>)
1796     {
1797         chomp($row);
1798         my ( $parameter, $value ) = split( /\t/, $row );
1799         if ($parameter =~/qtl_method/) {$parameter = 'Mapping method';}
1800         if ($parameter =~/qtl_model/) {$parameter = 'Mapping model';}
1801         if ($parameter =~/prob_method/) {$parameter = 'QTL genotype probability method';}
1802         if ($parameter =~/step_size/) {$parameter = 'Genome scan size (cM)';}
1803         if ($parameter =~/permu_level/) {$parameter = 'Permutation significance level';}
1804         if ($parameter =~/permu_test/) {$parameter = 'No. of permutations';}
1805         if ($parameter =~/prob_level/) {$parameter = 'QTL genotype signifance level';}
1806         if ($parameter =~/stat_no_draws/) {$parameter = 'No. of imputations';}
1807
1808         if ($value eq 'zero' || $value eq 'Marker Regression') {$ci = 'none';}
1809
1810         unless (($parameter =~/No. of imputations/ && !$value ) ||
1811                 ($parameter =~/QTL genotype probability/ && !$value ) ||
1812                 ($parameter =~/Permutation significance level/ && !$value)
1813                )
1814
1815         {
1816             push @stat, [map{$_} ($parameter, $value)];
1817
1818         }
1819
1820     }
1821
1822     my $sm;
1823     foreach my $st (@stat) {
1824         foreach my $i (@$st) {
1825             if ($i =~/zero/)  {
1826                 foreach my $s (@stat) {
1827                     foreach my $j (@$s) {
1828                         $j =~ s/Maximum Likelihood/Marker Regression/;
1829                         $ci = 'none';
1830                     }
1831                 }
1832             }
1833         }
1834     }
1835
1836
1837 my $permu_threshold_ref = $self->permu_values();
1838     my %permu_threshold     = %$permu_threshold_ref;
1839
1840     my @keys;
1841
1842     foreach my $key ( keys %permu_threshold )
1843     {
1844         if ( $key =~ m/^\d./ )
1845         {
1846             push @keys, $key;
1847         }
1848
1849     }
1850     my $lod1 = $permu_threshold{ $keys[0] };
1851    # my $lod2 = $permu_threshold{ $keys[1] };
1852
1853     if  (!$lod1)
1854     {
1855         $lod1 = qq |<i>Not calculated</i>|;
1856     }
1857
1858     push @stat,
1859     [
1860      map {$_} ('LOD threshold', $lod1)
1861     ];
1862
1863
1864     unless ($ci) {
1865         push @stat,
1866         [
1867          map {$_} ('Confidence interval', 'Based on 95% Bayesian Credible Interval')
1868         ];
1869     }
1870     push @stat,
1871     [
1872      map {$_} ('QTL software', "<a href=http://www.rqtl.org>R/QTL</a>")
1873     ];
1874     my $legend_data = columnar_table_html (
1875                                            headings    => [
1876                                                           ' ',
1877                                                           ' ',
1878
1879                                                          ],
1880                                            data        => \@stat,
1881                                           __alt_freq   => 2,
1882                                           __alt_width  => 1,
1883                                           __alt_offset => 3,
1884                                           __align      => 'l',
1885
1886                                           );
1887
1888
1889     return $legend_data;
1890
1891 }