---
title: "gene-methylation-expression-mixomics"
output: html_document
editor_options: 
  chunk_output_type: console
---

# Set up R Markdown document

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = normalizePath("../output/gene-methylation-expression-mixomics/")) #Set root directory
```

```{r}
getwd()
```

#Install packages

```{r}
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
#  BiocManager::install("mixOmics") #Install BiocManager and mixOmics
# install.packages("plotly")
# install.packages("RColorBrewer")
require(devtools)
#install_github("mixOmicsTeam/mixOmics", ref = github_pull("197"), force = TRUE) #Install dev version: https://github.com/mixOmicsTeam/mixOmics/issues/196 (used for running SPLS)
install_github("mixOmicsTeam/mixOmics", ref = github_pull("234"), force = TRUE) #Install dev version: https://github.com/mixOmicsTeam/mixOmics/issues/233 (used for network plot customization)
require(mixOmics)
require(tidyverse)
require(plotly)
require(RColorBrewer)
```

```{r}
sessionInfo()
```

# Import and reformat data

```{r}
metadata <- read.csv("data/adult-meta.csv", header = TRUE) #Import metadata
metadata <- metadata %>%
  dplyr::select(., Sample.ID, Treatment, Sex, TreatmentN) %>%
  rename(., sample_id = Sample.ID) %>% 
  mutate(SexTreatment = case_when(TreatmentN == 1 ~ "Control Male",
                                  TreatmentN == 2 ~ "Control Female",
                                  TreatmentN == 3 ~ "Exposed Male",
                                  TreatmentN == 4 ~ "Exposed Female")) #Retain necessary columns, rename sample_id, and create a new descriptive treatment column
head(metadata) #Confirm formatting
```

## Gene expression data

```{r}
geneExpression <- read.csv("data/gene_fpkm.csv", header = TRUE, row.names = metadata$sample_id)
geneExpressionMod <- geneExpression %>%
  dplyr::select(., -c(1:6)) #Remove extra information to have just FPKM for each gene
head(geneExpressionMod) #Confirm formatting
ncol(geneExpressionMod) #Count number of genes with data
```

```{r}
geneExpressionHist <- geneExpressionMod %>%
  t(.) %>%
  as.data.frame(.) %>%
  drop_na(.) %>%
  rowwise() %>%
  summarise(range = max(c_across(where(is.numeric))) - min(c_across(where(is.numeric)))) #Calculate the range of gene expression values across each sample

plot_ly(x = geneExpressionHist$range, type = "histogram") #Interactive plot of gene expression value ranges
```

```{r}
geneExpressionHist <- geneExpressionMod %>%
  t(.) %>%
  as.data.frame(.) %>%
  drop_na(.) %>%
  mutate(var = matrixStats::rowVars(as.matrix(.))) #Calculate the variance of gene expression values across each sample

range(geneExpressionHist$var)

plot_ly(x = geneExpressionHist$var, type = "histogram") #Interactive plot of gene expression value ranges
```

Based on the interactive histogram, I'll remove data where the range is < 10.

```{r}
#Transpose expression data
#Remove rows with missing data
#Rowwise range operation
#Remove rows with near zero variance
geneExpressionFilt <- geneExpressionMod %>%
  t(.) %>%
  as.data.frame(.) %>%
  drop_na(.) %>%
  mutate(var = matrixStats::rowVars(as.matrix(.))) %>%
  filter(., var > 5) %>%
  select(., !var) %>%
  t(.) %>%
  as.data.frame(.)
head(geneExpressionFilt) #Confirm formatting
sum(is.na(geneExpressionFilt)) #Confirm there are no NAs
ncol(geneExpressionFilt) #Number of genes remaining
```

```{r}
exp.pca <- pca(geneExpressionFilt) #PCA analysis of expression
plotIndiv(exp.pca, ind.names = FALSE,
          group = metadata$Treatment,
          pch = as.factor(metadata$Sex), col = c("#74C476", "#006D2C"),
          cex = 10, size.axis = 8, size.legend = 8,
          legend = TRUE, title = 'Gene Expression',
          legend.title = 'Treatment', legend.title.pch = 'Sex')
ggsave("output/gene-methylation-expression-mixomics/all-exp-pca.pdf", height = 8.5, width = 11)
```


## Gene methylation data

```{r}
geneMethylation <- read.csv("output/40-gene-methylaiton.csv", header = TRUE, row.names = 2)
geneMethylationMod <- geneMethylation %>%
  select(., -1) #Remove extra column
head(geneMethylationMod) #Confirm formatting
ncol(geneMethylationMod) #Count number of genes with data
```

```{r}
#Transpose methylation data
#Remove rows with missing data
#Rowwise range operation
#Remove rows without big differences in gene methylation
geneMethylationFilt <- geneMethylationMod %>%
  t(.) %>%
  as.data.frame(.) %>%
  drop_na(.) %>%
  mutate(var = matrixStats::rowVars(as.matrix(.))) %>%
  filter(., var > 5) %>%
  select(., !var) %>%
  t(.) %>%
  as.data.frame(.)
head(geneMethylationFilt) #Confirm formatting
sum(is.na(geneMethylationFilt)) #Confirm there are no NAs
ncol(geneMethylationFilt) #Number of genes remaining
```


```{r}
meth.pca <- pca(geneMethylationFilt) #PCA analysis of methylation
plotIndiv(meth.pca, ind.names = FALSE,
          group = metadata$Treatment,
          pch = as.factor(metadata$Sex), col = c("#74C476", "#006D2C"),
          cex = 10, size.axis = 8, size.legend = 8,
          legend = TRUE, title = 'Gene Methylation',
          legend.title = 'Treatment', legend.title.pch = 'Sex')
ggsave("output/gene-methylation-expression-mixomics/all-meth-pca.pdf", height = 8.5, width = 11)
```


# Sparse PLS

Followed [`mixOmics` manual](https://mixomicsteam.github.io/Bookdown/pls.html#inputs-and-outputs-1) for analysis methods.

## Run SPLS with all samples

I'm going to run an SPLS with default settings for all parameters to confirm there are large sex differences that should be accounted for in follow-up analyses. Keep top 100 loadings from each dataset. 

```{r}
geneResult.spls <- spls(X = geneMethylationFilt, Y = geneExpressionFilt, 
                        keepX = c(100, 100), keepY = c(100,100)) #Run SPLS
