# PCA plot significance

I have made my PCA plots for each hour of the differential expression analysis (Figure below).

I want to add statistical tests and significance to support the differential gene expression at each time point (suggestion from Dr. Kevin Wong - thank you!).

Kevin sent me his Physiology Analysis Rmd for a paper he’s currently working on, which has this PCA in it with PERMANOVAs included:

I want to also do a PERMANOVA for my PCA plots.

What I’ve got so far:

## Example of DDS object creation for hour 0 (control vs. wounded):

```
countdata_0 <- countsmatrix %>% dplyr::select(`0301C1`:`0306Z`)
metadata_0 = data.frame(sample=colnames(countdata_0),
condition = stringr::str_detect(pattern = ".*C.*",string = colnames(countdata_0)),
hour = stringr::str_replace(pattern = "(.).*",replacement="\\1",string = colnames(countdata_0)),
id = stringr::str_replace(pattern=".(...).*",replacement="\\1",string=colnames(countdata_0)))
#changing TRUE and FALSE for condition to Control and Wounded
metadata_0$condition[str_detect(metadata_0$condition,"TRUE")] <- "Control"
metadata_0$condition[str_detect(metadata_0$condition,"FALSE")] <- "Wounded"
metadata_0 <- metadata_0 %>% column_to_rownames("sample")
metadata_0$condition <- as.factor(metadata_0$condition)
#DESeq2 from a counts matrix; requires count data, metadata, and the experimental design
dds_0=DESeq2::DESeqDataSetFromMatrix(countData = countdata_0, colData = metadata_0, design = ~condition)
dds_0 <- estimateSizeFactors(dds_0)
#prefiltering recommended by DESeq2 Vignette (removes anything with reads below 10)
keep_0 <- rowSums(counts(dds_0)) >= 10
dds_0 <- dds_0[keep_0,]
dds_0<-DESeq(dds_0)
res_h0<-results(dds_0)
summary(res_h0,alpha=0.05)
#5 downregulated genes (3 outliers), no upregulated genes
#identifying significant genes
resSig_h0<-res_h0[which(res_h0$padj<0.05), ]
#annotate results to include Gene Function
anno_h0<-merge(as.data.frame(resSig_h0),Pdam_Gene_Names_Info,by='row.names',all=TRUE)
anno_h0<-na.omit(anno_h0)
```

## Principal component analysis plot and Scree plots

```
#need to transform the data for plotting
dds_vst0<- vst(dds_0,blind=FALSE)
plotPCA(dds_vst0)
#plotPCA is one function to do it, or you can use the "prcomp" function, which lets you create a scree plot
pca_h0 <- prcomp(t(assay(dds_vst0)))
fviz_eig(pca_h0)
```

## PCA figures for manuscript

```
#plotting the PCA in ggplot
pca12_0 <- plotPCA(dds_vst0,intgroup=c("condition"),returnData = TRUE)
pca_0 = ggplot(pca12_0, aes(PC1,PC2,shape=condition,color=condition)) +
geom_point(size=3) +
xlab(paste0("PC1 (49%)")) +
ylab(paste0("PC2 (24%)")) +
theme(legend.position="right") +
theme(text = element_text(size=12)) +
theme(legend.key.size = unit(0.5, "cm")) +
theme(legend.title=element_text(size=12)) +
scale_shape_discrete(name="Condition") +
scale_color_discrete(name="Condition") +
theme_classic(base_size=12,base_family="serif")
```

Now I’m having trouble with the PERMANOVA part.

## This is the section of code Kevin used:

