PCAPlot.R

# Amy Campbell 2016-2017
# PCA plotting script for chronobiome pilot study
#
# Fits and plots principal component analysis of
# activity and communication variables
# 
# Plots principal components generated by Seth Rhodes 
# of metabolite, protein, and microbe data

PlotPCA <- function(PCdata, title, filename) {
  # :param: PCdata - matrix of principle components for the dataset
  # :param: Title - string representing the title to be shared by the 
  #                 three side-by-side plots
  # :param: filename - file name, including ".eps" suffix, for the plot
  postscript(filename, width = 480, height = 480)
  plot1 <- qplot(x = PC2, y=PC3, data = PCdata, colour = Subject) +
    scale_colour_manual(values = colorblind_Palette) 
  plot2 <- qplot(x = PC1, y = PC3, data = PCdata, colour = Subject) +
    scale_colour_manual(values = colorblind_Palette) 
  plot3 <- qplot(x = PC1, y = PC2, data = PCdata, colour = Subject) +
    scale_colour_manual(values = colorblind_Palette) 
  grid.arrange(plot1, plot2, plot3, ncol = 3,
               top = title)
  dev.off()
}

ParseSubject <- function(row, half) {
    # Function to be applied to all rows of 
    # the dataframe in question
    # :param: row - the entire row of the dataframe
    # :param: half - indicates whether the "subject" information
    # is in the first or second half of the string
    timesubjectindex <- row["TimeSubjectIndex"]
    strlist <- (strsplit(toString(timesubjectindex), "_"))
    
    if (toString(half) == "2") {
        return(unlist(strlist)[2])
    } else {
        return(as.numeric(unlist(strlist)[1]))
    }
}

TimeList_window1 <- c(seq(922, 1232))
TimeList_window2 <- c(seq(1328, 1545))

TimeList <- c(TimeList_window1, TimeList_window2)


###############
# Load Packages
###############

# Read date of snapshot to use
args = commandArgs(TRUE)
checkpoint_date = args[1]

# Use checkpoint created by install script
library("checkpoint")
checkpoint(snapshotDate=checkpoint_date, checkpointLocation='.')

library("ggplot2")
library("gridExtra")


################################
# Communication and activity PCA
################################

communication.activity <- (read.csv("communication.activity.csv"))
communication.activity['subject'] <- 
  (apply(communication.activity, 1, ParseSubject, 2))

TimeSubj <- data.frame(communication.activity$TimeSubjectIndex)
subj <- data.frame(communication.activity$subject)
length <- as.numeric(dim(communication.activity)[2])
Matrix_For_PCA <- communication.activity[,3:(length-1)]
act.com.PCA <- prcomp(Matrix_For_PCA, scale.= TRUE)
act.com.pca.matrix <- data.frame(act.com.PCA$x)

act.com.pca.matrix["Subject"] <- subj
write.table(act.com.PCA$rotation, "Act_Com_Loadings.csv", sep = ",")


act.com.pca.matrix["TimeSubjectIndex"] <- TimeSubj

act.com.pca.matrix["TimesOnly"] <- 
  apply(act.com.pca.matrix, 1, ParseSubject, 1)

windowed_ActCom_PCA <- 
  subset(act.com.pca.matrix, TimesOnly %in% TimeList)

#Also perform a PCA analysis without the circadian stats from
#Lux, Mobility, and MobilityRadius. The variance correlation
#plots indicate all of these metrics tend to be correlated
#with each other (the scatter plots of the various circadian
#signal metrics show this particularly well). Here I'll repeat
#the PCA analysis and plotting without the additional circadian
#stats from Lux, Mobility, and MobilityRadius for comparison
#purposes.

act.com.PCA.no_circ_stats <- prcomp(dplyr::select(Matrix_For_PCA,
                                                  -dplyr::starts_with("Lux"),
                                                  -dplyr::starts_with("MobilityRadius"),
                                                  -c(Mobility.amplitude, Mobility.period,
                                                     Mobility.phase, circadian.signal.mob,
                                                     circadian.signal.mobR, circadian.signal.lux)),
                                    scale.= TRUE)
act.com.pca.matrix.no_circ_stats <- data.frame(act.com.PCA.no_circ_stats$x)

act.com.pca.matrix.no_circ_stats["Subject"] <- subj
write.table(act.com.PCA.no_circ_stats$rotation, "Act_Com_Loadings.no_circ_stats.csv", sep = ",")


act.com.pca.matrix.no_circ_stats["TimeSubjectIndex"] <- TimeSubj

act.com.pca.matrix.no_circ_stats["TimesOnly"] <- 
    apply(act.com.pca.matrix.no_circ_stats, 1, ParseSubject, 1)