```

```{r}
#pdf("../output/gene-methylation-expression-mixomics/yrv-geneSPLS-Full1.pdf", width = 11, height = 8.5)
plotIndiv(geneResult.spls, group = metadata$Sex,
          legend = TRUE,
          legend.title = 'Sex',
          ind.names = metadata$sample_id,
          title = "gene: sPLS") #Plot samples. Group by sex and add titles for treatment
#dev.off()
```

```{r}
#pdf("../output/gene-methylation-expression-mixomics/yrv-geneSPLS-Full2.pdf", width = 11, height = 8.5)
plotIndiv(geneResult.spls, group= metadata$Sex,
          pch = metadata$TreatmentN,
          rep.space = "XY-variate",  legend = TRUE,
          legend.title = "Sex", legend.title.pch = "Sex + Treatment",
          ind.names = FALSE,
          title = "gene: sPLS") #Plot samples. Group by sex and add symbols for sex + treatment
#dev.off()
```

```{r}
geneArrow <- plotArrow(geneResult.spls, group = metadata$SexTreatment, ind.names = FALSE,
                       legend = TRUE, legend.title = "Treatment + Sex") + 
  labs(x = "PLS Comp 1", y = "PLS Comp 2") #Group by treatment and sex
geneArrow$data
plot(geneArrow)
```


## Separate sex-specific data

```{r}
metadataFem <- metadata %>%
  filter(., Sex == "F") #Select female metadata
metadataMale <- metadata %>%
  filter(., Sex == "M") #Select male metadata
```

### Gene expression

```{r}
geneExpressionFem <- geneExpressionMod %>%
  filter(., grepl("F", rownames(.), fixed = TRUE)) %>%
  t(.) %>%
  as.data.frame(.) %>%
  drop_na(.) %>%
  mutate(var = matrixStats::rowVars(as.matrix(.))) %>%
  filter(., var > 5) %>%
  dplyr::select(., !var) %>%
  t(.) %>%
  as.data.frame(.) #Filter female samples only, remove genes with missing expression data (NA) for any sample, remove samples with low variation, and make sample IDs rownames
head(geneExpressionFem) #Confirm formatting
sum(is.na(geneExpressionFem)) #Confirm there are no NAs
ncol(geneExpressionFem) #Count genes with data for all samples
```

```{r}
geneExpressionMale <- geneExpressionMod %>%
  filter(., grepl("M", rownames(.), fixed = TRUE)) %>%
  t(.) %>%
  as.data.frame(.) %>%
  drop_na(.) %>%
  mutate(var = matrixStats::rowVars(as.matrix(.))) %>%
  filter(., var > 5) %>%
  dplyr::select(., !var) %>%
  t(.) %>%
  as.data.frame(.) #Filter male samples only, remove genes with missing expression data (NA) for any sample, remove samples with low variation, and make sample IDs rownames
head(geneExpressionMale) #Confirm formatting
sum(is.na(geneExpressionMale)) #Confirm there are no NAs
ncol(geneExpressionMale) #Count genes with data for all samples
```

### Gene methylation

```{r}
geneMethylationFem <- geneMethylationMod %>%
  filter(., grepl("F", rownames(.), fixed = TRUE)) %>% 
  t(.) %>%
  as.data.frame(.) %>%
  drop_na(.) %>%
  mutate(var = matrixStats::rowVars(as.matrix(.))) %>%
  filter(., var > 5) %>%
  dplyr::select(., !var) %>%
  t(.) %>%
  as.data.frame(.) #Filter female samples only
head(geneMethylationFem) #Confirm formatting
sum(is.na(geneMethylationFem)) #Confirm there are no NAs
ncol(geneMethylationFem) #Count genes with data for all samples
```

```{r}
geneMethylationMale <- geneMethylationMod %>%
  filter(., grepl("M", rownames(.), fixed = TRUE)) %>% 
  t(.) %>%
  as.data.frame(.) %>%
  drop_na(.) %>%
  mutate(var = matrixStats::rowVars(as.matrix(.))) %>%
  filter(., var > 5) %>%
  dplyr::select(., !var) %>%
  t(.) %>%
  as.data.frame(.) #Filter male samples only
head(geneMethylationMale) #Confirm formatting
sum(is.na(geneMethylationMale)) #Confirm there are no NAs
ncol(geneMethylationMale) #Count genes with data for all samples
```

## Female SPLS

### Tune parameters

#### Number of components

```{r}
geneResult.plsFem <- pls(geneMethylationFem, geneExpressionFem, ncomp = 4) #Run a PLS with a sufficient number of components
```

AH had an error at this point as of 15 Aug 2023: Error in Ypred[omit, , h] <- Y.hat[, , 1] : 
  number of items to replace is not a multiple of replacement length. This error continues in code below as well. 
  
Tried changing validation to "loo" instead of "Mfold". Didn't work. 

AH left off here. 





  
```{r}
genePerf.plsFem <- perf(geneResult.plsFem, validation = "Mfold", folds = 2,
                        progressBar = FALSE, nrepeat = 10) #Use perf for repeated k-fold cross-validation to determine the number of components that should be included
```


```{r}
#pdf(“../output/gene-methylation-expression-mixomics/yrv-geneSPLS-Fem-CompTune.pdf”, width = 11, height = 8.5)
plot(genePerf.plsFem$Q2.total) #Plot perf results
abline(h = 0.0975) #Only include components with value is ≤ 0.0975
#dev.off()
```

#### `keepX`: Number of variables per component

```{r}
list.keepX <- c(2:10, 15, 20, 25, 30) #Number of drivers to pull
tune.spls.RSS.FemX <- tune.spls(geneMethylationFem, geneExpressionFem, ncomp = 4,
                                test.keepX = list.keepX,
                                validation = "Mfold", folds = 2,
                                nrepeat = 10, progressBar = TRUE,
                                measure = "RSS") #Tune parameters based on Mean Absolute Value with 5 folds and 1000 repeats
