cgi-bin/phenome/qtl.pl

   1
   2 =head1 DESCRIPTION
   3 A QTL detail page.
   4
   5 =head1 AUTHOR
   6
   7 Isaak Y Tecle (iyt2@cornell.edu)
   8
   9 =cut
  10
  11 use strict;
  12 use warnings;
  13
  14
  15 use CGI qw //;
  16 use CXGN::People::PageComment;
  17 use CXGN::Phenome::Population;
  18 use CXGN::Phenome::UserTrait;
  19 use CXGN::Phenome::Qtl;
  20 use CXGN::Marker;
  21 use CXGN::Map;
  22 use CXGN::DB::Connection;
  23 use CXGN::Chado::Cvterm;
  24 use List::MoreUtils qw / uniq /;
  25 use List::Util qw / max min/;
  26 use File::Slurp qw / read_file /;
  27 use Number::Format;
  28 use SGN::Exception;
  29
  30
  31 use CatalystX::GlobalContext qw( $c );
  32
  33 my $cgi = CGI->new();
  34 my %params = $cgi->Vars();
  35
  36 our $pop_id    = $params{population_id};
  37 our $trait_id  = $params{term_id};
  38 our $lg        = $params{chr};
  39 our $p_m       = $params{peak_marker};
  40 our $lod       = $params{lod};
  41 our $qtl_image = $params{qtl};
  42
  43 if (    !$pop_id
  44      || !$trait_id
  45      || !$lg
  46      || !$p_m
  47      || !$qtl_image )
  48 {
  49     die 'QTL detail page error:  A required argument is missing';
  50 }
  51
  52
  53 our ($l_m,  $r_m);
  54 our $dbh         = CXGN::DB::Connection->new();
  55 our $pop         = CXGN::Phenome::Population->new($dbh, $pop_id);
  56 our $pop_name    = $pop->get_name();
  57 our $trait_name  = trait_name();
  58 our ($ci_table, $marker_details) = confidence_interval();
  59
  60 foreach my $k (keys %$marker_details) {
  61     $l_m = $k if ($marker_details->{$k}{orientation} eq  'left');
  62     $r_m = $k if ($marker_details->{$k}{orientation} eq  'right');
  63 }
  64
  65 my $genetic_link = genetic_map();
  66 my $cmv_link     = marker_positions();
  67 my $markers      = markers();
  68 my $legend       = legend();
  69 my $comment      = comment();
  70 my $download_qtl = download_qtl_region();
  71 $ci_table        = order_by_position();
  72
  73 $c->forward_to_mason_view( '/qtl/qtl/index.mas',
  74                            qtl_image    => $qtl_image,
  75                            pop_name     => $pop_name,
  76                            trait_name   => $trait_name,
  77                            cmv_link     => $cmv_link,
  78                            markers      => $markers,
  79                            marker_link  => $ci_table,
  80                            genetic_map  => $genetic_link,
  81                            legend       => $legend,
  82                            download     => $download_qtl,
  83                            comment      => $comment,
  84 );
  85
  86
  87
  88 =head2 marker_positions
  89
  90  Usage: $map_viewer = marker_positions();
  91  Desc: generates a link to the comparative map viewer page
  92        using the flanking markers and peak marker.
  93  Ret: a link to the map viewer page
  94  Args: None
  95  Side Effects:
  96  Example:
  97
  98 =cut
  99
 100
 101 sub marker_positions
 102 {
 103     my $mapv_id = $pop->mapversion_id();
 104     my $l_m_pos = $marker_details->{$l_m}{position};
 105     my $p_m_pos = $marker_details->{$p_m}{position};
 106     my $r_m_pos = $marker_details->{$r_m}{position};
 107
 108     my $fl_markers
 109         = qq |<a href="../cview/view_chromosome.pl?map_version_id=$mapv_id&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_m_pos&cM_end=$r_m_pos">Chromosome $lg ($l_m, $r_m)</a> |;
 110
 111     return $fl_markers;
 112 }
 113
 114 =head2 markers
 115
 116  Usage: $markers = markers();
 117  Desc: creates  marker objects
 118  Ret:  array ref of marker objects
 119  Args: None
 120  Side Effects:
 121  Example:
 122
 123 =cut
 124
 125
 126 sub markers
 127 {
 128     my @mrs = uniq ($l_m, $p_m, $r_m);
 129
 130     my @markers;
 131     foreach my $mr (@mrs)
 132     {
 133         if($mr)
 134         {
 135             my $marker = CXGN::Marker->new_with_name($dbh, $mr);
 136             if (!$marker)
 137             {
 138                 my @marker_id = CXGN::Marker::Tools::marker_name_to_ids( $dbh, $mr);
 139                 if ($marker_id[0])
 140                 {
 141                     $marker = CXGN::Marker->new($dbh, $marker_id[0]);
 142                 }
 143             }
 144
 145             push @markers, $marker if $marker;
 146         }
 147     }
 148
 149     return \@markers;
 150
 151 }
 152
 153 =head2 genetic_map
 154
 155  Usage: $population_map = genetic_map();
 156  Desc: generates a link to the genetic map of the
 157        population
 158  Ret: a link to the genetic map
 159  Args: None
 160  Side Effects:
 161  Example:
 162
 163 =cut
 164
 165
 166 sub genetic_map
 167 {
 168     my $mapv_id     = $pop->mapversion_id();
 169     my $map         = CXGN::Map->new( $dbh, { map_version_id => $mapv_id } );
 170     my $map_name    = $map->get_long_name();
 171     my $map_sh_name = $map->get_short_name();
 172     my $genetic_map
 173         = qq | <a href=/cview/map.pl?map_version_id=$mapv_id&hilite=$l_m+$p_m+$r_m>$map_name ($map_sh_name)</a>|;
 174
 175     return $genetic_map;
 176
 177 }
 178
 179 =head2 confidence_interval
 180
 181  Usage: ($marker_table, $marker_details) = confidence_interval();
 182  Desc: reads the confidence interval data for the QTL from a
 183        file containing the genome-wide confidence intervals
 184        and their lod profile. It calculates the left and right
 185        markers, their position values; and some interpretation of
 186        the data etc..
 187  Ret: an array ref of marker details table (for the viewer)
 188       and a ref to a hash of hash for the  marker details
 189       (for later access to the data)
 190  Args:
 191  Side Effects:
 192  Example:
 193
 194 =cut
 195
 196
 197
 198 sub confidence_interval
 199 {
 200     my $ci_lod_file = $pop->ci_lod_file( $c,
 201                                         $pop->cvterm_acronym( $trait_name )
 202         );
 203
 204     my (@marker_lods,  @all_lods, @all_positions, @marker_html);
 205     my %marker_details_of = ();
 206
 207     my @rows =  grep { /\t$lg\t/ } read_file( $ci_lod_file );
 208
 209     my $rnd  = Number::Format->new();
 210
 211     foreach my $row (@rows) {
 212         my ( $m, $m_chr, $m_pos, $m_lod ) = split (/\t/, $row);
 213         push @all_lods, $m_lod;
 214         push @all_positions, $m_pos;
 215
 216         my $marker = CXGN::Marker->new_with_name( $dbh, $m );
 217
 218         unless  ( !$marker )
 219         {
 220             push @marker_lods, $m_lod;
 221         }
 222     }
 223
 224     my $peak_marker_lod = $rnd->round(max( @marker_lods), 2 );
 225     my $highest_lod     = $rnd->round(max( @all_lods), 2 );
 226     my $right_position  = $rnd->round(max( @all_positions), 2 );
 227     my $left_position   = $rnd->round(min( @all_positions), 2 );
 228
 229     my ($peak_marker, $linkage_group, $peak_position) = split (/\t/, $rows[1]);
 230     $peak_position = $rnd->round($peak_position, 1);
 231
 232     foreach my $row ( @rows )
 233     {
 234         my ($m, $m_chr, $m_pos, $m_lod)  = split (/\t/, $row);
 235         $m_pos = $rnd->round( $m_pos, 2 );
 236         $m_lod = $rnd->round( $m_lod, 2 );
 237
 238         my $marker = CXGN::Marker->new_with_name( $dbh, $m );
 239
 240         unless ( !$marker )
 241         {
 242
 243             $marker_details_of{$m}{name}          = $m;
 244             $marker_details_of{$m}{linkage_group} = $m_chr;
 245
 246             $marker_details_of{$m}{position}  = !$m_pos && $m_pos ne '' ? '0.0'
 247                                               : $m_pos eq '' ? 'NA'
 248                                               : $m_pos
 249                                               ;
 250
 251             $marker_details_of{$m}{lod_score} = $m_lod;
 252
 253             if ($m eq $p_m) { $marker_details_of{$m}{orientation} = 'peak'; }
 254             if ($m_pos == $right_position) { $marker_details_of{$m}{orientation} = 'right'; }
 255             if ($m_pos == $left_position) { $marker_details_of{$m}{orientation} = 'left'; }
 256
 257             my $m_id    = $marker->marker_id();
 258             my $remark1 = "<i>Highest LOD score is $highest_lod at $peak_position cM</i>."  if $m_lod == $peak_marker_lod;
 259             my $remark2 = "<i>The closest marker to the peak position ($peak_position cM)</i>."  if $m eq $p_m;
 260
 261             push @marker_html,
 262             [
 263              map { $_ } (
 264                  qq | <a href="/search/markers/markerinfo.pl?marker_id=$m_id">$m</a>|,
 265                  $marker_details_of{$m}{position},
 266                  $m_lod,
 267                  $remark1 . $remark2,
 268              )
 269             ];
 270         }
 271     }
 272
 273
 274     return \@marker_html, \%marker_details_of;
 275
 276 }
 277
 278 =head2 trait_name
 279
 280  Usage: $trait_name = trait_name()
 281  Desc: returns the name of the QTL trait
 282  Ret: trait name
 283  Args: None
 284  Side Effects:
 285  Example:
 286
 287 =cut
 288
 289
 290 sub trait_name
 291 {
 292     my ( $term_obj, $term_name, $term_id );
 293     if ( $pop->get_web_uploaded() )
 294     {
 295         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $trait_id );
 296         $term_name = $term_obj->get_name();
 297         $term_id   = $term_obj->get_user_trait_id();
 298     }
 299     else
 300     {
 301         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $trait_id );
 302         $term_name = $term_obj->get_cvterm_name();
 303         $term_id   = $term_obj->get_cvterm_id();
 304     }
 305
 306     return $term_name;
 307 }
 308
 309 =head2 legend
 310
 311  Usage: $legend = legend();
 312  Desc: generates the appropriate legend describing the
 313        statistical methods and parameters used for the
 314        QTL analysis
 315  Ret: an array ref
 316  Args: None
 317  Side Effects:
 318  Example:
 319
 320 =cut
 321
 322
 323 sub legend
 324 {
 325     my $user_id;
 326     if ($c->user) {
 327         $user_id = $c->user->get_object->get_sp_person_id;
 328     } else {
 329         $user_id = $pop->get_sp_person_id();
 330     }
 331
 332     my $qtl            = CXGN::Phenome::Qtl->new( $user_id );
 333     my $user_stat_file = $qtl->get_stat_file( $c, $pop_id );
 334     my @stat;
 335     my $ci;
 336
 337     open my $uf, "<", $user_stat_file or die "$! reading $user_stat_file";
 338     while ( my $row = <$uf> )
 339     {
 340         chomp($row);
 341         my ( $parameter, $value ) = split( /\t/, $row );
 342
 343         if ( $parameter =~ /qtl_method/ )
 344         {
 345             $parameter = 'Mapping method';
 346         }
 347         if ( $parameter =~ /qtl_model/ )
 348         {
 349             $parameter = 'Mapping model';
 350         }
 351         if ( $parameter =~ /prob_method/ )
 352         {
 353             $parameter = 'QTL genotype probability method';
 354         }
 355         if ( $parameter =~ /step_size/ )
 356         {
 357             $parameter = 'Genome scan size (cM)';
 358         }
 359         if ( $parameter =~ /permu_level/ )
 360         {
 361             $parameter = 'Permutation significance level';
 362         }
 363         if ( $parameter =~ /permu_test/ ) {
 364             $parameter = 'No. of permutations';
 365         }
 366         if ( $parameter =~ /prob_level/ )
 367         {
 368             $parameter = 'QTL genotype significance level';
 369         }
 370         if ( $parameter =~ /stat_no_draws/ )
 371         {
 372             $parameter = 'No. of imputations';
 373         }
 374         if ( $value eq 'zero' || $value eq 'Marker Regression' )
 375         {
 376             $ci = 'none';
 377         }
 378
 379         unless ( ( $parameter =~ /No. of imputations/ && !$value )
 380             || ( $parameter =~ /QTL genotype probability/ && !$value )
 381             || ( $parameter =~ /Permutation significance level/ && !$value ) )
 382
 383         {
 384             push @stat, [ $parameter, $value ];
 385
 386         }
 387     }
 388
 389     foreach my $st ( @stat )
 390     {
 391         foreach my $i ( @$st )
 392         {
 393             if ( $i =~ /zero/ )
 394             {
 395                 foreach my $s ( @stat )
 396                 {
 397                     foreach my $j (@$s)
 398                     {
 399                         $j =~ s/Maximum Likelihood/Marker Regression/;
 400                         $ci = 'none';
 401                     }
 402                 }
 403             }
 404         }
 405     }
 406
 407     if ( !$lod )
 408     {
 409         $lod = qq |<i>Not calculated</i>|;
 410     }
 411
 412     push @stat, [ map {$_} ( 'LOD threshold', $lod ) ];
 413
 414     unless ($ci)
 415     {
 416         push @stat,
 417             [ map {$_} ( 'Confidence interval',
 418                          'Based on 95% Bayesian Credible Interval'
 419               )
 420             ];
 421     }
 422
 423     return \@stat;
 424
 425 }
 426
 427
 428 sub download_qtl_region
 429 {
 430 my $link = qq | <a href="https://www.eu-sol.wur.nl/marker2seq/marker2seq.do?marker1=$l_m&marker2=$r_m">View/download</a> genetic markers in the tomato F2.2000 reference map region (+5cM) matching the QTL region and gene models from the ITAG annotated tomato genome. |;
 431
 432 $link  .=   qq | <p><i>Courtesy of</i> <a href="http://www.eu-sol.wur.nl"><img src ="/img/eusol_logo_small.jpg"/></a></p> |;
 433
 434     return $link;
 435
 436 }
 437
 438
 439 =head2 comment
 440
 441  Usage: $comment = comment();
 442  Desc: generates the comment section html
 443  Ret: the comment html
 444  Args:
 445  Side Effects:
 446  Example:
 447
 448 =cut
 449
 450
 451 sub comment
 452 {
 453     my $comment;
 454     if ($pop_id)
 455     {
 456         my $page_comment_obj =
 457             CXGN::People::PageComment->new( $dbh, "population", $pop_id,
 458                                     "/phenome/qtl.pl?population_id=$pop_id" );
 459         $comment = $page_comment_obj->get_html();
 460     }
 461     return $comment;
 462
 463 }
 464
 465
 466 sub order_by_position {
 467     my ($marker_html, $markers_details) = confidence_interval();
 468     my @marker_html = sort { $a->[1] <=> $b->[1] }  @$marker_html;
 469     return \@marker_html;
 470 }