windowed_ActCom_PCA.no_circ_stats <- 
    subset(act.com.pca.matrix.no_circ_stats, TimesOnly %in% TimeList)


###########
# PCA Plots
###########
colorblind_Palette <- c("#000000", "#0072B2", "#56B4E9", "#F0E442",
                        "#D55E00", "#CC79A7")

# Plot activity and communication PC combinations individually
p_1 <- qplot(x=PC2, y=PC3, data=windowed_ActCom_PCA, colour=Subject) +
  scale_colour_manual(values=colorblind_Palette) 
p_2 <- qplot(x=PC1, y=PC3, data=windowed_ActCom_PCA, colour=Subject) +
  scale_colour_manual(values=colorblind_Palette)
p_3 <- qplot(x=PC1, y=PC2, data=windowed_ActCom_PCA, colour=Subject) +
  scale_colour_manual(values=colorblind_Palette)

postscript("Facets_Act_Com_PC2PC3.eps")
fp_1 <- p_1 + facet_grid(. ~Subject)
grid.arrange(fp_1, ncol = 1,
             top = paste0("Principal Components for Combined Activity",
                          "and Communication (Visits 1 and 2)"))
dev.off()

postscript("Facets_Act_Com_PC1PC3.eps")
fp_2 <- p_2 + facet_grid(. ~Subject) 
grid.arrange(fp_2, ncol = 1, 
             top = paste0("Principal Components for Combined Activity and",
                          " Communication (Visits 1 and 2)"))
dev.off()

postscript("Facets_Act_Com_PC1PC2.eps")
fp_3 <- p_3 + facet_grid(. ~Subject) 
grid.arrange(fp_3, ncol = 1, 
             top = paste0("Principal Components for Combined Activity and", 
                          " Communication (Visits 1 and 2)"))
dev.off()


# Activity and Energy PCA plots with PC combinations side-by-side
PlotPCA(act.com.pca.matrix, 
        title = paste0("Principal Components for Combined Activity and ",
                       "Communication (All Measurements)"), 
        filename = "act.com.PCA.PlotAll.eps")

PlotPCA(windowed_ActCom_PCA, 
        title= paste0("Principal Components for Combined Activity and ",
                      "Communication (Visits 1 and 2)"), 
        filename = "act.com.PCAPlot_48.eps")

#Regenerate PCA plots excluding the circadian stats from Lux, Mobility,
#and MobilityRadius

# Plot activity and communication PC combinations individually
p_1 <- qplot(x=PC2, y=PC3, data=windowed_ActCom_PCA.no_circ_stats, colour=Subject) +
    scale_colour_manual(values=colorblind_Palette) 
p_2 <- qplot(x=PC1, y=PC3, data=windowed_ActCom_PCA.no_circ_stats, colour=Subject) +
    scale_colour_manual(values=colorblind_Palette)
p_3 <- qplot(x=PC1, y=PC2, data=windowed_ActCom_PCA.no_circ_stats, colour=Subject) +
    scale_colour_manual(values=colorblind_Palette)

postscript("Facets_Act_Com_PC2PC3.no_circ_stats.eps")
fp_1 <- p_1 + facet_grid(. ~Subject)
grid.arrange(fp_1, ncol = 1,
             top = paste0("Principal Components for Combined Activity",
                          "and Communication (Visits 1 and 2)\n",
                          "Excluding circadian status from lux, mobility,",
                          "and mobility radius"))
dev.off()

postscript("Facets_Act_Com_PC1PC3.no_circ_stats.eps")
fp_2 <- p_2 + facet_grid(. ~Subject) 
grid.arrange(fp_2, ncol = 1, 
             top = paste0("Principal Components for Combined Activity and",
                          " Communication (Visits 1 and 2)\n",
                          "Excluding circadian status from lux, mobility,",
                          "and mobility radius"))
dev.off()

postscript("Facets_Act_Com_PC1PC2.no_circ_stats.eps")
fp_3 <- p_3 + facet_grid(. ~Subject) 
grid.arrange(fp_3, ncol = 1, 
             top = paste0("Principal Components for Combined Activity and", 
                          " Communication (Visits 1 and 2)\n",
                          "Excluding circadian status from lux, mobility,",
                          "and mobility radius"))
dev.off()


# Activity and Energy PCA plots with PC combinations side-by-side
PlotPCA(act.com.pca.matrix.no_circ_stats, 
        title = paste0("Principal Components for Combined Activity and ",
                       "Communication (All Measurements)\n",
                       "Excluding circadian status from lux, mobility,",
                       "and mobility radius"), 
        filename = "act.com.PCA.PlotAll.no_circ_stats.eps")