```

```{r}
#pdf(“../output/gene-methylation-expression-mixomics/yrv-geneSPLS-Fem-XTune.pdf”, width = 11, height = 8.5)
plot(tune.spls.MAE.FemX, legend.position = "topright") #Plot tuning results
#dev.off()
```

#### `keepY`: Number of variables per component

```{r}
list.keepY <- c(2:10, 15, 20, 25, 30) #Number of drivers to pull
tune.spls.MAE.FemY <- tune.spls(geneMethylationFem, geneExpressionFem, ncomp = FIX,
                                test.keepY = list.keepY,
                                validation = "Mfold", folds = 10,
                                nrepeat = 1000, progressBar = FALSE,
                                measure = "MAE") #Tune parameters based on Mean Absolute Value with 5 folds and 1000 repeats
```

```{r}
#pdf(“../output/gene-methylation-expression-mixomics/yrv-geneSPLS-Fem-YTune.pdf”, width = 11, height = 8.5)
plot(tune.spls.MAE.FemY, legend.position = "topright") #Plot tuning results
#dev.off()
```

### Run SPLS 

```{r}
geneResult.splsFem <- spls(X = geneMethylationFem, Y = geneExpressionFem,
                           keepX = c(25, 25), keepY = c(25,25)) #Run SPLS for female samples
```

#### Visualize results

```{r}
plotIndiv(geneResult.splsFem, group = metadataFem$Treatment,
          legend = TRUE,
          legend.title = 'Treatment',
          ind.names = metadataFem$sample_id,
          title = "gene: sPLS") #Plot samples. Group by treatment
ggsave("../output/gene-methylation-expression-mixomics/geneSPLS-fem.pdf", width = 11, height = 8.5)
```

```{r}
geneArrowFem <- plotArrow(geneResult.splsFem, group = metadataFem$Treatment, ind.names = FALSE,
                          legend = TRUE, legend.title = "Treatment") + 
  labs(x = "PLS Comp 1", y = "PLS Comp 2") #Group by treatment. Short arrows indicate good agreement found by the PLS between both data sets
geneArrowFem$data
plot(geneArrowFem)
ggsave("../output/gene-methylation-expression-mixomics/arrowPlot-fem.pdf", width = 11, height = 8.5)
```

```{r}
plotVar(geneResult.splsFem, cex = c(3,2), var.names = FALSE, cutoff = 0.8,
        title = "Female Correlation Circle Plot", legend = FALSE,
        X.label = "PLS Comp 1", Y.label = "PLS Comp 2") #Create correlation circle plot using 0.8 correlation cutoff. Variables that are close together are highly correlated. Variables with < 90º angle between them have a positive correlation, with variables having > 90º angle between them having a negative correlation. Variables at a right angle have 0 correlation.
ggsave("../output/gene-methylation-expression-mixomics/corrCircle-fem.pdf", width = 11, height = 8.5)
```
```{r}
femCircCoordinates <- plotVar(geneResult.splsFem, plot = FALSE) 
head(femCircCoordinates)
```

```{r}
cim(margins = c(10,10), geneResult.splsFem,
    xlab = "Expression", ylab = "Methylation",
    save = "pdf", name.save = "../output/gene-methylation-expression-mixomics/PLScim-fem") #Show Pearson correlations between X and Y variables on the first component. Cell color is based on values of the similarity matrix.
```

```{r}
femNetwork <- network(geneResult.splsFem, cutoff = 0.8,
                      row.names = FALSE, col.names = FALSE, 
                      color.node = c("grey50", "grey50"), shape.node = c("circle", "rectangle"), size.node = 0.2,
                      color.edge = rev(brewer.pal(10, "RdBu")),
                      save = "pdf", name.save = "../output/gene-methylation-expression-mixomics/network-fem") #Create and save network plot. Methylation = blue circles, expression = orange rectangles. Blue lines = neg corr, red lines = pos corr. Only use correlations > 0.8 cutoff. The network plot is split based on components
```

```{r}
write.graph(femNetwork$gR, "../output/gene-methylation-expression-mixomics/femNetwork.gml", format = "gml") #Export for cytoscape
```

### Component-level information

#### Component 1

```{r}
femNetworkComp1 <- network(geneResult.splsFem, cutoff = 0.8, comp = 1,
                           row.names = FALSE, col.names = FALSE, 
                           color.node = c("lightblue", "orange"), shape.node = c("circle", "rectangle"), size.node = 0.2,
                           color.edge = rev(brewer.pal(10, "RdBu")),
                           save = "pdf", name.save = "../output/gene-methylation-expression-mixomics/comp1Network-fem") #Create and save network plot for component 1. Methylation = blue circles, expression = orange rectangles. Blue lines = neg corr, red lines = pos corr. Only use correlations > 0.8 cutoff.
```

```{r}
femSelectedVariables1 <- selectVar(geneResult.splsFem, comp = 1) #Selected variables for first component
femSelectedVariables1$X$name #Selected genes from methylation data
femSelectedVariables1$Y$name #Selected genes from expression data
```

```{r}
range(femSelectedVariables1$X$value) #Selected genes from methylation data: -0.42825345  0.04797473
range(femSelectedVariables1$Y$value) #Selected genes from expression data: -0.3044497  0.4415122
```

```{r}
#pdf("comp1Loadings-fem.pdf", width = 11, height = 8.5)
plotLoadings(geneResult.splsFem, comp = 1, 
             size.name = rel(0.7),
             xlim = matrix(c(-0.45, 0.1, -0.35, 0.45), nrow = 2, byrow = TRUE),
             subtitle = c("Methylation", "Expression"), size.subtitle = 2) #Loadings for first component
#dev.off()
```

```{r}
rbind(as.data.frame(femSelectedVariables1$X$value),
      as.data.frame(femSelectedVariables1$Y$value)) %>% 
  mutate(., block = c(rep("X", times = 25),
                      rep("Y", times = 25))) %>%
  rownames_to_column(., var = "gene") %>%
  write.csv(., "../output/gene-methylation-expression-mixomics/comp1SelectedVariables-Fem.csv", quote = FALSE, row.names = FALSE)
