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,   $normal_dist
 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 . $normal_dist,
 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                 # my $log2 = $permu_threshold{ $p_keys[1] };
 787
 788                 $h_marker =
 789                     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|;
 790  #$h_marker =
 791 #qq |../cview/view_chromosome.pl?map_version_id=$mapversion&chr_nr=$lg&show_ruler=1&show_IL=&show_offsets=1&comp_map_version_id=&comp_chr=&color_model=&show_physical=&size=&show_zoomed=1&confidence=-2&hilite=$l_m+$p_m+$r_m&marker_type=&cM_start=$l_pos&cM_end=$r_pos |;
 792
 793                 $cache_tempimages->set( $key_h_marker, $h_marker, '30 days' );
 794             }
 795
 796             push @h_markers, $h_marker;
 797         }
 798
 799         my $count       = 0;
 800         my $old_chr_chr = 1;
 801         my (
 802              $chr_chr, $image,     $image_t,
 803              $thickbox,  $max_chr, $chr_chr_e, $marker_chr_e,
 804              $pos_chr_e, $lod_chr_e
 805            );
 806         my $chrs = ( scalar(@chromosomes) ) + 1;
 807
 808         for ( my $i = 1 ; $i < $chrs ; $i++ )
 809         {
 810             my ( @marker_chr, @chr_chr, @pos_chr, @lod_chr, @data, @m_html ) =
 811               ();
 812             my ( $marker_chr, $pos_chr, $lod_chr, $max_chr );
 813
 814             $h_marker = shift(@h_markers);
 815
 816             if ( ( $i == $old_chr_chr ) && ( $i != 12 ) )
 817             {
 818                 push @marker_chr, $marker_chr_e;
 819                 push @chr_chr,    $chr_chr_e;
 820                 $pos_chr_e = $round2->round($pos_chr_e);
 821                 push @pos_chr, $pos_chr_e;
 822                 $lod_chr = $round1->round($lod_chr_e);
 823                 push @lod_chr, $lod_chr_e;
 824             }
 825
 826             my $cache_qtl_plot = CXGN::Tools::WebImageCache->new();
 827             $cache_qtl_plot->set_basedir($basepath);
 828             $cache_qtl_plot->set_temp_dir( $tempfile_dir . "/temp_images" );
 829             $cache_qtl_plot->set_expiration_time(259200);
 830             $cache_qtl_plot->set_key(
 831                                 "qtlplot" . $i . "small" . $pop_id . $term_id );
 832             $cache_qtl_plot->set_force(0);
 833
 834             if ( !$cache_qtl_plot->is_valid() )
 835             {
 836
 837                 for ( my $j = 0 ; $j < @marker ; $j++ )
 838                 {
 839
 840                     $chr_chr = $chr[$j];
 841
 842                     if ( $i == $chr_chr )
 843                     {
 844                         $marker_chr = $marker[$j];
 845
 846                         $pos_chr = $pos[$j];
 847                         $lod_chr = $lod[$j];
 848
 849                         push @marker_chr, $marker_chr;
 850                         push @chr_chr,    $chr_chr;
 851                         $pos_chr = $round2->round($pos_chr);
 852                         push @pos_chr, $pos_chr;
 853                         $lod_chr = $round1->round($lod_chr);
 854                         push @lod_chr, $lod_chr;
 855
 856                         ( $chr_chr_e, $marker_chr_e, $pos_chr_e, $lod_chr_e ) =
 857                           ();
 858                     }
 859
 860                     elsif ( $i != $chr_chr )
 861                     {
 862
 863                         $chr_chr_e    = $chr[$j];
 864                         $marker_chr_e = $marker[$j];
 865                         $pos_chr_e    = $pos[$j];
 866                         $lod_chr_e    = $lod[$j];
 867                     }
 868
 869                 }
 870
 871                 @data = ( [ (@pos_chr) ], [@lod_chr] );
 872                 my $graph = new GD::Graph::lines( 110, 110 );
 873                 $graph->set_title_font('gdTinyFont');
 874                 $graph->set(
 875                              title             => " ",
 876                              x_label           => "Chr $i (cM)",
 877                              y_label           => "LOD",
 878                              y_max_value       => 10,
 879                              x_all_ticks       => 5,
 880                              y_all_ticks       => 1,
 881                              y_label_skip      => 1,
 882                              y_plot_values     => 1,
 883                              x_label_skip      => 5,
 884                              x_plot_values     => 1,
 885                              x_labels_vertical => 1,
 886                              textclr           => "black"
 887                            );
 888
 889                 $cache_qtl_plot->set_image_data( $graph->plot( \@data )->png );
 890
 891             }
 892
 893             $image      = $cache_qtl_plot->get_image_tag();
 894             $image_url  = $cache_qtl_plot->get_image_url();
 895            # $image_html = qq |<a href ="$h_marker&qtl=$image_url">$image</a>|;
 896
 897 ###########thickbox
 898             my $cache_qtl_plot_t = CXGN::Tools::WebImageCache->new();
 899             $cache_qtl_plot_t->set_basedir($basepath);
 900             $cache_qtl_plot_t->set_temp_dir( $tempfile_dir . "/temp_images" );
 901             $cache_qtl_plot_t->set_expiration_time(259200);
 902             $cache_qtl_plot_t->set_key(
 903                           "qtlplot_" . $i . "_thickbox_" . $pop_id . $term_id );
 904             $cache_qtl_plot_t->set_force(0);
 905
 906             if ( !$cache_qtl_plot_t->is_valid() )
 907             {
 908                 my @o_lod_chr = sort { $a <=> $b } @lod_chr;
 909                 $max_chr = pop(@o_lod_chr);
 910                 $max_chr = $max_chr + (0.5);
 911
 912                 my $graph_t = new GD::Graph::lines( 420, 420 );
 913                 $graph_t->set_title_font('gdTinyFont');
 914                 $graph_t->set(
 915                                title             => " ",
 916                                x_label           => "Chromosome $i (cM)",
 917                                y_label           => "LOD",
 918                                y_max_value       => $max_chr,
 919                                x_all_ticks       => 5,
 920                                y_all_ticks       => 1,
 921                                y_label_skip      => 1,
 922                                y_plot_values     => 1,
 923                                x_label_skip      => 5,
 924                                x_plot_values     => 1,
 925                                x_labels_vertical => 1,
 926                                textclr           => "black"
 927                              );
 928
 929                 $cache_qtl_plot_t->set_image_data(
 930                                                 $graph_t->plot( \@data )->png );
 931
 932             }
 933
 934             $image_t     = $cache_qtl_plot_t->get_image_tag();
 935             $image_t_url = $cache_qtl_plot_t->get_image_url();
 936
 937             $thickbox =
 938 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> |;
 939
 940             $qtl_image .= $thickbox;
 941             $title       = "  ";
 942             $old_chr_chr = $chr_chr;
 943         }
 944     }
 945
 946     return $qtl_image;
 947 }
 948
 949 =head2 infile_list
 950  Usage: my $file_in = $self->infile_list();
 951  Desc: returns an R input tempfile containing a tempfile
 952        holding the cvterm acronym, pop id, a filepath to the phenotype dataset file,
 953         a filepath to genotype dataset file, a filepath to the permuation file.
 954  Ret: an R input tempfile name (with abosulte path)
 955  Args:
 956  Side Effects:
 957  Example:
 958
 959
 960 =cut
 961
 962 sub infile_list
 963 {
 964
 965     my $self = shift;
 966     my $doc  = CXGN::Scrap::AjaxPage->new();
 967
 968     my ( $pop_id, $cvterm_id ) =
 969       $doc->get_encoded_arguments( "population_id", "cvterm_id" );
 970
 971     my $dbh = $self->get_dbh();
 972
 973     my ( $term_obj, $term_name, $term_id );
 974     my $population = $self->get_object();
 975
 976     if ( $population->get_web_uploaded() )
 977     {
 978         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $cvterm_id );
 979         $term_name = $term_obj->get_name();
 980         $term_id   = $term_obj->get_user_trait_id();
 981     }
 982     else
 983     {
 984         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $cvterm_id );
 985         $term_name = $term_obj->get_cvterm_name();
 986         $term_id   = $term_obj->get_cvterm_id();
 987     }
 988
 989     my $ac = $population->cvterm_acronym($term_name);
 990
 991     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
 992       $self->cache_temp_path();
 993
 994     my $prod_permu_file  = $self->permu_file();
 995     my $gen_dataset_file = $self->genotype_file();
 996     my $phe_dataset_file = $self->phenotype_file();
 997     my $crosstype_file   = $self->crosstype_file();
 998
 999     my $input_file_list_temp =