PlotPCA(windowed_ActCom_PCA.no_circ_stats, 
        title= paste0("Principal Components for Combined Activity and ",
                      "Communication (Visits 1 and 2)\n",
                      "Excluding circadian status from lux, mobility,",
                      "and mobility radius"), 
        filename = "act.com.PCAPlot_48.no_circ_stats.eps")




# Load and plot other PCAs (Generated by S. Rhodes)

#Saliva Metabolites
saliva.metabolites <- read.csv(
  file.path("chronobiome", "Seth_PCA", "PC_Scores",
            "SalivaMetabolites_BothVisits_ScoresFor10PCs.csv"))
colnames(saliva.metabolites) <- c("Visit", "Subject", "Hour", "PC1", "PC2",
                                  "PC3", "PC4", "PC5", "PC6", "PC7",
                                  "PC8", "PC9", "PC10")
PlotPCA(saliva.metabolites, 
        title = "Principal Components for Saliva Metabolites(Visits 1 and 2)",
        filename = "Saliva_Metabs_PCA.eps")

# Saliva Microbes
Saliva_Microbes <-
  read.csv(file.path("chronobiome", "Seth_PCA", "PC_Scores", 
                     "SalivaMicrobes_BothVisits_ScoresFor5PCs.csv"))
colnames(Saliva_Microbes) <-(c("X", "Subject", "Hour","PC1",
                               "PC2", "PC3", "PC4", "PC5"))
PlotPCA(Saliva_Microbes,
        title = "Principal Components for Saliva Microbes(Visits 1 and 2)",
        filename = "Saliva_Microbes_PCA.eps")

# Plasma Proteins
plasma.proteins <- read.csv(
  file.path("chronobiome", "Seth_PCA", "PC_Scores", 
            "PlasmaProteins_1stVisits_HCR009-12hrRemove_ScoresFor10PCs.csv"))
plasma.proteins <- na.omit(plasma.proteins)
colnames(plasma.proteins) <- c("X", "Subject", "Hour", "PC1", "PC2", "PC3",
                               "PC4", "PC5", "PC6", "PC7","PC8","PC9","PC10")
PlotPCA(plasma.proteins,
        title = "Principal Components for Plasma Proteins(Visit 1)",
        filename = "Plasma_Proteins_PCA.eps")

# Plasma Metabolites
plasma.metabolites <- read.csv(
  file.path("chronobiome", "Seth_PCA", "PC_Scores", 
            "PlasmaMetabolites_BothVisits_ScoresFor10PCs.csv"))
colnames(plasma.metabolites) <- c("Visit", "Subject", "Hour", "PC1",
                                  "PC2", "PC3", "PC4", "PC5", "PC6",
                                  "PC7", "PC8", "PC9", "PC10")
PlotPCA(plasma.metabolites,
        title = "Principal Components for Plasma Metabolites(Visits 1 and 2)",
        filename = "Plasma_Metabs_PCA.eps")

# Metabolites and microbes
metabolites.microbiome.both<- 
  read.csv(file.path("chronobiome", "Seth_PCA", "PC_Scores", "ScoresJan17",
                     "MetabsAndMicrobes_BothVisits_ScoresFor5PCs.csv"))
colnames(metabolites.microbiome.both) <- c("X", "Subject", "Hour",
                                           "PC1", "PC2", "PC3",
                                           "PC4","PC5")
PlotPCA(metabolites.microbiome.both,
        title = paste0("Principal Components for Metabolites and Microbes",
                       "Combined (Visits 1 and 2)"), 
        filename = "Metabolites_Microbes_BothVisits.eps"
)

# Metabolites and Microbes, Visit 1
metabolites.microbiome.visit1 <- 
  read.csv(file.path("chronobiome", "Seth_PCA", "PC_Scores", "ScoresJan17",
                     "MetabsAndMicrobes_Visit1_ScoresFor5PCs.csv"))
colnames(metabolites.microbiome.visit1) <-
  c("X", "Subject", "Hour", "PC1", "PC2", "PC3", "PC4", "PC5")
PlotPCA(metabolites.microbiome.visit1,
        title = paste0("Principal Components for Metabolites and Microbes",
                       "Combined (Visit 1)"),
        filename = "Metabolites_Microbes_Visit1.eps"
)

# Metabolites and Microbes, Visit 2
metabolites.microbiome.visit2 <- 
  read.csv(file.path("chronobiome", "Seth_PCA", "PC_Scores","ScoresJan17",
                     "MetabsAndMicrobes_Visit2_ScoresFor5PCs.csv"))
