Phenotype Data

Here we’ll load the phenotype data, and take a look at the first six lines. Each row corresponds to a grapevine accession from the USDA germplasm repository. Although 1200 accesssions were phenotyped by Chitwood et al (2012), only 122 had descriptions using Galet’s method. The first 13 columns are the digital traits, and the last 10 are the classical measurements.

chitwood = read.csv("Phenotypes/Chitwood_Galet.csv", row.names = 1)

# head(chitwood)

Spearman correlation between digital and classical ampelography metrics

Here, we’ll look at the correlation between digital and classical amelography measurements using Spearman’s Rank correlation.

sp.cor = cor(chitwood[1:13], chitwood[14:23], method = "spearman")

Plot the results

Now, we’ll plot the results as a nice heatmap.

# install.packages(gplots)
library(gplots)
heatmap.2(sp.cor, Rowv = F, Colv = F, dendrogram = "none", trace = "none", cellnote = round(sp.cor, 
    2), notecol = "black", notecex = 0.5, breaks = seq(-1, 1, length.out = 100), 
    margins = c(6, 6))

Hyperspectral data

The imaging system provides a stack of 235 greyscale images for each plant. Each image in the stack is a specific wavelength. For each image, we remove the background objects leaving just the “plant pixels”. Each plant pixel will be a single chandanel and the intensity of the channel should be proportional to the reflectance. So, a higher intensity means more reflectance for a given wavelength. For each image we obtain the mean pixel intensity (sum the pixel intensities and divide by the total number of plant pixels). This is equivalent to the mean reflectance for each plant. Here, we will read in the file of the mean pixel intensities for each plant.

dat = read.csv("Phenotypes/Hyperspec_pixel_density.csv")

# Check the first 10 columns of the phenotypic data head(dat[1:10])

# Plot the refectance for all wavelengths
dat.means = apply(dat[4:ncol(dat)], 2, mean)
plot(seq(from = 585, to = 1755, by = 5), dat.means, type = "l", ylab = "Reflectance", 
    xlab = "Wavelength")

Here, I will calculate the normalized difference vegetation index (NDVI). I will create a new dataframe named ‘NDVI’ that has the NDVI calculations, the plant ID, the treatment, and the day of imaging. I’ll create a plot of NDVI over the six time points (1 - 22 days of imaging).

# Make the data frame
NDVI.df = data.frame(ID = dat$ID, Treatment = dat$Treatment, DayOfImaging = dat$DayOfImaging, 
    NDVI = NDVI(dat))

# boxplot(NDVI~Treatment, data=NDVI.df, ylab='NDVI')

# Plot results install.packages('ggplot2')
library(ggplot2)

ggplot(NDVI.df, aes(x = DayOfImaging, y = NDVI, colour = Treatment)) + geom_smooth(method = "loess") + 
    xlab("Day of Imaging") + ylab("NDVI")

The indices defined above only capture a small portion of the total phenotype described by the hyperspectral data. Here, we will first utilize the whole hyperspectral range to predict plant water content. This data was collected at 24 days after imaging. For this part of the excersize we will utilize partial least squares regression in the package ‘pls’. Below I will provide a brief explaination of the methodology. For those interested, the ‘pls’ manual provides an in-depth explaination of the PLSR technique, and a few example datasets PLS.

# install.packages('pls')
library(pls)
# This dataset contains the reflacance values for all the wavelengths. These
# will be used as predictors for the percent of total biomass that is water
# ('Percent.Water').
dat2 = read.csv("Phenotypes/Hyperspec_PWC.csv")
dat2$ID = NULL
dat2$Treatment = NULL
dat2$TotalWater = NULL
# Make a training set and validation set to assess the prediction accuracy
# This is similar to what we did for genomic prediction, however here we
# will only sample the data once.
traindat <- dat2[1:50, ]
testdat <- dat2[51:87, ]