```

##### Annotate genes

```{r}
mRNATrack <- read.delim("../genome-features/C_virginica-3.0-mRNA.gff", header = FALSE) %>%
  separate(., V9, into = c("V10", "V11"), sep = "ID=rna-") %>%
  separate(., V11, into = c("V12", "V13"), sep = ";Parent=") %>%
  separate(., V13, into = c("V14", "V15"), sep = ";Dbxref=") %>%
  separate(., V15, into = c("V16", "V17"), sep = ";product=") %>%
  separate(., V17, into = c("V18", "V19"), sep = ";transcript_id=") %>%
  select(., c("V12", "V14", "V18")) %>%
  rename(., transcript = V12, gene = V14, product = V18) #Winnow down mRNA track
mRNATrack$gene <- gsub("gene-", "gene.", mRNATrack$gene) #Make gene naming consistent
head(mRNATrack) #Confirm changes
```

```{bash}
wget https://gannet.fish.washington.edu/spartina/paper-gonad-meth/analyses/2018-12-02-Gene-Enrichment-Analysis/_blast-GO-unfolded.sorted > ../genome-features/_blast-GO-unfolded.sorted
```

```{bash}
wget https://gannet.fish.washington.edu/spartina/project-gigas-oa-meth/output/GO-term-annotation/GO-GOslim.sorted > ../genome-features/GO-GOslim.sorted
```

```{r}
GOterms <- read.delim("../genome-features/_blast-GO-unfolded.sorted", header = FALSE, col.names = c("GO", "transcript"))
head(GOterms)
```

```{r}
GOslim <- read.delim("../genome-features/GO-GOslim.sorted", header = FALSE, col.names = c("GO", "GOterm", "GOslim", "process"))
head(GOslim)
```

```{r}
annotfemSelectedVariables1X <- femSelectedVariables1$X %>%
  as.data.frame(.) %>%
  rename(., gene = name) %>%
  left_join(., mRNATrack, by = "gene") %>% 
  left_join(., GOterms, by = "transcript") %>% 
  left_join(., GOslim, by = "GO") #Annotate genes
head(annotfemSelectedVariables1X) #Confirm annotation
```


```{r}
annotfemSelectedVariables1X %>%
  select(., c(gene, value.var, product)) %>%
  separate(., product, into = c("product", "transcriptVar"), sep = "%2C ") %>%
  select(., -transcriptVar) %>%
  distinct(.)
```


```{r}
annotfemSelectedVariables1XP <- annotfemSelectedVariables1X %>%
  dplyr::filter(., process == "P")
as.data.frame(table(annotfemSelectedVariables1XP$GOslim))
```

```{r}
annotfemSelectedVariables1XF <- annotfemSelectedVariables1X %>%
  dplyr::filter(., process == "F")
as.data.frame(table(annotfemSelectedVariables1XF$GOslim))
```

```{r}
annotfemSelectedVariables1Y <- femSelectedVariables1$Y %>%
  as.data.frame(.) %>%
  rename(., gene = name) %>%
  left_join(., mRNATrack, by = "gene") %>% 
  left_join(., GOterms, by = "transcript") %>% 
  left_join(., GOslim, by = "GO") #Annotate genes
head(annotfemSelectedVariables1Y) #Confirm annotation
```

```{r}
annotfemSelectedVariables1Y %>%
  select(., c(gene, value.var, product)) %>%
  separate(., product, into = c("product", "transcriptVar"), sep = "%2C ") %>%
  select(., -transcriptVar) %>%
  distinct(.)
```

```{r}
annotfemSelectedVariables1YP <- annotfemSelectedVariables1Y %>%
  dplyr::filter(., process == "P")
as.data.frame(table(annotfemSelectedVariables1YP$GOslim))
```

```{r}
annotfemSelectedVariables1YF <- annotfemSelectedVariables1Y %>%
  dplyr::filter(., process == "F")
as.data.frame(table(annotfemSelectedVariables1YF$GOslim))
```

#### Component 2

```{r}
femNetworkComp2 <- network(geneResult.splsFem, cutoff = 0.8, comp = 2,
                           row.names = FALSE, col.names = FALSE, 
                           color.node = c("lightblue", "orange"), shape.node = c("circle", "rectangle"), size.node = 0.2,
                           color.edge = rev(brewer.pal(10, "RdBu")),
                           save = "pdf", name.save = "../output/gene-methylation-expression-mixomics/comp2Network-fem") #Create and save network plot for component 1. Methylation = blue circles, expression = orange rectangles. Blue lines = neg corr, red lines = pos corr. Only use correlations > 0.8 cutoff.
```

```{r}
femSelectedVariables2 <- selectVar(geneResult.splsFem, comp = 2) #Selected variables for second component
femSelectedVariables2$X$name #Selected genes from methylation data
femSelectedVariables2$Y$name #Selected genes from expression data
```

```{r}
range(femSelectedVariables2$X$value) #Selected genes from methylation data
range(femSelectedVariables2$Y$value) #Selected genes from expression data
```

```{r}
rbind(as.data.frame(femSelectedVariables2$X$value),
      as.data.frame(femSelectedVariables2$Y$value)) %>% 
  mutate(., block = c(rep("X", times = 25),
                      rep("Y", times = 25))) %>%
  rownames_to_column(., var = "gene") %>% write.csv(., "../output/gene-methylation-expression-mixomics/comp2SelectedVariables-Fem.csv", quote = FALSE, row.names = FALSE)
```

```{r}
#pdf("comp2Loadings-fem.pdf", width = 11, height = 8.5)
plotLoadings(geneResult.splsFem, comp = 2, 
             size.name = rel(0.7),
             xlim = matrix(c(0, 0.6, 0, 0.4), nrow = 2, byrow = TRUE),
             subtitle = c("Methylation", "Expression"), size.subtitle = 2) #Loadings for first component
#dev.off()
```

##### Annotate genes

```{r}
annotfemSelectedVariables2X <- femSelectedVariables2$X %>%
  as.data.frame(.) %>%
  rename(., gene = name) %>%
  left_join(., mRNATrack, by = "gene") %>% 
  left_join(., GOterms, by = "transcript") %>% 
  left_join(., GOslim, by = "GO") #Annotate genes