colnames(metabolites.microbiome.visit2) <- c("X", "Subject", "Hour", "PC1",
                                             "PC2", "PC3", "PC4", "PC5")
PlotPCA(metabolites.microbiome.visit2,
        title = paste0("Principal Components for Metabolites",
                       " and Microbes Combined (Visit 2)"),
        filename = "Metabolites_Microbes_Visit2.eps"
)

# Tri-omics PCA on Metabolites, Proteins, Microbes
TriOmic_Visit1 <- read.csv(file.path("chronobiome", "Seth_PCA",
                                     "PC_Scores", "ScoresJan17",
                                     "TriOmic_Visit1_ScoresFor5PCs.csv"))
colnames(TriOmic_Visit1) <- c("X", "Subject", "Hour", "PC1", "PC2",
                              "PC3", "PC4", "PC5")
PlotPCA(TriOmic_Visit1,
        title = paste0("Principal Components for Metabolites,",
                       " Proteins and Microbes Combined (Visit 1)"),
        filename="TriOmic_Visit1.eps")

# Seth's updated Saliva Metabolites PCA excluding an outlier
saliva.metabolites.jan <- 
  read.csv(file.path("chronobiome", "Seth_PCA", "PC_Scores", "ScoresJan17",
                     "SalivaMetabolites_BothVisits_ScoresFor5PCs.csv"))
colnames(saliva.metabolites.jan) <- c("Visit", "Subject", "Hour",
                                      "PC1", "PC2", "PC3", "PC4", "PC5")
PlotPCA(saliva.metabolites.jan,
        title = "Principal Components for Saliva Metabolites(Visits 1 and 2)",
        filename = "Saliva_Metabs_PCA.jan17.eps" )

# Seth's updated Saliva Microbes PCA excluding an outlier
Saliva_Microbes.jan <- 
  read.csv(file.path("chronobiome", "Seth_PCA", "PC_Scores", "ScoresJan17", 
                     "SalivaMicrobes_BothVisits_ScoresFor5PCs.csv"))
colnames(Saliva_Microbes.jan) <- 
  (c("X", "Subject", "Hour", "PC1", "PC2", "PC3", "PC4", "PC5"))
PlotPCA(Saliva_Microbes.jan, 
        title = "Principal Components for Saliva Microbes(Visits 1 and 2)",
        filename = "Saliva_Microbes_PCA_Jan.eps")

# Seth's updated Plasma Proteins PCA excluding an outlier

plasma.proteins.jan <- 
  read.csv(file.path("chronobiome", "Seth_PCA", "PC_Scores", 
                     "ScoresJan17", "PlasmaProteins_Visit1_ScoresFor5PCs.csv"))
plasma.proteins.jan <- na.omit(plasma.proteins.jan)
colnames(plasma.proteins.jan) <- c("X", "Subject", "Hour",
                                   "PC1", "PC2", "PC3", "PC4", "PC5")
PlotPCA(plasma.proteins.jan, 
        title = "Principal Components for Plasma Proteins(Visit 1)", 
        filename = "Plasma_Proteins_PCA_Jan.eps")

# Seth's updated Plasma Metabolites PCA excluding an outlier
plasma.metabolites.jan <- 
  read.csv(file.path("chronobiome", "Seth_PCA", "PC_Scores", "ScoresJan17",
                     "PlasmaMetabolites_BothVisits_ScoresFor5PCs.csv"))
colnames(plasma.metabolites.jan) <- c("Visit", "Subject", "Hour", "PC1",
                                      "PC2", "PC3", "PC4", "PC5")
PlotPCA(plasma.metabolites.jan,
        title = paste0("Principal Components for ",
                       "Plasma Metabolites(Visits", 
                       " 1 and 2)"), 
        filename = "Plasma_Metabs_PCA_Jan.eps")

# Macronutrients plot (Raw values rather than PCA)
Macronutrients <- read.csv(
  file.path("chronobiome", "MacroNutrientData_SmallerChunks_BothVisits.csv"))
postscript("Macronutrients.eps")
n1 <- qplot(x = Protein, y = TotalFat, data = Macronutrients, colour = Subject) +
  scale_colour_manual(values = colorblind_Palette)
n2 <- qplot(x = Protein, y = Carbohydrate, data = Macronutrients, colour = Subject) +
  scale_colour_manual(values=colorblind_Palette)
n3 <- qplot(x = Carbohydrate, y = TotalFat, data = Macronutrients, colour = Subject) +
  scale_colour_manual(values=colorblind_Palette)
grid.arrange(n1, n2, n3, ncol = 3, top = "Macronutrients(Visits 1 and 2)")
dev.off()