# Now we'll fit the PLSR model. We're regressing the Percent water on the
# hyperspectral data (Percent.Water ~ a).
full.spec <- plsr(Percent.Water ~ ., data = traindat)
summary(full.spec)

Data:   X dimension: 50 242 
    Y dimension: 50 1
Fit method: kernelpls
Number of components considered: 49
TRAINING: % variance explained
               1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
X                97.32    98.97    99.26    99.33    99.38    99.43
               7 comps  8 comps  9 comps  10 comps  11 comps  12 comps
X                99.50    99.55    99.58     99.61     99.64     99.67
               13 comps  14 comps  15 comps  16 comps  17 comps  18 comps
X                 99.70     99.71     99.74     99.76     99.80     99.81
               19 comps  20 comps  21 comps  22 comps  23 comps  24 comps
X                 99.82     99.84     99.85     99.87     99.88     99.89
               25 comps  26 comps  27 comps  28 comps  29 comps  30 comps
X                  99.9     99.91     99.92     99.92     99.93     99.94
               31 comps  32 comps  33 comps  34 comps  35 comps  36 comps
X                 99.94     99.95     99.95     99.96     99.97     99.97
               37 comps  38 comps  39 comps  40 comps  41 comps  42 comps
X                 99.97     99.98     99.98     99.98     99.98     99.99
               43 comps  44 comps  45 comps  46 comps  47 comps  48 comps
X                 99.99     99.99     99.99       100       100       100
               49 comps
X                   100
 [ reached getOption("max.print") -- omitted 1 row ]

plot(RMSEP(full.spec), legendpos = "topright")

# This line shows the how well the model fits our data
plot(full.spec, ncomp = 20, asp = 1, line = TRUE)

# Now we predict the percent water for the testing dataset using the model
# fit using the training dataset
Y.tst = predict(full.spec, ncomp = 20, newdata = testdat)

Plot the results

plot(Y.tst, testdat$Percent.Water, xlab = "Percent Water Predicted", ylab = "Percent Water Observed")
abline(0, 1, col = "red")

Calculate the correlation between predicted values and observed values in the data set

cor(Y.tst, testdat$Percent.Water)

[1] 0.2602748

sub.hyperspec <- function(startBW, endBW) {
    if (sum(colnames(dat2) == paste0("bw_", startBW)) > 0) {
        start.col = which(colnames(dat2) == paste0("bw_", startBW))
    } else {
        print("Column name for first bandwidth does not exist")
    }
    if (sum(colnames(dat2) == paste0("bw_", endBW)) > 0) {
        end.col = which(colnames(dat2) == paste0("bw_", endBW))
    } else {
        print("Column name for last bandwidth does not exist")
    }
    tmp.df = as.data.frame(cbind(dat2$Percent.Water, dat2[, start.col:end.col]))
    colnames(tmp.df)[1] = "Percent.Water"
    return(tmp.df)
}

dat3 = sub.hyperspec(585, 605)

traindat2 <- dat3[1:50, ]
testdat2 <- dat3[51:87, ]

# Now we'll fit the PLSR model. We're regressing the Percent water on the
# hyperspectral data (Percent.Water ~ a).
full.spec <- plsr(Percent.Water ~ ., data = traindat2)
summary(full.spec)

Data:   X dimension: 50 5 
    Y dimension: 50 1
Fit method: kernelpls
Number of components considered: 5
TRAINING: % variance explained
               1 comps  2 comps  3 comps  4 comps  5 comps
X                91.53    94.02    95.49    97.94   100.00
Percent.Water    18.38    24.36    24.63    24.65    24.65

plot(RMSEP(full.spec), legendpos = "topright")

plot(full.spec, asp = 1, line = TRUE)

plot(Y.tst, testdat$Percent.Water, xlab = "Percent Water Predicted", ylab = "Percent Water Observed")
abline(0, 1, col = "red")

Calculate the correlation between predicted values and observed values in the data set

cor(Y.tst, testdat$Percent.Water)

[1] 0.2602748