1000       File::Temp->new(
1001                        TEMPLATE => "infile_list_${ac}_$pop_id-XXXXXX",
1002                        DIR      => $prod_temp_path,
1003                        UNLINK   => 0,
1004                      );
1005     my $file_in = $input_file_list_temp->filename();
1006
1007     my $file_cvin = File::Temp->new(
1008                                      TEMPLATE => 'cv_input-XXXXXX',
1009                                      DIR      => $prod_temp_path,
1010                                      UNLINK   => 0,
1011                                    );
1012     my $file_cv_in = $file_cvin->filename();
1013
1014     open CV, ">$file_cv_in" or die "can't open $file_cv_in: $!\n";
1015     print CV $ac;
1016     close CV;
1017
1018     my $file_in_list = join( "\t",
1019                              $file_cv_in,       "P$pop_id",
1020                              $gen_dataset_file, $phe_dataset_file,
1021                              $prod_permu_file, $crosstype_file);
1022
1023     open FI, ">$file_in" or die "can't open $file_in: $!\n";
1024     print FI $file_in_list;
1025     close FI;
1026
1027     return $file_in;
1028
1029 }
1030
1031 =head2 outfile_list
1032
1033  Usage: my ($file_out, $qtl_summary, $flanking_markers) = $self->outfile_list();
1034  Desc: returns an R output tempfile containing a tempfile supposed to hold the qtl
1035        mapping output and another tempfile for the qtl flanking markers
1036        and the qtl mapping output and qtl flanking markers files separately
1037        (convenient for reading their data when plotting the qtl)
1038  Ret: R output file names (with abosulte path)
1039  Args:
1040  Side Effects:
1041  Example:
1042
1043 =cut
1044
1045 sub outfile_list
1046 {
1047     my $self = shift;
1048
1049     my $doc = CXGN::Scrap::AjaxPage->new();
1050
1051     my ( $pop_id, $cvterm_id ) =
1052       $doc->get_encoded_arguments( "population_id", "cvterm_id" );
1053
1054     my $dbh = $self->get_dbh();
1055
1056     my ( $term_obj, $term_name, $term_id );
1057     my $population = $self->get_object();
1058
1059     if ( $population->get_web_uploaded() )
1060     {
1061         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $cvterm_id );
1062         $term_name = $term_obj->get_name();
1063         $term_id   = $term_obj->get_user_trait_id();
1064     }
1065     else
1066     {
1067         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $cvterm_id );
1068         $term_name = $term_obj->get_cvterm_name();
1069         $term_id   = $term_obj->get_cvterm_id();
1070     }
1071
1072     my $ac = $population->cvterm_acronym($term_name);
1073
1074     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1075       $self->cache_temp_path();
1076
1077     my $output_file_list_temp =
1078       File::Temp->new(
1079                        TEMPLATE => "outfile_list_${ac}_$pop_id-XXXXXX",
1080                        DIR      => $prod_temp_path,
1081                        UNLINK   => 0,
1082                      );
1083     my $file_out = $output_file_list_temp->filename();
1084
1085     my $qtl_temp = File::Temp->new(
1086                                  TEMPLATE => "qtl_summary_${ac}_$pop_id-XXXXXX",
1087                                  DIR      => $prod_temp_path,
1088                                  UNLINK   => 0
1089     );
1090     my $qtl_summary = $qtl_temp->filename;
1091
1092     my $marker_temp = File::Temp->new(
1093                             TEMPLATE => "flanking_markers_${ac}_$pop_id-XXXXXX",
1094                             DIR      => $prod_temp_path,
1095                             UNLINK   => 0
1096     );
1097
1098     my $flanking_markers = $marker_temp->filename;
1099
1100     my $file_out_list = join(
1101         "\t",
1102         $qtl_summary,
1103         $flanking_markers,
1104
1105         #$qtl_summary_file,
1106         #$flanking_markers_file
1107                             );
1108     open FO, ">$file_out" or die "can't open $file_out: $!\n";
1109     print FO $file_out_list;
1110     close FO;
1111
1112     return $file_out, $qtl_summary, $flanking_markers;
1113 }
1114
1115 =head2 cache_temp_path
1116
1117  Usage: my ($prod_cache_path, $prod_temp_path, $tempimages_path) = $self->cache_temp_path();
1118  Desc: creates the 'r_qtl' dir in the '/data/prod/tmp/' dir;
1119       'cache' and 'tempfiles' in the /data/prod/tmp/r_qtl/,
1120        and 'temp_images' in the /data/local/cxgn/sgn/documents/tempfiles'
1121  Ret: /data/prod/tmp/r_qtl/cache, /data/prod/tmp/r_qtl/tempfiles,
1122       /data/local/cxgn/sgn/documents/tempfiles/temp_images
1123  Args: none
1124  Side Effects:
1125  Example:
1126
1127 =cut
1128
1129 sub cache_temp_path
1130 {
1131     my $vh           = SGN::Context->new();
1132     my $basepath     = $vh->get_conf("basepath");
1133     my $tempfile_dir = $vh->get_conf("tempfiles_subdir");
1134
1135     my $tempimages_path =
1136       File::Spec->catfile( $basepath, $tempfile_dir, "temp_images" );
1137
1138     my $prod_temp_path = $vh->get_conf('r_qtl_temp_path');
1139     mkdir $prod_temp_path;
1140     my $prod_cache_path = "$prod_temp_path/cache";
1141     mkdir $prod_cache_path;
1142     $prod_temp_path = "$prod_temp_path/tempfiles";
1143     mkdir $prod_temp_path;
1144     -d $prod_temp_path
1145       or die "temp dir '$prod_temp_path' not found, and could not create!";
1146     -r $prod_temp_path or die "temp dir '$prod_temp_path' not readable!";
1147     -w $prod_temp_path or die "temp dir '$prod_temp_path' not writable!";
1148
1149     return $prod_cache_path, $prod_temp_path, $tempimages_path;
1150
1151 }
1152
1153 =head2 genotype_file
1154
1155  Usage: my $gen_file = $self->genotype_file();
1156  Desc: creates the genotype file in the /data/prod/tmp/r_qtl/cache,
1157        if it does not exist yet and caches it for R.
1158  Ret: genotype filename (with abosolute path)
1159  Args: none
1160  Side Effects:
1161  Example:
1162
1163 =cut
1164
1165 sub genotype_file
1166 {
1167     my $self       = shift;
1168     my $pop_id     = $self->get_object_id();
1169     my $population = $self->get_object();
1170
1171     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1172       $self->cache_temp_path();
1173     my $file_cache = Cache::File->new( cache_root => $prod_cache_path );
1174     $file_cache->purge();
1175
1176     my $key_gen          = "popid_" . $pop_id . "_genodata";
1177     my $gen_dataset_file = $file_cache->get($key_gen);
1178
1179     unless ($gen_dataset_file)
1180     {
1181         my $genodata     = $population->genotype_dataset();
1182         my $geno_dataset = ${$genodata};
1183
1184         my $filename = "genodata_" . $pop_id . ".csv";
1185         my $file     = "$prod_cache_path/$filename";
1186
1187         open OUT, ">$file" or die "can't open $file: !$\n";
1188         print OUT $geno_dataset;
1189         close OUT;
1190
1191         $file_cache->set( $key_gen, $file, '30 days' );
1192         $gen_dataset_file = $file_cache->get($key_gen);
1193     }
1194
1195     return $gen_dataset_file;
1196
1197 }
1198
1199 =head2 phenotype_file
1200
1201  Usage: my $gen_file = $self->phenotype_file();
1202  Desc: creates the phenotype file in the /data/prod/tmp/r_qtl/cache,
1203        if it does not exist yet and caches it for R.
1204  Ret: phenotype filename (with abosolute path)
1205  Args: none
1206  Side Effects:
1207  Example:
1208
1209 =cut
1210
1211 sub phenotype_file
1212 {
1213     my $self       = shift;
1214     my $pop_id     = $self->get_object_id();
1215     my $population = $self->get_object();
1216
1217     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1218       $self->cache_temp_path();
1219     my $file_cache = Cache::File->new( cache_root => $prod_cache_path );
1220
1221     my $key_phe          = "popid_" . $pop_id . "_phenodata";
1222     my $phe_dataset_file = $file_cache->get($key_phe);
1223
1224     unless ($phe_dataset_file)
1225     {
1226         my $phenodata     = $population->phenotype_dataset();
1227         my $pheno_dataset = ${$phenodata};
1228         my $filename      = "phenodata_" . $pop_id . ".csv";
1229
1230         my $file = "$prod_cache_path/$filename";
1231
1232         open OUT, ">$file" or die "can't open $file: !$\n";
1233         print OUT $pheno_dataset;
1234         close OUT;
1235
1236         $file_cache->set( $key_phe, $file, '30 days' );
1237         $phe_dataset_file = $file_cache->get($key_phe);
1238     }
1239
1240     return $phe_dataset_file;
1241
1242 }
1243
1244 =head2 crosstype_file
1245
1246  Usage: my $gen_file = $self->crosstype_file();
1247  Desc: creates the crosstype file in the /data/prod/tmp/r_qtl/temp,
1248
1249  Ret: crossotype filename (with abosolute path)
1250  Args: none
1251  Side Effects:
1252  Example:
1253
1254 =cut
1255
1256 sub crosstype_file
1257 {
1258     my $self       = shift;
1259     my $pop_id     = $self->get_object_id();
1260     my $population = $self->get_object();
1261
1262     my $cross_type = 'bc' if ($population->get_cross_type_id() == 2);
1263     $cross_type = 'f2' if ($population->get_cross_type_id() == 1);
1264
1265     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1266         $self->cache_temp_path();
1267
1268     my $cross_temp = File::Temp->new(
1269                                      TEMPLATE => "cross_type_${pop_id}-XXXXXX",
1270                                      DIR      => $prod_temp_path,
1271                                      UNLINK   => 0,
1272                                     );
1273
1274
1275     my $cross_file = $cross_temp->filename;
1276
1277     open CF, ">$cross_file" or die "can't open $cross_file: $!\n";
1278     print CF $cross_type;
1279     close FO;
1280
1281     return $cross_file;
1282
1283 }
1284
1285
1286
1287 =head2 run_r
1288
1289  Usage: my ($qtl_summary, $flanking_markers) = $self->run_r();
1290  Desc: run R in the cluster; copies permutation file from the /data/prod..
1291        to the tempimages dir; returns the R output files (with abosulate filepath) with qtl mapping data
1292        and flanking markers
1293  Ret:
1294  Args:  none
1295  Side Effects:
1296  Example:
1297
1298 =cut
1299
1300 sub run_r
1301 {
1302     my $self = shift;
1303
1304     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1305       $self->cache_temp_path();
1306     my $prod_permu_file = $self->permu_file();
1307     my $file_in         = $self->infile_list();
1308     my ( $file_out, $qtl_summary, $flanking_markers ) = $self->outfile_list();
1309     my $stat_file = $self->stat_files();
1310
1311     CXGN::Tools::Run->temp_base($prod_temp_path);
1312
1313     my ( $r_in_temp, $r_out_temp ) =
1314       map {
1315         my ( undef, $filename ) =
1316           tempfile(
1317                     File::Spec->catfile(
1318                                          CXGN::Tools::Run->temp_base(),
1319                                          "population_indls.pl-$_-XXXXXX",
1320                                        ),
1321                   );
1322         $filename
1323       } qw / in out /;
1324
1325     #copy our R commands into a cluster-accessible tempfile
1326     my $doc = CXGN::Scrap::AjaxPage->new();
1327
1328     {
1329         my $r_cmd_file = $doc->path_to('/cgi-bin/phenome/cvterm_qtl.r');
1330         copy( $r_cmd_file, $r_in_temp )
1331           or die "could not copy '$r_cmd_file' to '$r_in_temp'";
1332     }
1333
1334     # now run the R job on the cluster
1335     my $r_process = CXGN::Tools::Run->run_cluster(
1336         'R', 'CMD', 'BATCH',
1337         '--slave',
1338         "--args $file_in $file_out $stat_file",
1339         $r_in_temp,
1340         $r_out_temp,
1341         {
1342            working_dir => $prod_temp_path,
1343
1344            # don't block and wait if the cluster looks full
1345            max_cluster_jobs => 1_000_000_000,
1346         },
1347     );
1348
1349     sleep 1 while $r_process->alive;    #< wait for R to finish
1350                                         #unlink( $r_in_temp, $r_out_temp );
1351
1352     copy( $prod_permu_file, $tempimages_path )
1353       or die "could not copy '$prod_permu_file' to '$tempimages_path'";
1354
1355     return $qtl_summary, $flanking_markers;
1356
1357 }
1358
1359 =head2 permu_file
1360
1361  Usage: my $permu_file = $self->permu_file();
1362  Desc: creates the permutation file in the /data/prod/tmp/r_qtl/cache,
1363        if it does not exist yet and caches it for R.
1364  Ret: permutation filename (with abosolute path)
1365  Args: none
1366  Side Effects:
1367  Example:
1368
1369 =cut
1370
1371 sub permu_file
1372 {
1373     my $self = shift;
1374     my $doc  = CXGN::Scrap::AjaxPage->new();
1375     my ( $pop_id, $cvterm_id ) =
1376       $doc->get_encoded_arguments( "population_id", "cvterm_id" );
1377
1378     my $dbh = CXGN::DB::Connection->new();
1379
1380     my $population = CXGN::Phenome::Population->new( $dbh, $pop_id );
1381     my $pop_name = $population->get_name();
1382
1383     my ( $term_obj, $term_name, $term_id );
1384
1385     if ( $population->get_web_uploaded() )
1386     {
1387         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $cvterm_id );
1388         $term_name = $term_obj->get_name();
1389         $term_id   = $term_obj->get_user_trait_id();
1390     }
1391     else
1392     {
1393         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $cvterm_id );
1394         $term_name = $term_obj->get_cvterm_name();
1395         $term_id   = $term_obj->get_cvterm_id();
1396     }
1397
1398     my $ac = $population->cvterm_acronym($term_name);
1399
1400     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1401       $self->cache_temp_path();
1402
1403     my $file_cache = Cache::File->new( cache_root => $prod_cache_path );
1404
1405     my $key_permu = "$ac" . "_" . $pop_id . "_permu";
1406     my $filename  = "permu_" . $ac . "_" . $pop_id;
1407
1408     my $permu_file = $file_cache->get($key_permu);
1409
1410     unless ($permu_file)
1411     {
1412
1413         my $permu = undef;
1414
1415         my $permu_file = "$prod_cache_path/$filename";
1416
1417         open OUT, ">$permu_file" or die "can't open $permu_file: !$\n";
1418         print OUT $permu;
1419         close OUT;
1420
1421         $file_cache->set( $key_permu, $permu_file, '30 days' );
1422         $permu_file = $file_cache->get($key_permu);
1423     }
1424
1425     return $permu_file;
1426
1427 }
1428
1429 =head2 permu_values
1430
1431  Usage: my $permu_values = $self->permu_values();
1432  Desc: reads the permutation output from R,
1433        creates a hash with the probality level as key and LOD threshold as the value,
1434
1435  Ret: a hash ref of the permutation values
1436  Args: none
1437  Side Effects:
1438  Example:
1439
1440 =cut
1441
1442 sub permu_values
1443 {
1444     my $self            = shift;
1445     my $prod_permu_file = $self->permu_file();
1446
1447     my %permu_threshold = {};
1448
1449     my $permu_file = fileparse($prod_permu_file);
1450     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1451       $self->cache_temp_path();
1452     $permu_file = File::Spec->catfile( $tempimages_path, $permu_file );
1453
1454     my $round1 = Math::Round::Var->new(0.1);
1455
1456     open PERMUTATION, "<$permu_file"
1457       or die "can't open $permu_file: !$\n";
1458
1459     my $header = <PERMUTATION>;
1460
1461     while ( my $row = <PERMUTATION> )
1462     {
1463         my ( $significance, $lod_threshold ) = split( /\t/, $row );
1464         $lod_threshold = $round1->round($lod_threshold);
1465         $permu_threshold{$significance} = $lod_threshold;
1466     }
1467
1468     close PERMUTATION;
1469
1470     return \%permu_threshold;
1471
1472 }
1473
1474 =head2 permu_values_exist
1475
1476  Usage: my $permu_value = $self->permu_values_exist();
1477  Desc: checks if there is permutation value in the permutation file.
1478  Ret: undef or some value
1479  Args: none
1480  Side Effects:
1481  Example:
1482
1483 =cut
1484
1485 sub permu_values_exist
1486 {
1487     my $self            = shift;
1488     my $prod_permu_file = $self->permu_file();
1489
1490     my ( $size, $permu_file, $permu_data, $tempimages_path, $prod_cache_path,
1491          $prod_temp_path );
1492
1493     if ($prod_permu_file)
1494     {
1495
1496         $permu_file = fileparse($prod_permu_file);
1497         ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1498           $self->cache_temp_path();
1499     }
1500
1501     if ($permu_file)
1502     {
1503
1504         $permu_file = File::Spec->catfile( $tempimages_path, $permu_file );
1505     }
1506
1507     if ( -e $permu_file )
1508     {
1509
1510         open P, "<$permu_file" or die "can't open $permu_file: !$\n";
1511         my $h = <P>;
1512         while ( $permu_data = <P> )
1513         {
1514             last if ($permu_data);
1515
1516             #       #just checking if there is data in there
1517         }
1518         close P;
1519     }
1520
1521     if ($permu_data)
1522     {
1523         return 1;
1524     }
1525     else
1526     {
1527         return 0;
1528
1529     }
1530
1531 }
1532
1533 =head2 qtl_images_exist
1534
1535  Usage: my $qtl_images_ref = $self->qtl_images_exist();
1536  Desc: checks and returns a scalar ref if the qtl plots (with thickbox and their links to the comparative viewer) exist in the cache
1537  Ret: scalar ref to the images or undef
1538  Args: none
1539  Side Effects:
1540  Example:
1541
1542 =cut
1543
1544 sub qtl_images_exist
1545 {
1546     my $self = shift;
1547     my $doc  = CXGN::Scrap::AjaxPage->new();
1548
1549     my ( $pop_id, $cvterm_id ) =
1550       $doc->get_encoded_arguments( "population_id", "cvterm_id" );
1551
1552     my $dbh = $self->get_dbh();
1553
1554     my $population = $self->get_object();
1555     my $pop_name   = $population->get_name();
1556
1557     my @linkage_groups = $population->linkage_groups();
1558     @linkage_groups = sort ( { $a <=> $b } @linkage_groups );
1559
1560     my ( $term_obj, $term_name, $term_id );
1561
1562     if ( $population->get_web_uploaded() )
1563     {
1564         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $cvterm_id );
1565         $term_name = $term_obj->get_name();
1566         $term_id   = $term_obj->get_user_trait_id();
1567     }
1568     else
1569     {
1570         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $cvterm_id );
1571         $term_name = $term_obj->get_cvterm_name();
1572         $term_id   = $term_obj->get_cvterm_id();
1573     }
1574
1575     my $ac = $population->cvterm_acronym($term_name);
1576
1577     my $vh           = SGN::Context->new();
1578     my $basepath     = $vh->get_conf("basepath");
1579     my $tempfile_dir = $vh->get_conf("tempfiles_subdir");
1580
1581     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1582       $self->cache_temp_path();
1583
1584     my $cache_tempimages = Cache::File->new( cache_root => $tempimages_path );
1585     $cache_tempimages->purge();
1586
1587     my ( $qtl_image, $image, $image_t, $image_url, $image_html, $image_t_url,
1588          $thickbox, $title );
1589
1590     # my $chrs = scalar(@linkage_groups) + 1;
1591
1592   IMAGES: foreach my $lg (@linkage_groups)
1593     {
1594         my $cache_qtl_plot = CXGN::Tools::WebImageCache->new();
1595         $cache_qtl_plot->set_basedir($basepath);
1596         $cache_qtl_plot->set_temp_dir( $tempfile_dir . "/temp_images" );
1597
1598         my $key = "qtlplot" . $lg . "small" . $pop_id . $term_id;
1599         $cache_qtl_plot->set_key($key);
1600
1601         my $key_h_marker = "$ac" . "_pop_" . "$pop_id" . "_chr_" . $lg;
1602         my $h_marker     = $cache_tempimages->get($key_h_marker);
1603
1604         if ( $cache_qtl_plot->is_valid )
1605         {
1606             $image      = $cache_qtl_plot->get_image_tag();
1607             $image_url  = $cache_qtl_plot->get_image_url();
1608            # $image_html = qq |<a href ="$h_marker&$image_url">$image</a>|;
1609
1610         }
1611
1612         my $cache_qtl_plot_t = CXGN::Tools::WebImageCache->new();
1613         $cache_qtl_plot_t->set_basedir($basepath);
1614         $cache_qtl_plot_t->set_temp_dir( $tempfile_dir . "/temp_images" );
1615
1616         my $key_t = "qtlplot_" . $lg . "_thickbox_" . $pop_id . $term_id;
1617         $cache_qtl_plot_t->set_key($key_t);
1618
1619         if ( $cache_qtl_plot_t->is_valid )
1620         {
1621
1622             $image_t     = $cache_qtl_plot_t->get_image_tag();
1623             $image_t_url = $cache_qtl_plot_t->get_image_url();
1624
1625             $thickbox =
1626 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> |;
1627
1628             $qtl_image .= $thickbox;
1629             $title = "  ";
1630
1631
1632         }
1633         else
1634         {
1635             $qtl_image = undef;
1636             last IMAGES;
1637
1638         }
1639
1640     }
1641
1642     return $qtl_image;
1643
1644 }
1645
1646 # =head2 user_stat_file
1647
1648 #  Usage:
1649 #  Desc:
1650 #  Ret:
1651 #  Args:
1652 #  Side Effects:
1653 #  Example:
1654
1655 # =cut
1656
1657 # sub user_stat_file {
1658 #     my $self = shift;
1659 #     my $pop = $self->get_object();
1660 #     my $pop_id = $self->get_object_id();
1661 #     my $sp_person_id = $pop->get_sp_person_id();
1662 #     my $qtl  = CXGN::Phenome::Qtl->new($sp_person_id);
1663 #     #$qtl->set_population_id($pop_id);
1664
1665 #     my ($qtl_dir, $user_dir)  = $qtl->get_user_qtl_dir();
1666
1667 #     my $stat_file = "$user_dir/user_stat_pop_$pop_id.txt";
1668 #     print STDERR "stat_file: $stat_file";
1669
1670 #     if (-e $stat_file) {
1671 #       return $stat_file;
1672 #     } else {return 0;}
1673
1674 # }
1675
1676 =head2 stat_files
1677
1678  Usage: my $stat_param_files = $self->stat_files();
1679  Desc:  creates a master file containing individual files
1680         in /data/prod/tmp/r_qtl for each statistical parameter
1681         which are feed to  R.
1682  Ret:  an absolute path to the statistical parameter's
1683        master file (and individual files)
1684  Args: None
1685  Side Effects:
1686  Example:
1687
1688 =cut
1689
1690 sub stat_files
1691 {
1692     my $self           = shift;
1693     my $pop_id         = $self->get_object_id();
1694     my $pop            = $self->get_object();
1695     my $sp_person_id   = $pop->get_sp_person_id();
1696     my $qtl            = CXGN::Phenome::Qtl->new($sp_person_id);
1697     my $c = SGN::Context->new();
1698     my $user_stat_file = $qtl->get_stat_file($c, $pop_id);
1699
1700     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1701       $self->cache_temp_path();
1702
1703     open F, "<$user_stat_file" or die "can't open file: !$\n";
1704
1705     my $stat_files;
1706
1707     while (<F>)
1708     {
1709         my ( $parameter, $value ) = split( /\t/, $_ );
1710
1711         my $stat_temp = File::Temp->new(
1712                                       TEMPLATE => "${parameter}_$pop_id-XXXXXX",
1713                                       DIR      => $prod_temp_path,
1714                                       UNLINK   => 0
1715         );
1716         my $stat_file = $stat_temp->filename;
1717
1718         open SF, ">$stat_file" or die "can't open file: !$\n";
1719         print SF $value;
1720         close SF;
1721
1722         $stat_files .= $stat_file . "\t";
1723
1724     }
1725
1726     close F;
1727
1728     my $stat_param_files =
1729       $prod_temp_path . "/" . "stat_temp_files_pop_id_${pop_id}";
1730
1731     open STAT, ">$stat_param_files" or die "can't open file: !$\n";
1732     print STAT $stat_files;
1733     close STAT;
1734
1735     return $stat_param_files;
1736
1737 }
1738
1739 =head2 stat_param_hash
1740
1741  Usage: my %stat_param = $self->stat_param_hash();
1742  Desc: creates a hash (with the statistical parameters (as key) and
1743        their corresponding values) out of a tab delimited
1744        statistical parameters file.
1745  Ret: a hash statistics file
1746  Args: None
1747  Side Effects:
1748  Example:
1749
1750 =cut
1751
1752 sub stat_param_hash
1753 {
1754     my $self           = shift;
1755     my $pop_id         = $self->get_object_id();
1756     my $pop            = $self->get_object();
1757     my $sp_person_id   = $pop->get_sp_person_id();
1758     my $qtl            = CXGN::Phenome::Qtl->new($sp_person_id);
1759     my $c = SGN::Context->new();
1760     my $user_stat_file = $qtl->get_stat_file($c, $pop_id);
1761
1762     open F, "<$user_stat_file" or die "can't open file: !$\n";
1763
1764     my %stat_param;
1765
1766     while (<F>)
1767     {
1768         my ( $parameter, $value ) = split( /\t/, $_ );
1769
1770         $stat_param{$parameter} = $value;
1771
1772     }
1773
1774     return \%stat_param;
1775 }
1776
1777 sub submitter
1778 {
1779     my $self         = shift;
1780     my $population   = $self->get_object();
1781     my $sp_person_id = $population->get_sp_person_id();
1782     my $submitter    = CXGN::People::Person->new( $self->get_dbh(),
1783                                               $population->get_sp_person_id() );
1784     my $submitter_name =
1785       $submitter->get_first_name() . " " . $submitter->get_last_name();
1786     my $submitter_link =
1787 qq |<a href="/solpeople/personal-info.pl?sp_person_id=$sp_person_id">$submitter_name</a> |;
1788
1789     return $submitter, $submitter_link;
1790
1791 }
1792
1793 #move to qtl or population object
1794 sub legend {
1795     my $self = shift;
1796     my $pop = shift;
1797     my $sp_person_id   = $pop->get_sp_person_id();
1798     my $qtl            = CXGN::Phenome::Qtl->new($sp_person_id);
1799     my $stat_file = $qtl->get_stat_file($c, $pop->get_population_id());
1800     my @stat;
1801     my $ci;
1802
1803     open $_, "<", $stat_file or die "$! reading $stat_file\n";
1804     while (my $row = <$_>)
1805     {
1806         my ( $parameter, $value ) = split( /\t/, $row );
1807         if ($parameter =~/qtl_method/) {$parameter = 'Mapping method';}
1808         if ($parameter =~/qtl_model/) {$parameter = 'Mapping model';}
1809         if ($parameter =~/prob_method/) {$parameter = 'QTL genotype probability method';}
1810         if ($parameter =~/step_size/) {$parameter = 'Genome scan size (cM)';}
1811         if ($parameter =~/permu_level/) {$parameter = 'Permutation significance level';}
1812         if ($parameter =~/permu_test/) {$parameter = 'No. of permutations';}
1813         if ($parameter =~/prob_level/) {$parameter = 'QTL genotype signifance level';}
1814
1815         if ($value eq 'zero' || $value eq 'Marker Regression') {$ci = 'none';}
1816
1817         unless (($parameter=~/no_draws/ && $value==' ') ||
1818                ($parameter =~/QTL genotype probability/ && $value==' ')
1819                )
1820         {
1821             push @stat, [map{$_} ($parameter, $value)];
1822         }
1823
1824     }
1825
1826     my $sm;
1827     foreach my $st (@stat) {
1828         foreach my $i (@$st) {
1829             if ($i =~/zero/)  {
1830                 foreach my $s (@stat) {
1831                     foreach my $j (@$s) {
1832                         $j =~ s/Maximum Likelihood/Marker Regression/;
1833                         $ci = 'none';
1834                     }
1835                 }
1836             }
1837         }
1838     }
1839
1840
1841 my $permu_threshold_ref = $self->permu_values();
1842     my %permu_threshold     = %$permu_threshold_ref;
1843
1844     my @keys;
1845
1846     foreach my $key ( keys %permu_threshold )
1847     {
1848         if ( $key =~ m/^\d./ )
1849         {
1850             push @keys, $key;
1851         }
1852
1853     }
1854     my $lod1 = $permu_threshold{ $keys[0] };
1855     my $lod2 = $permu_threshold{ $keys[1] };
1856
1857     if  (!$lod1)
1858     {
1859         $lod1 = qq |<i>Not calculated</i>|;
1860     }
1861
1862     push @stat,
1863     [
1864      map {$_} ('LOD threshold', $lod1)
1865     ];
1866
1867
1868     unless ($ci) {
1869         push @stat,
1870         [
1871          map {$_} ('Confidence interval', 'Based on 95% Bayesian Credible Interval')
1872         ];
1873     }
1874     push @stat,
1875     [
1876      map {$_} ('QTL software', "<a href=http://www.rqtl.org>R/QTL</a>")
1877     ];
1878     my $legend_data = columnar_table_html (
1879                                            headings    => [
1880                                                           ' ',
1881                                                           ' ',
1882
1883                                                          ],
1884                                            data        => \@stat,
1885                                           __alt_freq   => 2,
1886                                           __alt_width  => 1,
1887                                           __alt_offset => 3,
1888                                           __align      => 'l',
1889
1890                                           );
1891
1892
1893     return $legend_data;
1894
1895 }