improve pca score plot accessions trial color coding

R/solGS/pca.r

@@ -1,8 +1,8 @@
  #SNOPSIS
 
- #runs population structure analysis using PCA from SNPRelate, a bioconductor R package
+ #runs population structure analysis using PCA on genotype or phenotype data
 
  #AUTHOR
  # Isaak Y Tecle (iyt2@cornell.edu)
 
 
@@ -17,16 +17,16 @@
 library(phenoAnalysis)
 library(ggplot2)
 
 allArgs <- commandArgs()
 
 outputFile  <- grep("output_files", allArgs, value = TRUE)
 outputFiles <- scan(outputFile, what = "character")
 
 inputFile  <- grep("input_files", allArgs, value = TRUE)
 inputFiles <- scan(inputFile, what = "character")
 
 scoresFile       <- grep("pca_scores", outputFiles, value = TRUE)
 screeDataFile    <- grep("pca_scree_data", outputFiles, value = TRUE)
 screeFile        <- grep("pca_scree_plot", outputFiles, value = TRUE)
 loadingsFile     <- grep("pca_loadings", outputFiles, value = TRUE)
 varianceFile     <- grep("pca_variance", outputFiles, value = TRUE)
@@ -60,13 +60,16 @@
 dataType <- ifelse(isTRUE(pcaDataFile[1]), 'genotype', 'phenotype')
 
 if (dataType == 'genotype') {
-    if (length(inputFiles) > 1) {
-        genoData <- genoDataFilter::combineGenoData(inputFiles)
+    genoFiles <- grep("genotype_data", inputFiles, value = TRUE)
+    trialMembershipFile <- grep("pca_trial_membership", inputFiles, value = TRUE)
+
+    if (length(genoFiles) > 1) {
+        genoData <- genoDataFilter::combineGenoData(genoFiles)
         genoMetaData   <- genoData$trial
         genoData$trial <- NULL
 
     } else {
-        genoDataFile <- grep("genotype_data", inputFiles,  value = TRUE)
+        genoDataFile <- genoFiles
         genoData     <- fread(genoDataFile,
                               header = TRUE,
                               na.strings = c("NA", " ", "--", "-", "."))
@@ -85,6 +88,14 @@
         genoData <- data.frame(genoData)
         genoData <- column_to_rownames(genoData, 'V1')
 
+        if (length(trialMembershipFile) == 1) {
+            trialMembership <- fread(trialMembershipFile, header = TRUE)
+            genoMetaData <- trialMembership$trial[
+                match(rownames(genoData), trialMembership$accession)
+            ]
+            genoMetaData[is.na(genoMetaData)] <- "unassigned"
+        }
+
     }
 } else if (dataType == 'phenotype') {
 

js/source/legacy/solGS/gebvs.js

@@ -132,10 +132,8 @@ jQuery(document).ready(function () {
                 var gebvsDownloadLinks = solGS.gebvs.createGebvsDownloadLinks(res);
                 var geneticValuesDownloadLinks = solGS.gebvs.createGeneticValuesDownloadLinks(res);
                 var downloadLinks = `${gebvsDownloadLinks} | ${geneticValuesDownloadLinks}`;
-                console.log(`Calling getGebvsData`)
 
                 solGS.gebvs.getGebvsData().done(function (res) {
-                    console.log(`getGebvsData res: ${JSON.stringify(res)}`)
                     solGS.gebvs.plotGebvs(res.gebvs_data, downloadLinks);
                 });
             });

js/source/legacy/solGS/heatMap.js

@@ -49,8 +49,6 @@ solGS.heatmap = {
       jQuery("#" + scatterPlotDivId).before("<div id=" + scatterPlotMsgDivId + "></div>");
     }
 
-
-
     var parsedScatterInputData = null;
     if (scatterInputData) {
         parsedScatterInputData = JSON.parse(scatterInputData);
@@ -65,11 +63,11 @@ solGS.heatmap = {
     }
 
     var nLabels = parsedHeatmaprInputData.labels.length;
+    var fontSize = "0.85em";
 
     if (nLabels >= 100) {
       height = 600;
       width = 600;
-      fontSize = "0.75em";
     } else {
       height = 500;
       width = 500;
@@ -89,8 +87,6 @@ solGS.heatmap = {
     var nral = "#98AFC7"; //blue gray
     var txtColor = "#523CB5";
 
-    var fontSize = "0.95em";
-
     var corr = [];
     var coefs = [];
 
@@ -172,7 +168,7 @@ solGS.heatmap = {
       .attr("id", heatmapPlotDivId)
       .attr("transform", "translate(0, 0)");
 
-    corrplot
+    yAxisLabels = corrplot
       .append("g")
       .attr("class", "y_axis")
       .attr("transform", `translate(${pad.left}, ${pad.top})`)
@@ -184,7 +180,7 @@ solGS.heatmap = {
       .attr("fill", txtColor)
       .style("font-size", fontSize);
 
-    corrplot
+    var xAxisLabels = corrplot
       .append("g")
       .attr("class", "x_axis")
       .attr("transform", `translate(${pad.left}, ${pad.top + height})`)
@@ -197,6 +193,18 @@ solGS.heatmap = {
       .attr("transform", "rotate(-90)")
       .attr("fill", txtColor)
       .style("font-size", fontSize);
+
+    var longestXAxisLabelWidth = 0;
+    xAxisLabels.each(function () {
+      longestXAxisLabelWidth = Math.max(
+        longestXAxisLabelWidth,
+        this.getComputedTextLength()
+      );
+    });
+
+    pad.bottom = Math.max(pad.bottom, Math.ceil(longestXAxisLabelWidth + 25));
+    totalH = height + pad.top + pad.bottom;
+    svg.attr("height", totalH);
 
     corrplot
       .selectAll()
@@ -402,7 +410,22 @@ solGS.heatmap = {
       if (!heatmapPlotDivId.match("#")) {
         heatmapPlotDivId = "#" + heatmapPlotDivId;
       }
-      jQuery(heatmapPlotDivId).append('<p style="margin: 20px 0px 0px 40px">' + downloadLinks + "</p>");
+
+    var longestYAxisLabelWidth = 0;
+    yAxisLabels.each(function () {
+      longestYAxisLabelWidth = Math.max(
+        longestYAxisLabelWidth,
+        this.getComputedTextLength()
+      );
+    });
+
+    var yAxisLabelsLeft = Math.max(
+      0,
+      Math.ceil(pad.left - 10 - longestYAxisLabelWidth)
+    );
+      jQuery(heatmapPlotDivId).append(
+        `<p style="margin: 20px 0 0 ${yAxisLabelsLeft}px">${downloadLinks}</p>`
+      );
     }
   },
 

js/source/legacy/solGS/pca.js

@@ -202,6 +202,11 @@ solGS.pca = {
         if (selectedPopDiv) {
             var selectedPopData = selectedPopDiv.dataset;
             var selectedPop = JSON.parse(selectedPopData.selectedPop);
+
+            if (runPcaElemId.match(/save_pcs/)) {
+                return selectedPop;
+            }
+
             pcaPopId = selectedPop.pca_pop_id;
 
             var pcaArgs = selectedPopData.selectedPop;
@@ -932,6 +937,10 @@ solGS.pca = {
             popName = plotData.list_name;
         }
 
+        if (plotData.pca_pop_name) {
+            popName = plotData.pca_pop_name;
+        }
+
         popName = popName
             ? popName + " (" + plotData.data_type + ")"
             : " (" + plotData.data_type + ")";
@@ -972,7 +981,7 @@ solGS.pca = {
 
                     groupName = "common to: " + groupName.join(", ");
                 } else {
-                    groupName = trialsNames[id];
+                    groupName = trialsNames[id] || id;
                 }
                 legendValues.push([cnt, id, groupName]);
                 cnt++;

lib/SGN/Controller/solGS/pca.pm

@@ -160,7 +160,7 @@ sub prepare_pca_output_response {
 
             if ($scores) {
                 $res = {
-                    "scores"          =>  $scores,
+                    "scores"          => $scores,
                     "variances"       => $variances,
                     "scores_file"     => $scores_file,
                     "variances_file"  => $variances_file,
@@ -171,6 +171,7 @@ sub prepare_pca_output_response {
                     "status"          => 'success',
                     "cached"          => 1,
                     "pca_pop_id"      => $c->stash->{pca_pop_id},
+                    "pca_pop_name"    => $c->stash->{pca_pop_name},
                     "file_id"         => $file_id,
                     "list_id"         => $c->stash->{list_id},
                     "trials_names"    => $trials_names,
@@ -434,21 +435,50 @@ sub pca_input_files {
 
 }
 
+sub pca_trial_membership_file {
+    my ( $self, $c ) = @_;
+
+    my $dataset_id = $c->stash->{dataset_id};
+
+    my $rows = $c->controller('solGS::Search')->model($c)
+      ->get_dataset_accession_trial_memberships($dataset_id);
+
+    if (!@$rows) {
+        return;
+    }
+
+    my $file_id = $c->stash->{file_id};
+    my $name = "pca_trial_membership_${file_id}";
+    my $tmp_dir = $self->pca_temp_dir($c);
+    my $file = $c->controller('solGS::Files')->create_tempfile($tmp_dir, $name);
+
+    my $headers = "accession" . "\t" . "accession_id" . "\t" . "trial" . "\n";
+    my @lines = ( $headers );
+    push @lines, map { join( "\t", @$_ ) . "\n" } @$rows;
+    write_file( $file, { binmode => ':utf8' }, @lines );
+
+    $c->stash->{pca_trial_membership_file} = $file;
+
+    return $file;
+}
+
 sub pca_geno_input_files {
     my ( $self, $c ) = @_;
 
     my $data_type = $c->stash->{data_type};
     my $files = [];
 
     if ( $data_type =~ /genotype/i ) {
-        if ( $c->req->referer =~
-            /solgs\/selection\/|solgs\/combined\/model\/\d+\/selection\// )
-        {
+        if ( $c->req->referer =~ /solgs\/selection\/|solgs\/combined\/model\/\d+\/selection\// ) {
             $self->training_selection_geno_files($c);
         }
 
-        $files =
-          $c->stash->{genotype_files_list} || $c->stash->{genotype_file_name};
+        $files = $c->stash->{genotype_files_list} || $c->stash->{genotype_file_name};
+
+        if ( $c->stash->{data_structure} && $c->stash->{data_structure} =~ /dataset/ ) {
+            my $membership_file = $self->pca_trial_membership_file($c);
+            $files .= "\t$membership_file" if $membership_file;
+        }
     }
 
     $files = join( "\t", @$files ) if reftype($files) eq 'ARRAY';

lib/SGN/Model/solGS/solGS.pm

@@ -2128,6 +2128,52 @@ sub get_dataset_data {
     return $dataset_data;
 }
 
+sub get_dataset_accession_trial_memberships {
+    my ( $self, $dataset_id ) = @_;
+
+    my $data          = $self->get_dataset_data($dataset_id);
+    my $categories    = $data->{categories} || {};
+    my $accession_ids = $categories->{accessions} || [];
+    my $trial_ids     = $categories->{trials} || [];
+
+    if (!@$accession_ids || !@$trial_ids) {
+        return [];
+    }
+
+    my $accession_placeholders = join( ', ', ('?') x @$accession_ids );
+    my $trial_placeholders     = join( ', ', ('?') x @$trial_ids );
+    my $query = qq{
+        SELECT accession.stock_id, accession.uniquename, axt.trial_id
+        FROM stock AS accession
+        LEFT JOIN accessionsXtrials AS axt
+          ON axt.accession_id = accession.stock_id
+         AND axt.trial_id IN ($trial_placeholders)
+        WHERE accession.stock_id IN ($accession_placeholders)
+        ORDER BY accession.stock_id, axt.trial_id
+    };
+
+    my $sth = $self->schema->storage->dbh->prepare($query);
+    $sth->execute( @$trial_ids, @$accession_ids );
+
+    my %memberships;
+    while ( my ( $accession_id, $accession_name, $trial_id ) =
+        $sth->fetchrow_array() )
+    {
+        $memberships{$accession_id}{name} = $accession_name;
+        push @{ $memberships{$accession_id}{trials} }, $trial_id
+          if defined $trial_id;
+    }
+
+    my @rows;
+    foreach my $accession_id ( sort { $a <=> $b } keys %memberships ) {
+        my @trials = uniq( @{ $memberships{$accession_id}{trials} || [] } );
+        my $trial_group = @trials ? join( '-', sort { $a <=> $b } @trials ) : 'unassigned';
+        push @rows, [ $memberships{$accession_id}{name}, $accession_id, $trial_group ];
+    }
+
+    return \@rows;
+}
+
 sub get_dataset_plots_list {
     my ( $self, $dataset_id ) = @_;
 

t/selenium2/solgs/pca.t

@@ -21,53 +21,39 @@ my $solgs_data = SGN::Test::solGSData->new({
 
 my $cache_dir = $solgs_data->base_analyses_cache_dir();
 my $pca_dir = catdir( $cache_dir, 'pca' );
-print STDERR "\ncache_dir-- $cache_dir\n";
 
 my $accessions_list = $solgs_data->load_accessions_list();
-
-# my $accessions_list = $solgs_data->get_list_details('accessions');
 my $accessions_list_name = $accessions_list->{list_name};
 my $accessions_list_id   = 'list_' . $accessions_list->{list_id};
-print STDERR "\naccessions list: $accessions_list_name -- $accessions_list_id\n";
-my $plots_list = $solgs_data->load_plots_list();
 
-# my $plots_list =  $solgs_data->get_list_details('plots');
+my $plots_list = $solgs_data->load_plots_list();
 my $plots_list_name = $plots_list->{list_name};
 my $plots_list_id   = 'list_' . $plots_list->{list_id};
 
-print STDERR "\nadding trials list\n";
 my $trials_list = $solgs_data->load_trials_list();
-
-# my $trials_list =  $solgs_data->get_list_details('trials');
 my $trials_list_name = $trials_list->{list_name};
 my $trials_list_id   = 'list_' . $trials_list->{list_id};
-print STDERR "\nadding trials dataset\n";
 
-# my $trials_dt =  $solgs_data->get_dataset_details('trials');
 my $trials_dt      = $solgs_data->load_trials_dataset();
 my $trials_dt_name = $trials_dt->{dataset_name};
 my $trials_dt_id   = 'dataset_' . $trials_dt->{dataset_id};
-print STDERR "\nadding accessions dataset\n";
 
-# my $accessions_dt =  $solgs_data->get_dataset_details('accessions');
 my $accessions_dt      = $solgs_data->load_accessions_dataset();
 my $accessions_dt_name = $accessions_dt->{dataset_name};
 my $accessions_dt_id   = 'dataset_' . $accessions_dt->{dataset_id};
 
-print STDERR "\nadding plots dataset\n";
-
-# my $plots_dt =  $solgs_data->get_dataset_details('plots');
 my $plots_dt      = $solgs_data->load_plots_dataset();
 my $plots_dt_name = $plots_dt->{dataset_name};
 my $plots_dt_id   = 'dataset_' . $plots_dt->{dataset_id};
 
-
 my @test_trials_ids = @{$solgs_data->trials_ids()};
 
 remove_tree($cache_dir, {safe => 1});
 
 $d->while_logged_in_as("submitter", sub {
 
+
+    ###### standalone pca analysis page #####
     $d->get_ok('/pca/analysis', 'pca home page');
     sleep(5);
 
@@ -283,6 +269,9 @@ $d->while_logged_in_as("submitter", sub {
     remove_tree($cache_dir, {safe => 1});
     sleep(5);
 
+
+    ########## trial detail page pca analysis ##########
+
     $d->get_ok('/breeders/trial/' . $test_trials_ids[0], 'trial detail home page');
     sleep(10);
     my $analysis_tools = $d->find_element('Analysis Tools', 'partial_link_text', 'toogle analysis tools');
@@ -363,6 +352,8 @@ $d->while_logged_in_as("submitter", sub {
     remove_tree($cache_dir, {safe => 1});
     sleep(5);
 
+
+    ####### solgs pca analysis ########
     $d->get_ok('/solgs', 'solgs homepage');
     sleep(10);