head(annotfemSelectedVariables2X) #Confirm annotation
```

```{r}
annotfemSelectedVariables2X %>%
  select(., c(gene, value.var, product)) %>%
  separate(., product, into = c("product", "transcriptVar"), sep = "%2C ") %>%
  select(., -transcriptVar) %>%
  distinct(.)
```

```{r}
annotfemSelectedVariables2XP <- annotfemSelectedVariables2X %>%
  dplyr::filter(., process == "P")
as.data.frame(table(annotfemSelectedVariables2XP$GOslim))
```

```{r}
annotfemSelectedVariables2XF <- annotfemSelectedVariables2X %>%
  dplyr::filter(., process == "F")
as.data.frame(table(annotfemSelectedVariables2XF$GOslim))
```

```{r}
annotfemSelectedVariables2Y <- femSelectedVariables2$Y %>%
  as.data.frame(.) %>%
  rename(., gene = name) %>%
  left_join(., mRNATrack, by = "gene") %>% 
  left_join(., GOterms, by = "transcript") %>% 
  left_join(., GOslim, by = "GO") #Annotate genes
head(annotfemSelectedVariables2Y) #Confirm annotation
```

```{r}
annotfemSelectedVariables2Y %>%
  select(., c(gene, value.var, product)) %>%
  separate(., product, into = c("product", "transcriptVar"), sep = "%2C ") %>%
  select(., -transcriptVar) %>%
  distinct(.)
```

```{r}
annotfemSelectedVariables2YP <- annotfemSelectedVariables2Y %>%
  dplyr::filter(., process == "P")
as.data.frame(table(annotfemSelectedVariables2YP$GOslim))
```

```{r}
annotfemSelectedVariables2YF <- annotfemSelectedVariables2Y %>%
  dplyr::filter(., process == "F")
as.data.frame(table(annotfemSelectedVariables2YF$GOslim))
```
#### Create table with loading values and annotations

```{r}
rbind(as.data.frame(femSelectedVariables1$X$value),
      as.data.frame(femSelectedVariables1$Y$value),
      as.data.frame(femSelectedVariables2$X$value),
      as.data.frame(femSelectedVariables2$Y$value)) %>% 
  mutate(., block = c(rep("X", times = 25),
                      rep("Y", times = 25),
                      rep("X", times = 25),
                      rep("Y", times = 25))) %>%
  mutate(., component = c(rep("1", times = 50),
                          rep("2", times = 50))) %>%
  rownames_to_column(., var = "gene") %>%
  left_join(., mRNATrack, by = "gene") %>% 
  left_join(., GOterms, by = "transcript") %>% 
  left_join(., GOslim, by = "GO") %>%
  arrange(., component, block, value.var, process, transcript, GO) %>%
  distinct(.) %>%
  dplyr::group_by(gene, value.var, block) %>%
  mutate(transcript = paste(transcript, collapse = ";")) %>%
  mutate(product = paste(product, collapse = ";")) %>%
  mutate(GO = paste(GO, collapse = ";")) %>%
  mutate(GOterm = paste(GOterm, collapse = ";")) %>%
  mutate(GOslim = paste(GOslim, collapse = ";")) %>%
  mutate(process = paste(process, collapse = ";")) %>%
  ungroup(.) %>%
  distinct(.) %>%
  mutate(transcript = lapply(strsplit(transcript, split = ";"), unique)) %>%
  mutate(transcript = sapply(transcript, toString)) %>%
  mutate(transcript = gsub(transcript, pattern = ", ", replacement = ";")) %>%
  mutate(transcript = na_if(x = transcript, y = "NA")) %>%
  mutate(product = lapply(strsplit(product, split = ";"), unique)) %>%
  mutate(product = sapply(product, toString)) %>%
  mutate(product = gsub(product, pattern = ", ", replacement = ";")) %>%
  mutate(product = na_if(x = product, y = "NA")) %>%
  mutate(GO = lapply(strsplit(GO, split = ";"), unique)) %>%
  mutate(GO = sapply(GO, toString)) %>%
  mutate(GO = gsub(GO, pattern = ", ", replacement = ";")) %>%
  mutate(GO = na_if(x = GO, y = "NA")) %>%
  mutate(GOterm = lapply(strsplit(GOterm, split = ";"), unique)) %>%
  mutate(GOterm = sapply(GOterm, toString)) %>%
  mutate(GOterm = gsub(GOterm, pattern = ", ", replacement = ";")) %>%
  mutate(GOterm = na_if(x = GOterm, y = "NA")) %>%
  mutate(GOslim = lapply(strsplit(GOslim, split = ";"), unique)) %>%
  mutate(GOslim = sapply(GOslim, toString)) %>%
  mutate(GOslim = gsub(GOslim, pattern = ", ", replacement = ";")) %>%
  mutate(GOslim = na_if(x = GOslim, y = "NA")) %>%
  mutate(process = lapply(strsplit(process, split = ";"), unique)) %>%
  mutate(process = sapply(process, toString)) %>%
  mutate(process = gsub(process, pattern = ", ", replacement = ";")) %>%
  mutate(process = na_if(x = process, y = "NA")) %>%
  distinct(.) %>%
  write.table("selectedVariables-Annot-Fem.txt", quote = FALSE, sep = "\t", row.names = FALSE) #Combine driver gene information for both X and Y datasets for the first component. Add block (X or Y) and component (1 or 2) information for gene methylation or expression datasets. Join with mRNA, GO term, and GO Slim information for annotations. Arrange by component, block, loading value, process, transcript, and GO, then retain only distinct entries. Group by gene, loading values, and block to collapse rows across all descriptive columns such that there is only one row per driver gene. Ungroup and retain only distinct rows. Use a series of apply commands to split strings in each column by ; then remove any duplicate values. Replace commas with semicolon delimiters and make sure NAs are identified. Retain distinct rows. Save as a txt file
