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 = $population->genotype_file($c);
 993     my $phe_dataset_file = $population->phenotype_file($c);
 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 $basepath     = $c->get_conf("basepath");
1126     my $tempfile_dir = $c->get_conf("tempfiles_subdir");
1127     my $prod_temp_path = $c->get_conf('r_qtl_temp_path');
1128
1129     mkdir $prod_temp_path;
1130     my $prod_cache_path = "$prod_temp_path/cache";
1131     mkdir $prod_cache_path;
1132     $prod_temp_path = "$prod_temp_path/tempfiles";
1133     mkdir $prod_temp_path;
1134     -d $prod_temp_path
1135       or die "temp dir '$prod_temp_path' not found, and could not create!";
1136     -r $prod_temp_path or die "temp dir '$prod_temp_path' not readable!";
1137     -w $prod_temp_path or die "temp dir '$prod_temp_path' not writable!";
1138
1139     my $tempimages_path =
1140       File::Spec->catfile( $basepath, $tempfile_dir, "temp_images" );
1141
1142     return $prod_cache_path, $prod_temp_path, $tempimages_path;
1143
1144 }
1145
1146
1147
1148 =head2 crosstype_file
1149
1150  Usage: my $cross_file = $self->crosstype_file();
1151  Desc: creates the crosstype file in the /data/prod/tmp/r_qtl/temp,
1152
1153  Ret: crosstype filename (with absolute path)
1154  Args: none
1155  Side Effects:
1156  Example:
1157
1158 =cut
1159
1160 sub crosstype_file
1161 {
1162     my $self       = shift;
1163     my $pop_id     = $self->get_object_id();
1164     my $population = $self->get_object();
1165
1166     my $cross_type = 'bc' if ($population->get_cross_type_id() == 2);
1167     $cross_type = 'f2' if ($population->get_cross_type_id() == 1);
1168
1169     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1170         $self->cache_temp_path();
1171
1172     my $cross_temp = File::Temp->new(
1173                                      TEMPLATE => "cross_type_${pop_id}-XXXXXX",
1174                                      DIR      => $prod_temp_path,
1175                                      UNLINK   => 0,
1176                                     );
1177
1178
1179     my $cross_file = $cross_temp->filename;
1180
1181     open CF, ">$cross_file" or die "can't open $cross_file: $!\n";
1182     print CF $cross_type;
1183     close FO;
1184
1185     return $cross_file;
1186
1187 }
1188
1189
1190
1191 =head2 run_r
1192
1193  Usage: my ($qtl_summary, $flanking_markers) = $self->run_r();
1194  Desc: run R in the cluster; copies permutation file from the /data/prod..
1195        to the tempimages dir; returns the R output files (with abosulate filepath) with qtl mapping data
1196        and flanking markers
1197  Ret:
1198  Args:  none
1199  Side Effects:
1200  Example:
1201
1202 =cut
1203
1204 sub run_r
1205 {
1206     my $self = shift;
1207
1208     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1209       $self->cache_temp_path();
1210     my $prod_permu_file = $self->permu_file();
1211     my $file_in         = $self->infile_list();
1212     my ( $file_out, $qtl_summary, $flanking_markers ) = $self->outfile_list();
1213     my $stat_file = $self->stat_files();
1214
1215     CXGN::Tools::Run->temp_base($prod_temp_path);
1216
1217     my ( $r_in_temp, $r_out_temp ) =
1218       map {
1219         my ( undef, $filename ) =
1220           tempfile(
1221                     File::Spec->catfile(
1222                                          CXGN::Tools::Run->temp_base(),
1223                                          "population_indls.pl-$_-XXXXXX",
1224                                        ),
1225                   );
1226         $filename
1227       } qw / in out /;
1228
1229     #copy our R commands into a cluster-accessible tempfile
1230     my $doc = CXGN::Scrap::AjaxPage->new();
1231
1232     {
1233         my $r_cmd_file = $doc->path_to('/cgi-bin/phenome/cvterm_qtl.r');
1234         copy( $r_cmd_file, $r_in_temp )
1235           or die "could not copy '$r_cmd_file' to '$r_in_temp'";
1236     }
1237
1238     # now run the R job on the cluster
1239     my $r_process = CXGN::Tools::Run->run_cluster(
1240         'R', 'CMD', 'BATCH',
1241         '--slave',
1242         "--args $file_in $file_out $stat_file",
1243         $r_in_temp,
1244         $r_out_temp,
1245         {
1246            working_dir => $prod_temp_path,
1247
1248            # don't block and wait if the cluster looks full
1249            max_cluster_jobs => 1_000_000_000,
1250         },
1251     );
1252
1253     sleep 1 while $r_process->alive;    #< wait for R to finish
1254                                         #unlink( $r_in_temp, $r_out_temp );
1255
1256     copy( $prod_permu_file, $tempimages_path )
1257       or die "could not copy '$prod_permu_file' to '$tempimages_path'";
1258
1259     return $qtl_summary, $flanking_markers;
1260
1261 }
1262
1263 =head2 permu_file
1264
1265  Usage: my $permu_file = $self->permu_file();
1266  Desc: creates the permutation file in the /data/prod/tmp/r_qtl/cache,
1267        if it does not exist yet and caches it for R.
1268  Ret: permutation filename (with abosolute path)
1269  Args: none
1270  Side Effects:
1271  Example:
1272
1273 =cut
1274
1275 sub permu_file
1276 {
1277     my $self = shift;
1278     my $doc  = CXGN::Scrap::AjaxPage->new();
1279     my ( $pop_id, $cvterm_id ) =
1280       $doc->get_encoded_arguments( "population_id", "cvterm_id" );
1281
1282     my $dbh = CXGN::DB::Connection->new();
1283
1284     my $population = CXGN::Phenome::Population->new( $dbh, $pop_id );
1285     my $pop_name = $population->get_name();
1286
1287     my ( $term_obj, $term_name, $term_id );
1288
1289     if ( $population->get_web_uploaded() )
1290     {
1291         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $cvterm_id );
1292         $term_name = $term_obj->get_name();
1293         $term_id   = $term_obj->get_user_trait_id();
1294     }
1295     else
1296     {
1297         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $cvterm_id );
1298         $term_name = $term_obj->get_cvterm_name();
1299         $term_id   = $term_obj->get_cvterm_id();
1300     }
1301
1302     my $ac = $population->cvterm_acronym($term_name);
1303
1304     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1305       $self->cache_temp_path();
1306
1307     my $file_cache = Cache::File->new( cache_root => $prod_cache_path );
1308
1309     my $key_permu = "$ac" . "_" . $pop_id . "_permu";
1310     my $filename  = "permu_" . $ac . "_" . $pop_id;
1311
1312     my $permu_file = $file_cache->get($key_permu);
1313
1314     unless ($permu_file)
1315     {
1316
1317         my $permu = undef;
1318
1319         my $permu_file = "$prod_cache_path/$filename";
1320
1321         open OUT, ">$permu_file" or die "can't open $permu_file: !$\n";
1322         print OUT $permu;
1323         close OUT;
1324
1325         $file_cache->set( $key_permu, $permu_file, '30 days' );
1326         $permu_file = $file_cache->get($key_permu);
1327     }
1328
1329     return $permu_file;
1330
1331 }
1332
1333 =head2 permu_values
1334
1335  Usage: my $permu_values = $self->permu_values();
1336  Desc: reads the permutation output from R,
1337        creates a hash with the probality level as key and LOD threshold as the value,
1338
1339  Ret: a hash ref of the permutation values
1340  Args: none
1341  Side Effects:
1342  Example:
1343
1344 =cut
1345
1346 sub permu_values
1347 {
1348     my $self            = shift;
1349     my $prod_permu_file = $self->permu_file();
1350
1351     my %permu_threshold = {};
1352
1353     my $permu_file = fileparse($prod_permu_file);
1354     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1355       $self->cache_temp_path();
1356     $permu_file = File::Spec->catfile( $tempimages_path, $permu_file );
1357
1358     my $round1 = Math::Round::Var->new(0.1);
1359
1360     open PERMUTATION, "<$permu_file"
1361       or die "can't open $permu_file: !$\n";
1362
1363     my $header = <PERMUTATION>;
1364
1365     while ( my $row = <PERMUTATION> )
1366     {
1367         my ( $significance, $lod_threshold ) = split( /\t/, $row );
1368         $lod_threshold = $round1->round($lod_threshold);
1369         $permu_threshold{$significance} = $lod_threshold;
1370     }
1371
1372     close PERMUTATION;
1373
1374     return \%permu_threshold;
1375
1376 }
1377
1378 =head2 permu_values_exist
1379
1380  Usage: my $permu_value = $self->permu_values_exist();
1381  Desc: checks if there is permutation value in the permutation file.
1382  Ret: undef or some value
1383  Args: none
1384  Side Effects:
1385  Example:
1386
1387 =cut
1388
1389 sub permu_values_exist
1390 {
1391     my $self            = shift;
1392     my $prod_permu_file = $self->permu_file();
1393
1394     my ( $size, $permu_file, $permu_data, $tempimages_path, $prod_cache_path,
1395          $prod_temp_path );
1396
1397     if ($prod_permu_file)
1398     {
1399
1400         $permu_file = fileparse($prod_permu_file);
1401         ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1402           $self->cache_temp_path();
1403     }
1404
1405     if ($permu_file)
1406     {
1407
1408         $permu_file = File::Spec->catfile( $tempimages_path, $permu_file );
1409     }
1410
1411     if ( -e $permu_file )
1412     {
1413
1414         open P, "<$permu_file" or die "can't open $permu_file: !$\n";
1415         my $h = <P>;
1416         while ( $permu_data = <P> )
1417         {
1418             last if ($permu_data);
1419
1420             #       #just checking if there is data in there
1421         }
1422         close P;
1423     }
1424
1425     if ($permu_data)
1426     {
1427         return 1;
1428     }
1429     else
1430     {
1431         return 0;
1432
1433     }
1434
1435 }
1436
1437 =head2 qtl_images_exist
1438
1439  Usage: my $qtl_images_ref = $self->qtl_images_exist();
1440  Desc: checks and returns a scalar ref if the qtl plots (with thickbox and their links to the comparative viewer) exist in the cache
1441  Ret: scalar ref to the images or undef
1442  Args: none
1443  Side Effects:
1444  Example:
1445
1446 =cut
1447
1448 sub qtl_images_exist
1449 {
1450     my $self = shift;
1451     my $doc  = CXGN::Scrap::AjaxPage->new();
1452
1453     my ( $pop_id, $cvterm_id ) =
1454       $doc->get_encoded_arguments( "population_id", "cvterm_id" );
1455
1456     my $dbh = $self->get_dbh();
1457
1458     my $population = $self->get_object();
1459     my $pop_name   = $population->get_name();
1460
1461     my @linkage_groups = $population->linkage_groups();
1462     @linkage_groups = sort ( { $a <=> $b } @linkage_groups );
1463
1464     my ( $term_obj, $term_name, $term_id );
1465
1466     if ( $population->get_web_uploaded() )
1467     {
1468         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $cvterm_id );
1469         $term_name = $term_obj->get_name();
1470         $term_id   = $term_obj->get_user_trait_id();
1471     }
1472     else
1473     {
1474         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $cvterm_id );
1475         $term_name = $term_obj->get_cvterm_name();
1476         $term_id   = $term_obj->get_cvterm_id();
1477     }
1478
1479     my $ac = $population->cvterm_acronym($term_name);
1480
1481     my $vh           = SGN::Context->new();
1482     my $basepath     = $vh->get_conf("basepath");
1483     my $tempfile_dir = $vh->get_conf("tempfiles_subdir");
1484
1485     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1486       $self->cache_temp_path();
1487
1488     my $cache_tempimages = Cache::File->new( cache_root => $tempimages_path );
1489     $cache_tempimages->purge();
1490
1491     my ( $qtl_image, $image, $image_t, $image_url, $image_html, $image_t_url,
1492          $thickbox, $title );
1493
1494
1495   IMAGES: foreach my $lg (@linkage_groups)
1496     {
1497         my $cache_qtl_plot = CXGN::Tools::WebImageCache->new();
1498         $cache_qtl_plot->set_basedir($basepath);
1499         $cache_qtl_plot->set_temp_dir( $tempfile_dir . "/temp_images" );
1500
1501         my $key = "qtlplot" . $lg . "small" . $pop_id . $term_id;
1502         $cache_qtl_plot->set_key($key);
1503
1504         my $key_h_marker = "$ac" . "_pop_" . "$pop_id" . "_chr_" . $lg;
1505         my $h_marker     = $cache_tempimages->get($key_h_marker);
1506
1507         if ( $cache_qtl_plot->is_valid )
1508         {
1509             $image      = $cache_qtl_plot->get_image_tag();
1510             $image_url  = $cache_qtl_plot->get_image_url();
1511         }
1512
1513         my $cache_qtl_plot_t = CXGN::Tools::WebImageCache->new();
1514         $cache_qtl_plot_t->set_basedir($basepath);
1515         $cache_qtl_plot_t->set_temp_dir( $tempfile_dir . "/temp_images" );
1516
1517         my $key_t = "qtlplot_" . $lg . "_thickbox_" . $pop_id . $term_id;
1518         $cache_qtl_plot_t->set_key($key_t);
1519
1520         if ( $cache_qtl_plot_t->is_valid )
1521         {
1522
1523             $image_t     = $cache_qtl_plot_t->get_image_tag();
1524             $image_t_url = $cache_qtl_plot_t->get_image_url();
1525
1526             $thickbox =
1527 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> |;
1528
1529             $qtl_image .= $thickbox;
1530             $title = "  ";
1531
1532
1533         }
1534         else
1535         {
1536             $qtl_image = undef;
1537             last IMAGES;
1538
1539         }
1540
1541     }
1542
1543     return $qtl_image;
1544
1545 }
1546
1547 # =head2 user_stat_file
1548
1549 #  Usage:
1550 #  Desc:
1551 #  Ret:
1552 #  Args:
1553 #  Side Effects:
1554 #  Example:
1555
1556 # =cut
1557
1558 # sub user_stat_file {
1559 #     my $self = shift;
1560 #     my $pop = $self->get_object();
1561 #     my $pop_id = $self->get_object_id();
1562 #     my $sp_person_id = $pop->get_sp_person_id();
1563 #     my $qtl  = CXGN::Phenome::Qtl->new($sp_person_id);
1564 #     #$qtl->set_population_id($pop_id);
1565
1566 #     my ($qtl_dir, $user_dir)  = $qtl->get_user_qtl_dir();
1567
1568 #     my $stat_file = "$user_dir/user_stat_pop_$pop_id.txt";
1569 #     print STDERR "stat_file: $stat_file";
1570
1571 #     if (-e $stat_file) {
1572 #       return $stat_file;
1573 #     } else {return 0;}
1574
1575 # }
1576
1577 =head2 stat_files
1578
1579  Usage: my $stat_param_files = $self->stat_files();
1580  Desc:  creates a master file containing individual files
1581         in /data/prod/tmp/r_qtl for each statistical parameter
1582         which are feed to  R.
1583  Ret:  an absolute path to the statistical parameter's
1584        master file (and individual files)
1585  Args: None
1586  Side Effects:
1587  Example:
1588
1589 =cut
1590
1591 sub stat_files
1592 {
1593     my $self           = shift;
1594     my $pop_id         = $self->get_object_id();
1595     my $pop            = $self->get_object();
1596     my $sp_person_id   = $pop->get_sp_person_id();
1597     my $qtl            = CXGN::Phenome::Qtl->new($sp_person_id);
1598     my $c = SGN::Context->new();
1599     my $user_stat_file = $qtl->get_stat_file($c, $pop_id);
1600
1601     my ( $prod_cache_path, $prod_temp_path, $tempimages_path ) =
1602       $self->cache_temp_path();
1603
1604     open F, "<$user_stat_file" or die "can't open file: !$\n";
1605
1606     my $stat_files;
1607
1608     while (<F>)
1609     {
1610         my ( $parameter, $value ) = split( /\t/, $_ );
1611
1612         my $stat_temp = File::Temp->new(
1613                                       TEMPLATE => "${parameter}_$pop_id-XXXXXX",
1614                                       DIR      => $prod_temp_path,
1615                                       UNLINK   => 0
1616         );
1617         my $stat_file = $stat_temp->filename;
1618
1619         open SF, ">$stat_file" or die "can't open file: !$\n";
1620         print SF $value;
1621         close SF;
1622
1623         $stat_files .= $stat_file . "\t";
1624
1625     }
1626
1627     close F;
1628
1629     my $stat_param_files =
1630       $prod_temp_path . "/" . "stat_temp_files_pop_id_${pop_id}";
1631
1632     open STAT, ">$stat_param_files" or die "can't open file: !$\n";
1633     print STAT $stat_files;
1634     close STAT;
1635
1636     return $stat_param_files;
1637
1638 }
1639
1640 =head2 stat_param_hash
1641
1642  Usage: my %stat_param = $self->stat_param_hash();
1643  Desc: creates a hash (with the statistical parameters (as key) and
1644        their corresponding values) out of a tab delimited
1645        statistical parameters file.
1646  Ret: a hash statistics file
1647  Args: None
1648  Side Effects:
1649  Example:
1650
1651 =cut
1652
1653 sub stat_param_hash
1654 {
1655     my $self           = shift;
1656     my $pop_id         = $self->get_object_id();
1657     my $pop            = $self->get_object();
1658     my $sp_person_id   = $pop->get_sp_person_id();
1659     my $qtl            = CXGN::Phenome::Qtl->new($sp_person_id);
1660     my $c = SGN::Context->new();
1661     my $user_stat_file = $qtl->get_stat_file($c, $pop_id);
1662
1663     open F, "<$user_stat_file" or die "can't open file: !$\n";
1664
1665     my %stat_param;
1666
1667     while (<F>)
1668     {
1669         my ( $parameter, $value ) = split( /\t/, $_ );
1670
1671         $stat_param{$parameter} = $value;
1672
1673     }
1674
1675     return \%stat_param;
1676 }
1677
1678 sub submitter
1679 {
1680     my $self         = shift;
1681     my $population   = $self->get_object();
1682     my $sp_person_id = $population->get_sp_person_id();
1683     my $submitter    = CXGN::People::Person->new( $self->get_dbh(),
1684                                               $population->get_sp_person_id() );
1685     my $submitter_name =
1686       $submitter->get_first_name() . " " . $submitter->get_last_name();
1687     my $submitter_link =
1688 qq |<a href="/solpeople/personal-info.pl?sp_person_id=$sp_person_id">$submitter_name</a> |;
1689
1690     return $submitter, $submitter_link;
1691
1692 }
1693
1694 #move to qtl or population object
1695 sub legend {
1696     my $self = shift;
1697     my $pop = shift;
1698     my $sp_person_id   = $pop->get_sp_person_id();
1699     my $qtl            = CXGN::Phenome::Qtl->new($sp_person_id);
1700     my $stat_file = $qtl->get_stat_file($c, $pop->get_population_id());
1701     my @stat;
1702     my $ci;
1703
1704     open my $sf, "<", $stat_file or die "$! reading $stat_file\n";
1705     while (my $row = <$sf>)
1706     {
1707         chomp($row);
1708         my ( $parameter, $value ) = split( /\t/, $row );
1709         if ($parameter =~/qtl_method/) {$parameter = 'Mapping method';}
1710         if ($parameter =~/qtl_model/) {$parameter = 'Mapping model';}
1711         if ($parameter =~/prob_method/) {$parameter = 'QTL genotype probability method';}
1712         if ($parameter =~/step_size/) {$parameter = 'Genome scan size (cM)';}
1713         if ($parameter =~/permu_level/) {$parameter = 'Permutation significance level';}
1714         if ($parameter =~/permu_test/) {$parameter = 'No. of permutations';}
1715         if ($parameter =~/prob_level/) {$parameter = 'QTL genotype signifance level';}
1716         if ($parameter =~/stat_no_draws/) {$parameter = 'No. of imputations';}
1717
1718         if ($value eq 'zero' || $value eq 'Marker Regression') {$ci = 'none';}
1719
1720         unless (($parameter =~/No. of imputations/ && !$value ) ||
1721                 ($parameter =~/QTL genotype probability/ && !$value ) ||
1722                 ($parameter =~/Permutation significance level/ && !$value)
1723                )
1724
1725         {
1726             push @stat, [map{$_} ($parameter, $value)];
1727
1728         }
1729
1730     }
1731
1732     my $sm;
1733     foreach my $st (@stat) {
1734         foreach my $i (@$st) {
1735             if ($i =~/zero/)  {
1736                 foreach my $s (@stat) {
1737                     foreach my $j (@$s) {
1738                         $j =~ s/Maximum Likelihood/Marker Regression/;
1739                         $ci = 'none';
1740                     }
1741                 }
1742             }
1743         }
1744     }
1745
1746
1747 my $permu_threshold_ref = $self->permu_values();
1748     my %permu_threshold     = %$permu_threshold_ref;
1749
1750     my @keys;
1751
1752     foreach my $key ( keys %permu_threshold )
1753     {
1754         if ( $key =~ m/^\d./ )
1755         {
1756             push @keys, $key;
1757         }
1758
1759     }
1760     my $lod1 = $permu_threshold{ $keys[0] };
1761    # my $lod2 = $permu_threshold{ $keys[1] };
1762
1763     if  (!$lod1)
1764     {
1765         $lod1 = qq |<i>Not calculated</i>|;
1766     }
1767
1768     push @stat,
1769     [
1770      map {$_} ('LOD threshold', $lod1)
1771     ];
1772
1773
1774     unless ($ci) {
1775         push @stat,
1776         [
1777          map {$_} ('Confidence interval', 'Based on 95% Bayesian Credible Interval')
1778         ];
1779     }
1780     push @stat,
1781     [
1782      map {$_} ('QTL software', "<a href=http://www.rqtl.org>R/QTL</a>")
1783     ];
1784     my $legend_data = columnar_table_html (
1785                                            headings    => [
1786                                                           ' ',
1787                                                           ' ',
1788
1789                                                          ],
1790                                            data        => \@stat,
1791                                           __alt_freq   => 2,
1792                                           __alt_width  => 1,
1793                                           __alt_offset => 3,
1794                                           __align      => 'l',
1795
1796                                           );
1797
1798
1799     return $legend_data;
1800
1801 }