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         push @markers, CXGN::Marker->new_with_name($dbh, $mr);
 134     }
 135
 136     return \@markers;
 137
 138 }
 139
 140 =head2 genetic_map
 141
 142  Usage: $population_map = genetic_map();
 143  Desc: generates a link to the genetic map of the
 144        population
 145  Ret: a link to the genetic map
 146  Args: None
 147  Side Effects:
 148  Example:
 149
 150 =cut
 151
 152
 153 sub genetic_map
 154 {
 155     my $mapv_id     = $pop->mapversion_id();
 156     my $map         = CXGN::Map->new( $dbh, { map_version_id => $mapv_id } );
 157     my $map_name    = $map->get_long_name();
 158     my $map_sh_name = $map->get_short_name();
 159     my $genetic_map
 160         = qq | <a href=/cview/map.pl?map_version_id=$mapv_id&hilite=$l_m+$p_m+$r_m>$map_name ($map_sh_name)</a>|;
 161
 162     return $genetic_map;
 163
 164 }
 165
 166 =head2 confidence_interval
 167
 168  Usage: ($marker_table, $marker_details) = confidence_interval();
 169  Desc: reads the confidence interval data for the QTL from a
 170        file containing the genome-wide confidence intervals
 171        and their lod profile. It calculates the left and right
 172        markers, their position values; and some interpretation of
 173        the data etc..
 174  Ret: an array ref of marker details table (for the viewer)
 175       and a ref to a hash of hash for the  marker details
 176       (for later access to the data)
 177  Args:
 178  Side Effects:
 179  Example:
 180
 181 =cut
 182
 183
 184
 185 sub confidence_interval
 186 {
 187     my $ci_lod_file = $pop->ci_lod_file( $c,
 188                                         $pop->cvterm_acronym( $trait_name )
 189         );
 190
 191     my (@marker_lods,  @all_lods, @all_positions, @marker_html);
 192     my %marker_details_of = ();
 193
 194     my @rows =  grep { /\t$lg\t/ } read_file( $ci_lod_file );
 195
 196     my $rnd  = Number::Format->new();
 197
 198     foreach my $row (@rows) {
 199         my ( $m, $m_chr, $m_pos, $m_lod ) = split (/\t/, $row);
 200         push @all_lods, $m_lod;
 201         push @all_positions, $m_pos;
 202
 203         my $marker = CXGN::Marker->new_with_name( $dbh, $m );
 204
 205         unless  ( !$marker )
 206         {
 207             push @marker_lods, $m_lod;
 208         }
 209     }
 210
 211     my $peak_marker_lod = $rnd->round(max( @marker_lods), 2 );
 212     my $highest_lod     = $rnd->round(max( @all_lods), 2 );
 213     my $right_position  = $rnd->round(max( @all_positions), 2 );
 214     my $left_position   = $rnd->round(min( @all_positions), 2 );
 215
 216     my ($peak_marker, $linkage_group, $peak_position) = split (/\t/, $rows[1]);
 217     $peak_position = $rnd->round($peak_position, 1);
 218
 219     foreach my $row ( @rows )
 220     {
 221         my ($m, $m_chr, $m_pos, $m_lod)  = split (/\t/, $row);
 222         $m_pos = $rnd->round( $m_pos, 1 );
 223         $m_lod = $rnd->round( $m_lod, 2 );
 224
 225         my $marker = CXGN::Marker->new_with_name( $dbh, $m );
 226
 227         unless ( !$marker )
 228         {
 229
 230             $marker_details_of{$m}{name}          = $m;
 231             $marker_details_of{$m}{linkage_group} = $m_chr;
 232
 233             $marker_details_of{$m}{position}  = !$m_pos && $m_pos ne '' ? '0.0'
 234                                               : $m_pos eq '' ? 'NA'
 235                                               : $m_pos
 236                                               ;
 237
 238             $marker_details_of{$m}{lod_score} = $m_lod;
 239
 240             if ($m eq $p_m) { $marker_details_of{$m}{orientation} = 'peak'; }
 241             if ($m_pos == $right_position) { $marker_details_of{$m}{orientation} = 'right'; }
 242             if ($m_pos == $left_position) { $marker_details_of{$m}{orientation} = 'left'; }
 243
 244             my $m_id    = $marker->marker_id();
 245             my $remark1 = "<i>Highest LOD score is $highest_lod at $peak_position cM</i>."  if $m_lod == $peak_marker_lod;
 246             my $remark2 = "<i>The closest marker to the peak position ($peak_position cM)</i>."  if $m eq $p_m;
 247
 248             push @marker_html,
 249             [
 250              map { $_ } (
 251                  qq | <a href="/search/markers/markerinfo.pl?marker_id=$m_id">$m</a>|,
 252                  $marker_details_of{$m}{position},
 253                  $m_lod,
 254                  $remark1 . $remark2,
 255              )
 256             ];
 257         }
 258     }
 259
 260
 261     return \@marker_html, \%marker_details_of;
 262
 263 }
 264
 265 =head2 trait_name
 266
 267  Usage: $trait_name = trait_name()
 268  Desc: returns the name of the QTL trait
 269  Ret: trait name
 270  Args: None
 271  Side Effects:
 272  Example:
 273
 274 =cut
 275
 276
 277 sub trait_name
 278 {
 279     my ( $term_obj, $term_name, $term_id );
 280     if ( $pop->get_web_uploaded() )
 281     {
 282         $term_obj  = CXGN::Phenome::UserTrait->new( $dbh, $trait_id );
 283         $term_name = $term_obj->get_name();
 284         $term_id   = $term_obj->get_user_trait_id();
 285     }
 286     else
 287     {
 288         $term_obj  = CXGN::Chado::Cvterm->new( $dbh, $trait_id );
 289         $term_name = $term_obj->get_cvterm_name();
 290         $term_id   = $term_obj->get_cvterm_id();
 291     }
 292
 293     return $term_name;
 294 }
 295
 296 =head2 legend
 297
 298  Usage: $legend = legend();
 299  Desc: generates the appropriate legend describing the
 300        statistical methods and parameters used for the
 301        QTL analysis
 302  Ret: an array ref
 303  Args: None
 304  Side Effects:
 305  Example:
 306
 307 =cut
 308
 309
 310 sub legend
 311 {
 312     my $user_id;
 313     if ($c->user) {
 314         $user_id = $c->user->get_object->get_sp_person_id;
 315     } else {
 316         $user_id = $pop->get_sp_person_id();
 317     }
 318
 319     my $qtl            = CXGN::Phenome::Qtl->new( $user_id );
 320     my $user_stat_file = $qtl->get_stat_file( $c, $pop_id );
 321     my @stat;
 322     my $ci;
 323
 324     open my $uf, "<", $user_stat_file or die "$! reading $user_stat_file";
 325     while ( my $row = <$uf> )
 326     {
 327         chomp($row);
 328         my ( $parameter, $value ) = split( /\t/, $row );
 329
 330         if ( $parameter =~ /qtl_method/ )
 331         {
 332             $parameter = 'Mapping method';
 333         }
 334         if ( $parameter =~ /qtl_model/ )
 335         {
 336             $parameter = 'Mapping model';
 337         }
 338         if ( $parameter =~ /prob_method/ )
 339         {
 340             $parameter = 'QTL genotype probability method';
 341         }
 342         if ( $parameter =~ /step_size/ )
 343         {
 344             $parameter = 'Genome scan size (cM)';
 345         }
 346         if ( $parameter =~ /permu_level/ )
 347         {
 348             $parameter = 'Permutation significance level';
 349         }
 350         if ( $parameter =~ /permu_test/ ) {
 351             $parameter = 'No. of permutations';
 352         }
 353         if ( $parameter =~ /prob_level/ )
 354         {
 355             $parameter = 'QTL genotype significance level';
 356         }
 357         if ( $parameter =~ /stat_no_draws/ )
 358         {
 359             $parameter = 'No. of imputations';
 360         }
 361         if ( $value eq 'zero' || $value eq 'Marker Regression' )
 362         {
 363             $ci = 'none';
 364         }
 365
 366         unless ( ( $parameter =~ /No. of imputations/ && !$value )
 367             || ( $parameter =~ /QTL genotype probability/ && !$value )
 368             || ( $parameter =~ /Permutation significance level/ && !$value ) )
 369
 370         {
 371             push @stat, [ $parameter, $value ];
 372
 373         }
 374     }
 375
 376     foreach my $st ( @stat )
 377     {
 378         foreach my $i ( @$st )
 379         {
 380             if ( $i =~ /zero/ )
 381             {
 382                 foreach my $s ( @stat )
 383                 {
 384                     foreach my $j (@$s)
 385                     {
 386                         $j =~ s/Maximum Likelihood/Marker Regression/;
 387                         $ci = 'none';
 388                     }
 389                 }
 390             }
 391         }
 392     }
 393
 394     if ( !$lod )
 395     {
 396         $lod = qq |<i>Not calculated</i>|;
 397     }
 398
 399     push @stat, [ map {$_} ( 'LOD threshold', $lod ) ];
 400
 401     unless ($ci)
 402     {
 403         push @stat,
 404             [ map {$_} ( 'Confidence interval',
 405                          'Based on 95% Bayesian Credible Interval'
 406               )
 407             ];
 408     }
 409
 410     return \@stat;
 411
 412 }
 413
 414
 415 sub download_qtl_region
 416 {
 417 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. |;
 418
 419 $link  .=   qq | <p><i>Courtesy of</i> <a href="http://www.eu-sol.wur.nl"><img src ="/img/eusol_logo_small.jpg"/></a></p> |;
 420
 421     return $link;
 422
 423 }
 424
 425
 426 =head2 comment
 427
 428  Usage: $comment = comment();
 429  Desc: generates the comment section html
 430  Ret: the comment html
 431  Args:
 432  Side Effects:
 433  Example:
 434
 435 =cut
 436
 437
 438 sub comment
 439 {
 440     my $comment;
 441     if ($pop_id)
 442     {
 443         my $page_comment_obj =
 444             CXGN::People::PageComment->new( $dbh, "population", $pop_id,
 445                                     "/phenome/qtl.pl?population_id=$pop_id" );
 446         $comment = $page_comment_obj->get_html();
 447     }
 448     return $comment;
 449
 450 }
 451
 452
 453 sub order_by_position {
 454     my ($marker_html, $markers_details) = confidence_interval();
 455     my @marker_html = sort { $a->[1] <=> $b->[1] }  @$marker_html;
 456     return \@marker_html;
 457 }