```

```{bash}
head selectedVariables-Annot-Fem.txt
```

## Male SPLS

### Tune parameters

#### Number of components

```{r}
geneResult.plsMale <- pls(geneMethylationMale, geneExpressionMale, ncomp = 4) #Run a PLS with a sufficient number of components
```

```{r}
genePerf.plsMale <- perf(geneResult.plsMale, validation = "Mfold", folds = 2,
                         progressBar = FALSE, nrepeat = 10) #Use perf for repeated k-fold cross-validation to determine the number of components that should be included
```

```{r}
#pdf(“../output/gene-methylation-expression-mixomics/yrv-geneSPLS-Male-CompTune.pdf”, width = 11, height = 8.5)
plot(genePerf.plsMale$Q2.total) #Plot perf results
abline(h = 0.0975) #Only include components with value is ≤ 0.0975
#dev.off()
```

#### `keepX`: Number of variables per component

```{r}
list.keepX <- c(2:10, 15, 20, 25, 30) #Number of drivers to pull
tune.spls.RSS.MaleX <- tune.spls(geneMethylationMale, geneExpressionMale, ncomp = FIX,
                                 test.keepX = list.keepX,
                                 validation = "Mfold", folds = 2,
                                 nrepeat = 10, progressBar = TRUE,
                                 measure = "RSS") #Tune parameters based on Mean Absolute Value with 5 folds and 1000 repeats
```

```{r}
#pdf(“../output/gene-methylation-expression-mixomics/yrv-geneSPLS-Male-XTune.pdf”, width = 11, height = 8.5)
plot(tune.spls.MAE.MaleX, legend.position = "topright") #Plot tuning results
#dev.off()
```

#### `keepY`: Number of variables per component

```{r}
list.keepY <- c(2:10, 15, 20, 25, 30) #Number of drivers to pull
tune.spls.MAE.MaleY <- tune.spls(geneMethylationMale, geneExpressionMale, ncomp = FIX,
                                 test.keepY = list.keepY,
                                 validation = "Mfold", folds = 10,
                                 nrepeat = 1000, progressBar = FALSE,
                                 measure = "MAE") #Tune parameters based on Mean Absolute Value with 5 folds and 1000 repeats
```

```{r}
#pdf(“../output/gene-methylation-expression-mixomics/yrv-geneSPLS-Male-YTune.pdf”, width = 11, height = 8.5)
plot(tune.spls.MAE.MaleY, legend.position = "topright") #Plot tuning results
#dev.off()
```

### Run SPLS 

```{r}
geneResult.splsMale <- spls(X = geneMethylationMale, Y = geneExpressionMale,
                            keepX = c(25, 25), keepY = c(25,25)) #Run SPLS for Maleale samples
```

#### Visualize results

```{r}
plotIndiv(geneResult.splsMale, group = metadataMale$Treatment,
          legend = TRUE,
          legend.title = 'Treatment',
          ind.names = metadataMale$sample_id,
          title = "gene: sPLS") #Plot samples. Group by treatment
ggsave("../output/gene-methylation-expression-mixomics/geneSPLS-Male.pdf", width = 11, height = 8.5)
```

```{r}
geneArrowMale <- plotArrow(geneResult.splsMale, group = metadataMale$Treatment, ind.names = FALSE,
                           legend = TRUE, legend.title = "Treatment") + 
  labs(x = "PLS Comp 1", y = "PLS Comp 2") #Group by treatment. Short arrows indicate good agreement found by the PLS between both data sets
geneArrowMale$data
plot(geneArrowMale)
ggsave("../output/gene-methylation-expression-mixomics/arrowPlot-Male.pdf", width = 11, height = 8.5)
```

```{r}
plotVar(geneResult.splsMale, cex = c(3,2), var.names = FALSE, cutoff = 0.8,
        title = "Male Correlation Circle Plot", legend = FALSE,
        X.label = "PLS Comp 1", Y.label = "PLS Comp 2") #Create correlation circle plot using 0.8 correlation cutoff. Variables that are close together are highly correlated. Variables with < 90º angle between them have a positive correlation, with variables having > 90º angle between them having a negative correlation. Variables at a right angle have 0 correlation.
ggsave("../output/gene-methylation-expression-mixomics/corrCircle-Male.pdf", width = 11, height = 8.5)
```

```{r}
MaleCircCoordinates <- plotVar(geneResult.splsMale, plot = FALSE) 
head(MaleCircCoordinates)
```

```{r}
cim(margins = c(10,10), geneResult.splsMale,
    xlab = "Expression", ylab = "Methylation",
    save = "pdf", name.save = "../output/gene-methylation-expression-mixomics/PLScim-Male") #Show Pearson correlations between X and Y variables on the first component. Cell color is based on values of the similarity matrix.
```

```{r}
MaleNetwork <- network(geneResult.splsMale, cutoff = 0.8,
                       row.names = FALSE, col.names = FALSE, 
                      color.node = c("grey50", "grey50"), shape.node = c("circle", "rectangle"), size.node = 0.2,
                      color.edge = rev(brewer.pal(10, "RdBu")),
                       save = "pdf", name.save = "../output/gene-methylation-expression-mixomics/network-Male") #Create and save network plot. Methylation = blue circles, expression = orange rectangles. Blue lines = neg corr, red lines = pos corr. Only use correlations > 0.8 cutoff. The network plot is split based on components
```

```{r}
write.graph(MaleNetwork$gR, "../output/gene-methylation-expression-mixomics/MaleNetwork.gml", format = "gml") #Export for cytoscape
```

### Component-level information

#### Component 1

```{r}
MaleNetworkComp1 <- network(geneResult.splsMale, cutoff = 0.8, comp = 1,
                            row.names = FALSE, col.names = FALSE, 
                            color.node = c("lightblue", "orange"), shape.node = c("circle", "rectangle"), size.node = 0.2,
                            color.edge = rev(brewer.pal(10, "RdBu")),
                            save = "pdf", name.save = "../output/gene-methylation-expression-mixomics/comp1Network-Male") #Create and save network plot for component 1. Methylation = blue circles, expression = orange rectangles. Blue lines = neg corr, red lines = pos corr. Only use correlations > 0.8 cutoff.