```
master <- join_all(dfs, by=c("Fragment.ID", "Day", "Group"))
master.pca <- master
coral_info<-master.pca[c(2,3)] #columns 2 and 3 are "Day" and "Group" in his metadata file
#Examine PERMANOVA results.
# scale data
vegan <- scale(master.pca[c(4:18)]) #these columns he created as calculations for physiology metrics (i.e. ChlC2.ugcell, Carb.mgcell, Protein.mgcell)
# PerMANOVA
permanova<-adonis(vegan ~ Day*Group, data = master.pca, method='eu')
z_pca<-permanova$aov.tab
z_pca
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Day 2 174.38 87.189 9.0242 0.26421 0.001 ***
## Group 2 78.41 39.205 4.0577 0.11880 0.002 **
## Day:Group 4 59.39 14.848 1.5368 0.08999 0.069 .
## Residuals 36 347.82 9.662 0.52700
## Total 44 660.00 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

I think I need to make a dataset like this

Sample ID | Condition | Pdam0000001 |
---|---|---|

301 | Control | 6 |

302 | Control | 18 |

303 | Control | 12 |

304 | Wounded | 12 |

305 | Wounded | 20 |

306 | Wounded | 43 |

Ok, when I rearranged the data to look like that, the permanova didn’t work.

# PERMANOVA

```
#first need to restructure results matrix so that it includes each row as a sample with the condition, and then each column is a gene
hour0_countsmatrix <- assay(dds_0)
t(hour0_countsmatrix) -> hour0_countsmatrix
#we use the dds object instead of the original countdata_0 matrix because this has filtered out genes with low counts (<10)
PCA.h0.countsdata <- merge(metadata_0, hour0_countsmatrix, by='row.names', all=TRUE)
dim(PCA.h0.countsdata)
# scale data
vegan <- scale(PCA.h0.countsdata[c(5:19451)]) #we just want to scale the gene counts
# PERMANOVA
permanova<-adonis2(vegan ~ condition, data = PCA.h0.countsdata, method='eu', na.rm = TRUE, nperm = 999)
```

So I think maybe it has to be formatted a different way.

Sample ID | Condition | Gene | Count |
---|---|---|---|

301 | Control | Pdam0001 | 6 |

302 | Control | Pdam0001 | 18 |

303 | Control | Pdam0001 | 12 |

304 | Wounded | Pdam0001 | 12 |

305 | Wounded | Pdam0001 | 20 |

306 | Wounded | Pdam0001 | 43 |

And then the design for the PERMANOVA would actually be count ~ gene*condition

Let’s try that:

```
#first need to restructure results matrix so that it includes each row as a sample with the condition, and then each column is a gene
#we use the dds object instead of the original countdata_0 matrix because this has filtered out genes with low counts (<10)
hour0_countsmatrix <- assay(dds_0)
head(hour0_countsmatrix) #first gene is pdam_00021773
tail(hour0_countsmatrix) #last gene is pdam_00025493
t(hour0_countsmatrix) %>% as.data.frame() -> hour0_countsmatrix
PCA.h0.countsdata <- merge(metadata_0, hour0_countsmatrix, by='row.names', all=TRUE)
PCA.h0.countsdata %>%
pivot_longer(cols = pdam_00021773:pdam_00025493, names_to = "Gene ID", values_to = "Count") %>% as.matrix()->PCA.h0.countsdata_longformat
# PERMANOVA
permanova<-adonis2(Count ~ `Gene ID`*condition, data = PCA.h0.countsdata_longformat, method='eu')
```

Doesn’t work. Error:Error in eval(YVAR, environment(formula), globalenv()) : object ‘Count’ not found

This is what the long format data frame looks like:

I’m so confused. I thought maybe i would want Gene and Condition and Gene:Condition as the terms for the PERMANOVA. I need to look up adonis2 and what it accepts in the formula.

From adonis CRAN file “The left-hand side (LHS) of the formula must be either a community data matrix or a dissimilarity matrix, e.g., from vegdist or dist.”

So that’s why Kevin used “vegan” as his input variable in the formula, which is a matrix of just his dependent variables.

Adonis/Vegan is essentially running some sort of dissimilarity assessment, so it needs a matrix to run it on. So you can’t put all the counts into one column.

I think I’m back to square one, where I actually got a result (133-150). I think it worked but it didn’t give me significant results. I ran it for each hour as well, here is that code and screenshots of the results below.

# PERMANOVA for hour 0

```
#we use the dds object instead of the original countdata_0 matrix because this has filtered out genes with low counts (<10)
hour0_countsmatrix <- assay(dds_0)
t(hour0_countsmatrix) -> hour0_countsmatrix
PCA.h0.countsdata <- merge(metadata_0, hour0_countsmatrix, by='row.names', all=TRUE)
# scale data
vegan <- scale(PCA.h0.countsdata[c(5:19451)]) #we just want to scale the gene counts
# PERMANOVA
permanova<-adonis2(vegan ~ condition, data = PCA.h0.countsdata, method='eu', na.rm=TRUE, nperm = 999)
permanova
```

# PERMANOVA for hour 1

```
#we use the dds object instead of the original countdata_0 matrix because this has filtered out genes with low counts (<10)
hour1_countsmatrix <- assay(dds_1)
t(hour1_countsmatrix) -> hour1_countsmatrix
PCA.h1.countsdata <- merge(metadata_1, hour1_countsmatrix, by='row.names', all=TRUE)
# scale data
vegan <- scale(PCA.h1.countsdata[c(5:19536)]) #we just want to scale the gene counts
# PERMANOVA
permanova<-adonis2(vegan ~ condition, data = PCA.h1.countsdata, method='eu', na.rm=TRUE, nperm = 999)
permanova
```

# PERMANOVA for hour 2

```
#we use the dds object instead of the original countdata_0 matrix because this has filtered out genes with low counts (<10)
hour2_countsmatrix <- assay(dds_2)
t(hour2_countsmatrix) -> hour2_countsmatrix
PCA.h2.countsdata <- merge(metadata_2, hour2_countsmatrix, by='row.names', all=TRUE)
# scale data
vegan <- scale(PCA.h2.countsdata[c(5:19779)]) #we just want to scale the gene counts
# PERMANOVA
permanova<-adonis2(vegan ~ condition, data = PCA.h2.countsdata, method='eu', na.rm=TRUE, nperm = 999)
permanova
```

# PERMANOVA for hour 4

```
#we use the dds object instead of the original countdata_0 matrix because this has filtered out genes with low counts (<10)
hour4_countsmatrix <- assay(dds_4)
t(hour4_countsmatrix) -> hour4_countsmatrix
PCA.h4.countsdata <- merge(metadata_4, hour4_countsmatrix, by='row.names', all=TRUE)
# scale data
vegan <- scale(PCA.h4.countsdata[c(5:19773)]) #we just want to scale the gene counts
# PERMANOVA
permanova<-adonis2(vegan ~ condition, data = PCA.h4.countsdata, method='eu', na.rm=TRUE, nperm = 999)
permanova
```

I also tried running it without the scale() function for hour 4 and I got a less significant result.

I’m going to abandon this analysis for the time being because I don’t think this is super necessary for the PCA plots anyways.