```

```{r}
MaleSelectedVariables1 <- selectVar(geneResult.splsMale, comp = 1) #Selected variables for first component
MaleSelectedVariables1$X$name #Selected genes from methylation data
MaleSelectedVariables1$Y$name #Selected genes from expression data
```

```{r}
range(MaleSelectedVariables1$X$value) #Selected genes from methylation data
range(MaleSelectedVariables1$Y$value) #Selected genes from expression data
```

```{r}
rbind(as.data.frame(MaleSelectedVariables1$X$value),
      as.data.frame(MaleSelectedVariables1$Y$value)) %>% 
  mutate(., block = c(rep("X", times = 25),
                      rep("Y", times = 25))) %>%
  rownames_to_column(., var = "gene") %>% write.csv(., "../output/gene-methylation-expression-mixomics/comp1SelectedVariables-Male.csv", quote = FALSE, row.names = FALSE)
```

```{r}
#pdf("comp1Loadings-male.pdf", width = 11, height = 8.5)
plotLoadings(geneResult.splsMale, comp = 1, 
             size.name = rel(0.7),
             xlim = matrix(c(-0.2, 0.35, 0, 0.45), nrow = 2, byrow = TRUE),
             subtitle = c("Methylation", "Expression"), size.subtitle = 2) #Loadings for first component
#dev.off()
```

##### Annotate genes

```{r}
annotmaleSelectedVariables1X <- MaleSelectedVariables1$X %>%
  as.data.frame(.) %>%
  rename(., gene = name) %>%
  left_join(., mRNATrack, by = "gene") %>% 
  left_join(., GOterms, by = "transcript") %>% 
  left_join(., GOslim, by = "GO") #Annotate genes
head(annotmaleSelectedVariables1X) #Confirm annotation
```

```{r}
annotmaleSelectedVariables1X %>%
  select(., c(gene, value.var, product)) %>%
  separate(., product, into = c("product", "transcriptVar"), sep = "%2C ") %>%
  select(., -transcriptVar) %>%
  distinct(.)
```

```{r}
annotmaleSelectedVariables1XP <- annotmaleSelectedVariables1X %>%
  dplyr::filter(., process == "P")
as.data.frame(table(annotmaleSelectedVariables1XP$GOslim))
```

```{r}
annotmaleSelectedVariables1XF <- annotmaleSelectedVariables1X %>%
  dplyr::filter(., process == "F")
as.data.frame(table(annotmaleSelectedVariables1XF$GOslim))
```

```{r}
annotmaleSelectedVariables1Y <- MaleSelectedVariables1$Y %>%
  as.data.frame(.) %>%
  rename(., gene = name) %>%
  left_join(., mRNATrack, by = "gene") %>% 
  left_join(., GOterms, by = "transcript") %>% 
  left_join(., GOslim, by = "GO") #Annotate genes
head(annotmaleSelectedVariables1Y) #Confirm annotation
```

```{r}
annotmaleSelectedVariables1Y %>%
  select(., c(gene, value.var, product)) %>%
  separate(., product, into = c("product", "transcriptVar"), sep = "%2C ") %>%
  select(., -transcriptVar) %>%
  distinct(.)
```

```{r}
annotmaleSelectedVariables1YP <- annotmaleSelectedVariables1Y %>%
  dplyr::filter(., process == "P")
as.data.frame(table(annotmaleSelectedVariables1YP$GOslim))
```

```{r}
annotmaleSelectedVariables1YF <- annotmaleSelectedVariables1Y %>%
  dplyr::filter(., process == "F")
as.data.frame(table(annotmaleSelectedVariables1YF$GOslim))
```

#### Component 2

```{r}
MaleNetworkComp2 <- network(geneResult.splsMale, cutoff = 0.8, comp = 2,
                            row.names = FALSE, col.names = FALSE, 
                            color.node = c("lightblue", "orange"), shape.node = c("circle", "rectangle"), size.node = 0.2,
                            color.edge = rev(brewer.pal(10, "RdBu")),
                            save = "pdf", name.save = "../output/gene-methylation-expression-mixomics/comp2Network-Male") #Create and save network plot for component 1. Methylation = blue circles, expression = orange rectangles. Blue lines = neg corr, red lines = pos corr. Only use correlations > 0.8 cutoff.
```

```{r}
MaleSelectedVariables2 <- selectVar(geneResult.splsMale, comp = 2) #Selected variables for second component
MaleSelectedVariables2$X$name #Selected genes from methylation data
MaleSelectedVariables2$Y$name #Selected genes from expression data
```

```{r}
range(MaleSelectedVariables2$X$value) #Selected genes from methylation data
range(MaleSelectedVariables2$Y$value) #Selected genes from expression data
```

```{r}
rbind(as.data.frame(MaleSelectedVariables2$X$value),
      as.data.frame(MaleSelectedVariables2$Y$value)) %>% 
  mutate(., block = c(rep("X", times = 25),
                      rep("Y", times = 25))) %>%
  rownames_to_column(., var = "gene") %>%
  write.csv(., "../output/gene-methylation-expression-mixomics/comp2SelectedVariables-Male.csv", quote = FALSE, row.names = FALSE)
```

```{r}
#pdf("comp2Loadings-male.pdf", width = 11, height = 8.5)
plotLoadings(geneResult.splsMale, comp = 2, 
             size.name = rel(0.7),
             xlim = matrix(c(0, 0.45, 0, 0.4), nrow = 2, byrow = TRUE),
             subtitle = c("Methylation", "Expression"), size.subtitle = 2) #Loadings for first component
#dev.off()
```

##### Annotate genes

```{r}
annotmaleSelectedVariables2X <- MaleSelectedVariables2$X %>%
  as.data.frame(.) %>%
  rename(., gene = name) %>%
  left_join(., mRNATrack, by = "gene") %>% 
  left_join(., GOterms, by = "transcript") %>% 
  left_join(., GOslim, by = "GO") #Annotate genes
head(annotmaleSelectedVariables2X) #Confirm annotation
```

```{r}
annotmaleSelectedVariables2X %>%
  select(., c(gene, value.var, product)) %>%
  separate(., product, into = c("product", "transcriptVar"), sep = "%2C ") %>%
  select(., -transcriptVar) %>%
  distinct(.)
```

```{r}
annotmaleSelectedVariables2XP <- annotmaleSelectedVariables2X %>%
  dplyr::filter(., process == "P")
as.data.frame(table(annotmaleSelectedVariables2XP$GOslim))
```

```{r}
annotmaleSelectedVariables2XF <- annotmaleSelectedVariables2X %>%
  dplyr::filter(., process == "F")
as.data.frame(table(annotmaleSelectedVariables2XF$GOslim))
```

```{r}
annotmaleSelectedVariables2Y <- MaleSelectedVariables2$Y %>%
  as.data.frame(.) %>%
  rename(., gene = name) %>%
  left_join(., mRNATrack, by = "gene") %>% 
  left_join(., GOterms, by = "transcript") %>% 
  left_join(., GOslim, by = "GO") #Annotate genes
head(annotmaleSelectedVariables2Y) #Confirm annotation
```

```{r}
annotmaleSelectedVariables2Y %>%
  select(., c(gene, value.var, product)) %>%
  separate(., product, into = c("product", "transcriptVar"), sep = "%2C ") %>%
  select(., -transcriptVar) %>%
  distinct(.)
```

```{r}
annotmaleSelectedVariables2YP <- annotmaleSelectedVariables2Y %>%
  dplyr::filter(., process == "P")
as.data.frame(table(annotmaleSelectedVariables2YP$GOslim))
```

```{r}
annotmaleSelectedVariables2YF <- annotmaleSelectedVariables2Y %>%
  dplyr::filter(., process == "F")
as.data.frame(table(annotmaleSelectedVariables2YF$GOslim))
```

#### Create table with loading values and annotations

```{r}
rbind(as.data.frame(MaleSelectedVariables1$X$value),
      as.data.frame(MaleSelectedVariables1$Y$value),
      as.data.frame(MaleSelectedVariables2$X$value),
      as.data.frame(MaleSelectedVariables2$Y$value)) %>% 
  mutate(., block = c(rep("X", times = 25),
                      rep("Y", times = 25),
                      rep("X", times = 25),
                      rep("Y", times = 25))) %>%
  mutate(., component = c(rep("1", times = 50),
                          rep("2", times = 50))) %>%
  rownames_to_column(., var = "gene") %>%
  left_join(., mRNATrack, by = "gene") %>% 
  left_join(., GOterms, by = "transcript") %>% 
  left_join(., GOslim, by = "GO") %>%
  arrange(., component, block, value.var, process, transcript, GO) %>%
  distinct(.) %>%
  dplyr::group_by(gene, value.var, block) %>%
  mutate(transcript = paste(transcript, collapse = ";")) %>%
  mutate(product = paste(product, collapse = ";")) %>%
  mutate(GO = paste(GO, collapse = ";")) %>%
  mutate(GOterm = paste(GOterm, collapse = ";")) %>%
  mutate(GOslim = paste(GOslim, collapse = ";")) %>%
  mutate(process = paste(process, collapse = ";")) %>%
  ungroup(.) %>%
  distinct(.) %>%
  mutate(transcript = lapply(strsplit(transcript, split = ";"), unique)) %>%
  mutate(transcript = sapply(transcript, toString)) %>%
  mutate(transcript = gsub(transcript, pattern = ", ", replacement = ";")) %>%
  mutate(transcript = na_if(x = transcript, y = "NA")) %>%
  mutate(product = lapply(strsplit(product, split = ";"), unique)) %>%
  mutate(product = sapply(product, toString)) %>%
  mutate(product = gsub(product, pattern = ", ", replacement = ";")) %>%
  mutate(product = na_if(x = product, y = "NA")) %>%
  mutate(GO = lapply(strsplit(GO, split = ";"), unique)) %>%
  mutate(GO = sapply(GO, toString)) %>%
  mutate(GO = gsub(GO, pattern = ", ", replacement = ";")) %>%
  mutate(GO = na_if(x = GO, y = "NA")) %>%
  mutate(GOterm = lapply(strsplit(GOterm, split = ";"), unique)) %>%
  mutate(GOterm = sapply(GOterm, toString)) %>%
  mutate(GOterm = gsub(GOterm, pattern = ", ", replacement = ";")) %>%
  mutate(GOterm = na_if(x = GOterm, y = "NA")) %>%
  mutate(GOslim = lapply(strsplit(GOslim, split = ";"), unique)) %>%
  mutate(GOslim = sapply(GOslim, toString)) %>%
  mutate(GOslim = gsub(GOslim, pattern = ", ", replacement = ";")) %>%
  mutate(GOslim = na_if(x = GOslim, y = "NA")) %>%
  mutate(process = lapply(strsplit(process, split = ";"), unique)) %>%
  mutate(process = sapply(process, toString)) %>%
  mutate(process = gsub(process, pattern = ", ", replacement = ";")) %>%
  mutate(process = na_if(x = process, y = "NA")) %>%
  distinct(.) %>%
  write.table("selectedVariables-Annot-Male.txt", quote = FALSE, sep = "\t", row.names = FALSE) #Combine driver gene information for both X and Y datasets for the first component. Add block (X or Y) and component (1 or 2) information for gene methylation or expression datasets. Join with mRNA, GO term, and GO Slim information for annotations. Arrange by component, block, loading value, process, transcript, and GO, then retain only distinct entries. Group by gene, loading values, and block to collapse rows across all descriptive columns such that there is only one row per driver gene. Ungroup and retain only distinct rows. Use a series of apply commands to split strings in each column by ; then remove any duplicate values. Replace commas with semicolon delimiters and make sure NAs are identified. Retain distinct rows. Save as a txt file
```

```{bash}
head selectedVariables-Annot-Male.txt
```


# THINGS TO DO BEFORE YOU DO ANYTHING YOU FOOL:

- perf error still present even after rm samples w nzv + using dev package
- remove outliers from arrow plot? check mapping stats, library size, etc. 
- why is there much lower # fpkm genes for females? did i filter too much?
- tune sex-specific parameters
- run sex-specific spls
- extract VIPs from spls: loading value > 1
- add categorial isoform information
- add numerical CV information
