---
output:
  word_document: default
  pdf_document: default
  html_document: default
---
#---
title: "Supplemental Data 2 - statistical analysis file for neuronal somata data in Davies et al: APOE isoform does not influence trans-synaptic spread of tau pathology in a mouse model"
output: html_document
#---

### Packages
Load packages required for analysis
```{r}
library(readr)
library(table1)
library(ggplot2)
library(ggResidpanel)
library(RColorBrewer)
library(tidyverse)
library(dplyr)
library(car)
library(FSA)
library(rcompanion)
library(forecast)
library(gridExtra)
library(corrplot)
library(lmerTest)
library(emmeans)
library(ggforce)
```


### Set working directory and import dataframe 
Dataframe columns are assigned as factors (f)/integers (i)/double (d) etc
```{r setup, include=FALSE}
# Set working directory (wd) where your spreadsheet is stored.
#setwd("~/path")

# Load in dataframe (uses "readr")
MouseFF <- read_csv("SuppData1_Davies_etal_SomaData.csv", col_types = cols (
  BlindedID = 'f', 
  UnblindedID = "f", 
  Genotype = 'f',
  Sex = 'f',
  AgeGroup = "f", 
  Age_Inj = 'd',
  Age_Cull = 'd',
  Incubation = 'd',
  Brain_Region = 'f',
  TotalGFP = "d",
  TotalGFP_DAPI = "d",
  TotalTau13 = "d",
  TotalGFP_Tau13 = "d",
  TotalTau13_GFP = "d",
  TotalGFP_Tau13_DAPI = "d",
  TotalTau13_GFP_DAPI = "d",
  PercentGFP_DAPI_ColocTau13 = "d",
  PS_GFP = "d",
  PS_GFP_DAPI = "d",
  PS_Tau13 = "d",
  PS_GFP_Tau13 = "d",
  PS_Tau13_GFP = "d",
  PS_GFP_Tau13_DAPI = "d",
  PS_Tau13_GFP_DAPI = "d",
  TotalRecipient = "d",
  PS_Recipient = "d",
  RecipientNormGFP_DAPI = "d",
  RecipientNormGFP_DAPI_100 = "d",
  NumDonor_1Recipient = 'd',
  RecipientNormGFP_Tau13_DAPI = "d",
  RecipientNormGFP_Tau13_DAPI_100 = "d"))

# Are some NAs in dataframe - in the DonorCellTau13 DonorCellALZ50 DonorCellTau13ALZ50 columns. These are necessary - will have to remember to remove them when analyse this part. 
```


```{r}
# Check structure of dataframe
str(MouseFF)

# Create a genotype ordered variable (as not sure how R will plot them)
MouseFF$GenotypeOrder <- ordered(MouseFF$Genotype, levels=c("APOEKO", "APOE2", "APOE3", "APOE4"))

# Create a age ordered variable (as not sure how R will plot them) - in correct order
levels(MouseFF$AgeGroup)

# Concatenate Genotype and AgeGroup - to use to determine how bars are filled in plots later
MouseFF$GenotypeAge <- paste(MouseFF$Genotype, MouseFF$AgeGroup)
MouseFF$GenotypeAge <- as.factor(MouseFF$GenotypeAge)

# Order new GenotypeAgeGroup so that bars are filled in the correct order
MouseFF$GenotypeAgeOrder <- ordered(MouseFF$GenotypeAge, levels=c("APOEKO Adult", "APOE2 Adult", "APOE3 Adult", "APOE4 Adult", "APOEKO Aged", "APOE2 Aged", "APOE3 Aged", "APOE4 Aged"))
```


### Remove 81TSJ
This mouse was left for longer than others due to lockdown - needs to be removed
```{r}
# Subset based BlindedID
MouseFF2 <- subset(MouseFF, BlindedID!="81TSJ") 
```


### Remove mice that do not have enough donor cells
Some mice have very few donor (GFP+/DAPI+) cells. 

Used spreadsheet 'Donor cell number to identify 1 recipient cell' (based on recipientcell summary table)
For each group, calculated the average number of donro cells required to identify 1 recipient cell.

APOEKO Adult is the highest at 219.298246
APOE2 Adult have the lowest at 49

Could also calculate the mean number of recipient for 1 donor for whole dataset = 93 donor cells required to see 1 recipient cell

Could also apply group specific cut offs? --> 22/94 removed - But seriously effects n of APOE2 Adult
If I am normalising to the number of GFP cells then is this accounted for? i.e is it necessary to do this? 

Decided to go with the least conservative estimate- remove all mice with < 49 donor cells

```{r}
# Create 'cut' dataframe. Subset based on number of donor cells (GFP_DAPI) in TOTAL measurement - if apply it to whole dataframe will remove regions that have low number of donor cells 

# Work out which mice we want to remove = 7 mice removed
MiceToRemoveKOAdult <- subset(MouseFF2, Brain_Region=="Total" & Genotype =="APOEKO" & AgeGroup =="Adult" & TotalGFP_DAPI < 219.298246) # would remove 9/13
MiceToRemove2Adult <- subset(MouseFF2, Brain_Region=="Total" & Genotype =="APOE2" & AgeGroup =="Adult" & TotalGFP_DAPI < 49.7512437810945) # would remove 4/13
MiceToRemove3Adult <- subset(MouseFF2, Brain_Region=="Total" & Genotype =="APOE3" & AgeGroup =="Adult" & TotalGFP_DAPI < 124.069478908189) # would remove 4/14
MiceToRemove4Adult <- subset(MouseFF2, Brain_Region=="Total" & Genotype =="APOE4" & AgeGroup =="Adult" & TotalGFP_DAPI < 83.3333333333333) # would remove 0/12

MiceToRemoveKOAged <- subset(MouseFF2, Brain_Region=="Total" & Genotype =="APOEKO" & AgeGroup =="Aged" & TotalGFP_DAPI < 150.15015015015) # would remove 0/6
MiceToRemove2Aged <- subset(MouseFF2, Brain_Region=="Total" & Genotype =="APOE2" & AgeGroup =="Aged" & TotalGFP_DAPI < 118.623962040332) # would remove 2/10
MiceToRemove3Aged <- subset(MouseFF2, Brain_Region=="Total" & Genotype =="APOE3" & AgeGroup =="Aged" & TotalGFP_DAPI < 110.987791342952) # would remove 0/14
MiceToRemove4Aged <- subset(MouseFF2, Brain_Region=="Total" & Genotype =="APOE4" & AgeGroup =="Aged" & TotalGFP_DAPI < 67.5675675675676) # would remove 3/12


MiceToRemove <- subset(MouseFF, TotalGFP_DAPI < 93.5103796521414) 


MouseFFTotal <- subset(MouseFF2, Brain_Region=="Total")
MiceToRemove <- subset(MouseFF, TotalGFP_DAPI < 49) # Use this one

MouseFF3 <- subset(MouseFF2, !(BlindedID %in% c("14", "88", "29", "28", "29", "103", "31", "32", "34", "45", "46", "53", "75")))
```

### Replace NAs with 0
```{r}
MouseFF3[is.na(MouseFF3)] <- 0
```

### Replace negative values with 0 (artefact of segmentation)
```{r}
MouseFF3[MouseFF3<0] <- 0
```

### Split dataframe into separate dataframes
One with just the total measurements from each mouse, the second with the measurements from the EC (Super & Deep) split by brain region. These will be needed for different parts of analyses. 
```{r}
# Subset based on brain region variabe - Total (all regions summed)
MouseFFTotal <- subset(MouseFF3, Brain_Region=="Total") 
write_csv(MouseFFTotal, "MouseFFTotal.csv")


# Subset based on brain region variabe - ECTotal (ECSuper and ECDeep summed)
MouseFFECTotal <- subset(MouseFF3, Brain_Region=="ECTotal") 

# Subset based on brain region variabe - Regions (Total and ECTotal removed)
MouseFFRegions <- subset(MouseFF3, Brain_Region!="Total") 
MouseFFRegions <- subset(MouseFFRegions, Brain_Region!="ECTotal") 

# Subset based on brain region variabe - Regions (Total and ECTotal removed)
MouseFFECTotalRegions <- subset(MouseFF3, Brain_Region!="Total") 
MouseFFECTotalRegions <- subset(MouseFFECTotalRegions, Brain_Region!="ECDeep") 
MouseFFECTotalRegions <- subset(MouseFFECTotalRegions, Brain_Region!="ECSuper") 

  
# Even after subsetting (based on Region factor levels) to get split dataframes, the factor levels are retained in the dataframe:
str(MouseFFRegions)
levels(MouseFFRegions$Brain_Region) # "CA1"     "CA23"    "DG"      "EC"      "ECSuper" "ECDeep"  "Total"

# To fix, just apply factor() to variable again after subsetting
MouseFFTotal$Brain_Region <- factor(MouseFFTotal$Brain_Region)
MouseFFECTotal$Brain_Region <- factor(MouseFFECTotal$Brain_Region)
MouseFFRegions$Brain_Region <- factor(MouseFFRegions$Brain_Region)
MouseFFECTotalRegions$Brain_Region <- factor(MouseFFECTotalRegions$Brain_Region)

levels(MouseFFTotal$Brain_Region)
levels(MouseFFECTotal$Brain_Region)
levels(MouseFFRegions$Brain_Region)
levels(MouseFFECTotalRegions$Brain_Region)

# Order brain regions so that bars are filled in the correct order
MouseFFRegions$Brain_RegionOrder <- ordered(MouseFFRegions$Brain_Region, levels=c("ECSuper", "ECDeep", "DG", "CA23", "CA1"))
levels(MouseFFRegions$Brain_RegionOrder)

MouseFFECTotalRegions$Brain_RegionOrder <- ordered(MouseFFECTotalRegions$Brain_Region, levels=c("ECTotal", "DG", "CA23", "CA1"))
levels(MouseFFECTotalRegions$Brain_RegionOrder)
```

### Summary table of Mouse info
Use the MouseFFTotal dataframe as this has one row per mouse
```{r}
# Can start from here and load "clean" dataset after exlusion critera applied
#read_csv("MouseFFTotal.csv")
# Mouse info summary
# Assign headers that you want in your table to the variables you want in your table - can include various other things
label(MouseFFTotal$AgeGroup) <- "Age Group"
label(MouseFFTotal$Age_Inj) <- "Age at Injection (months)"
label(MouseFFTotal$Age_Cull) <- "Age at Cull (months)"
label(MouseFFTotal$Incubation) <- "Length of incubation (months)"
label(MouseFFTotal$TotalGFP) <- "Total number of donor cells (GFP+)"
label(MouseFFTotal$TotalGFP_DAPI) <- "Total number of donor cells (GFP+/DAPI+)"
label(MouseFFTotal$TotalGFP_Tau13_DAPI) <- "Total number of donor cells (GFP+/DAPI+/Tau13+)"
label(MouseFFTotal$PercentGFP_DAPI_ColocTau13) <- "Percent of donor cells (GFP+/DAPI+) positive for Tau13 (%)"

label(MouseFFTotal$PS_GFP) <- "Average number of donor cells (GFP+) per section"
label(MouseFFTotal$PS_GFP_DAPI) <- "Average number of donor cells (GFP+/DAPI+) per section"
label(MouseFFTotal$PS_GFP_Tau13_DAPI) <- "Average number of donor cells (GFP+/DAPI+/Tau13+) per section"

label(MouseFFTotal$TotalRecipient) <- "Total number of recipient cells (GFP-/Tau13+)"
label(MouseFFTotal$PS_Recipient) <- "Average number of recipient cells (GFP-/Tau13+) per section"
label(MouseFFTotal$RecipientNormGFP_DAPI_100) <- "Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+) (%)"
label(MouseFFTotal$RecipientNormGFP_Tau13_DAPI_100) <- "Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+/Tau13+) (%)"
label(MouseFFTotal$NumDonor_1Recipient) <- "Number of donor cells (GFP+/DAPI+) to see 1 recipient cell (GFP-/DAPI+)"

# Then plot the table
info1 <- table1(~ Sex + Age_Inj + Age_Cull + Incubation | AgeGroup + GenotypeOrder, data = MouseFFTotal)
info1

info2 <- table1(~ Sex + Age_Inj + Age_Cull + Incubation | GenotypeOrder + AgeGroup, data = MouseFFTotal)
info2

donorcell <- table1(~ TotalGFP + TotalGFP_DAPI + TotalGFP_Tau13_DAPI + PS_GFP + PS_GFP_DAPI + PS_GFP_Tau13_DAPI | AgeGroup + GenotypeOrder, data = MouseFFTotal)
donorcell

donorcellcoloc <- table1(~ PercentGFP_DAPI_ColocTau13 | AgeGroup + GenotypeOrder , data = MouseFFTotal)
donorcellcoloc

recipientcells <- table1(~ TotalRecipient + PS_Recipient + RecipientNormGFP_DAPI_100 + RecipientNormGFP_Tau13_DAPI_100 | AgeGroup + GenotypeOrder, data = MouseFFTotal)
recipientcells
```



### Check whether groups are age-matched
### Age at injection
```{r}
# Plot age at injection data - both ages

ggplot(MouseFFTotal, aes(x = GenotypeOrder, y = Age_Inj, fill = GenotypeOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Age at Injection") + # title of the plot
  scale_y_continuous("Age at Injection (months)", breaks = c(5,6,7,8,9,10,11,12,13,14,15,16), limits = c(5,16.5)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  # scale_fill_brewer(palette = "Blues") + 
  scale_fill_manual(values = c("#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)

#"#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"

ggsave("1_Age at Injection_All.png", width = 30, height = 20, units = "cm", dpi = 300)



# Plot age at injection data - adult only
ggplot(subset(MouseFFTotal, AgeGroup %in% c("Adult")), aes(x = GenotypeOrder, y = Age_Inj, fill = GenotypeOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Age at Injection - Adult") + # title of the plot
  scale_y_continuous("Age at Injection (months)", limits = c(5,6.1)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Blues")
  scale_fill_manual(values = c("#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"))  

ggsave("1_Age at Injection_Adult.png", width = 30, height = 20, units = "cm", dpi = 300)



# Plot age at injection data - aged only
ggplot(subset(MouseFFTotal, AgeGroup %in% c("Aged")), aes(x = GenotypeOrder, y = Age_Inj, fill = GenotypeOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Age at Injection - Aged") + # title of the plot
  scale_y_continuous("Age at Injection (months)", limits = c(15,15.8)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  # scale_fill_brewer(palette = "Blues")
  scale_fill_manual(values = c("#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"))

ggsave("1_Age at Injection_Aged.png", width = 30, height = 20, units = "cm", dpi = 300)
```

```{r}
# Run *two-way ANOVA* to assess whether there is an age difference - two-way this time - is this more useful? Can't get data to meet assumptions, so can't use.
AgeInjAOV <- aov(Age_Inj ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(AgeInjAOV)

# Check that meets the assumptions
# Check for normality of residuals
hist(resid(AgeInjAOV))
shapiro.test(resid(AgeInjAOV)) # p-value = 0.006881

# Levene's Test for Homogeneity of Variance (center = median)
leveneTest(Age_Inj ~ Genotype * AgeGroup, data = MouseFFTotal) # p-value = 0.5473

# Could also run a linear model and assess residuals
AgeInjLM <- lm(Age_Inj ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(AgeInjLM)

resid_panel(AgeInjLM) # Look okay

shapiro.test(resid(AgeInjLM)) # p-value = 0.02877

anova(AgeInjLM)
emmeans(AgeInjLM, pairwise ~ Genotype)

# Significant main effect of Age Group and Genotype
# But no significant effects in post-hoc test

# Tukey transform the data to see if helps with assumptions (uses rcompanion)
MouseFFTotal$Age_Inj_tukey = 
  transformTukey(MouseFFTotal$Age_Inj,
                 plotit=FALSE)

AgeInjAOV_tukey <- aov(Age_Inj_tukey ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(AgeInjAOV_tukey)

shapiro.test(resid(AgeInjAOV_tukey)) # p-value = 9.252e-11
leveneTest(Age_Inj_tukey ~ Genotype * AgeGroup, data = MouseFFTotal) # p-value = 0.0004054


# Try linear model
AgeInjLM_tukey <- lm(Age_Inj_tukey ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(AgeInjLM_tukey)

resid_panel(AgeInjLM_tukey) # Looks even worse

shapiro.test(resid(AgeInjLM_tukey)) # p-value = 1.234e-09

# Box-Cox transformation - log likelihood method
lam.l = BoxCox.lambda(MouseFFTotal$Age_Inj, method = "loglik") 
MouseFFTotal$Age_Inj_BCL = BoxCox(MouseFFTotal$Age_Inj, lam.l)

AgeInjAOV_BCL <- aov(Age_Inj_BCL ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(AgeInjAOV_BCL)

shapiro.test(resid(AgeInjAOV_BCL)) # p-value = 1.586e-07
leveneTest(Age_Inj_BCL ~ Genotype * AgeGroup, data = MouseFFTotal) # p-value = 0.0002706 

# Try linear model
AgeInjLM_BCL <- lm(Age_Inj_BCL ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(AgeInjLM_BCL)

resid_panel(AgeInjLM_BCL) # Looks even worse

shapiro.test(resid(AgeInjLM_BCL)) # p-value = 3.493e-06


# Box-Cox transformation - guerrero method
lam.g = BoxCox.lambda(MouseFFTotal$Age_Inj, method = "guerrero") 
MouseFFTotal$Age_Inj_BCG = BoxCox(MouseFFTotal$Age_Inj, lam.g)

AgeInjAOV_BCG <- aov(Age_Inj_BCG ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(AgeInjAOV_BCG)

shapiro.test(resid(AgeInjAOV_BCG)) # p-value = 3.428e-08
leveneTest(Age_Inj_BCG ~ Genotype * AgeGroup, data = MouseFFTotal) # p-value = 0.0001011

# Try linear model
AgeInjLM_BCG <- lm(Age_Inj_BCG ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(AgeInjLM_BCG)

resid_panel(AgeInjLM_BCG) # Looks even worse

shapiro.test(resid(AgeInjLM_BCG)) # p-value = 1.838e-07

# Transformations don't help with meeting assumptions of 2-way ANOVA
# LM maybe okay - but main effects and post-hoc don't make sense
```

```{r}
# Run one-way ANOVA to assess whether there is an age difference on each independent group (Adult and Aged)

# Adult
AgeInjAOVAdult <- aov(Age_Inj ~ Genotype, data = (subset(MouseFFTotal, AgeGroup %in% c("Adult"))))

summary(AgeInjAOVAdult)

# Check that meets the assumptions
# Check for normality of residuals
shapiro.test(resid(AgeInjAOVAdult)) # p-value = 0.04838

# Levene's Test for Homogeneity of Variance (center = median)
leveneTest(Age_Inj ~ Genotype, data = (subset(MouseFFTotal, AgeGroup %in% c("Adult")))) # p-value = 0.8084


# Residuals are not normally distributed so will use a Kruskall-Wallis 
AgeInjKWAdult <- kruskal.test(Age_Inj ~ Genotype, data = (subset(MouseFFTotal, AgeGroup %in% c("Adult"))))
AgeInjKWAdult # p-value = 0.1069


# No significant difference between genotyped in adult group

# Aged
AgeInjAOVAged <- aov(Age_Inj ~ Genotype, data = (subset(MouseFFTotal, AgeGroup %in% c("Aged"))))

summary(AgeInjAOVAged)

# Check that meets the assumptions
# Check for normality of residuals
shapiro.test(resid(AgeInjAOVAged)) # p-value = 1.149e-06

# Levene's Test for Homogeneity of Variance (center = median)
leveneTest(Age_Inj ~ Genotype, data = (subset(MouseFFTotal, AgeGroup %in% c("Aged")))) # p-value = 0.8522


# Residuals are not normally distributed so will use a Kruskall-Wallis 
AgeInjKWAged <- kruskal.test(Age_Inj ~ Genotype, data = (subset(MouseFFTotal, AgeGroup %in% c("Aged"))))
AgeInjKWAged # p-value = 0.7886

# No significant difference between gebotypes in aged group
```


### Age at cull
```{r}
# Plot age at cull data - both ages

ggplot(MouseFFTotal, aes(x = GenotypeOrder, y = Age_Cull, fill = GenotypeOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Age at Cull") + # title of the plot
  scale_y_continuous("Age at Cull (months)", breaks = c(8,9,10,11,12,13,14,15,16,17,18, 19), limits = c(8,19.5)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  # scale_fill_brewer(palette = "Blues") + 
  scale_fill_manual(values = c("#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)

#"#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"

ggsave("2_Age at Cull_All.png", width = 30, height = 20, units = "cm", dpi = 300)


# Plot age at cull data - adult only
ggplot(subset(MouseFFTotal, AgeGroup %in% c("Adult")), aes(x = GenotypeOrder, y = Age_Cull, fill = GenotypeOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Age at Cull - Adult") + # title of the plot
  scale_y_continuous("Age at Cull (months)", limits = c(8.2,9.3)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Blues")
  scale_fill_manual(values = c("#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"))  

ggsave("2_Age at Cull_Adult.png", width = 30, height = 20, units = "cm", dpi = 300)


# Plot age at cull data - aged only
ggplot(subset(MouseFFTotal, AgeGroup %in% c("Aged")), aes(x = GenotypeOrder, y = Age_Cull, fill = GenotypeOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Age at Cull - Aged") + # title of the plot
  scale_y_continuous("Age at Cull (months)", limits = c(18.3,19.1)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  # scale_fill_brewer(palette = "Blues")
  scale_fill_manual(values = c("#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"))

ggsave("2_Age at Cull_Aged.png", width = 30, height = 20, units = "cm", dpi = 300)
```

```{r}
# Run *two-way ANOVA* to assess whether there is an age difference - two-way this time - is this more useful? Can't get data to meet assumptions, so can't use.
AgeCullAOV <- aov(Age_Cull ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(AgeCullAOV)

# Check that meets the assumptions
# Check for normality of residuals
hist(resid(AgeCullAOV))
shapiro.test(resid(AgeCullAOV)) # p-value = 0.009657

# Levene's Test for Homogeneity of Variance (center = median)
leveneTest(Age_Cull ~ Genotype * AgeGroup, data = MouseFFTotal) # p-value = 0.3965

# Could also run a linear model and assess residuals
AgeCullLM <- lm(Age_Cull ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(AgeCullLM)

resid_panel(AgeCullLM) # Look okay

shapiro.test(resid(AgeCullLM)) # p-value = 0.009657

anova(AgeCullLM)
emmeans(AgeCullLM, pairwise ~ Genotype)

# Significant main effect of Age Group and Genotype
# But no significant effects in post-hoc test

# Tukey transform the data to see if helps with assumptions (uses rcompanion)
MouseFFTotal$Age_Cull_tukey = 
  transformTukey(MouseFFTotal$Age_Cull,
                 plotit=FALSE)

AgeCullAOV_tukey <- aov(Age_Cull_tukey ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(AgeInjAOV_tukey)

shapiro.test(resid(AgeCullAOV_tukey)) # p-value = 1.54e-10
leveneTest(Age_Cull_tukey ~ Genotype * AgeGroup, data = MouseFFTotal) # p-value = 0.000156

# Box-Cox transformation - log likelihood method
lam.l = BoxCox.lambda(MouseFFTotal$Age_Cull, method = "loglik") 
MouseFFTotal$Age_Cull_BCL = BoxCox(MouseFFTotal$Age_Cull, lam.l)

AgeCullAOV_BCL <- aov(Age_Cull_BCL ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(AgeCullAOV_BCL)

shapiro.test(resid(AgeCullAOV_BCL)) # p-value = 7.55e-07
leveneTest(Age_Cull_BCL ~ Genotype * AgeGroup, data = MouseFFTotal) # p-value = 0.0003983  


# Box-Cox transformation - guerrero method
lam.g = BoxCox.lambda(MouseFFTotal$Age_Cull, method = "guerrero") 
MouseFFTotal$Age_Cull_BCG = BoxCox(MouseFFTotal$Age_Cull, lam.g)

AgeCullAOV_BCG <- aov(Age_Cull_BCG ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(AgeCullAOV_BCG)

shapiro.test(resid(AgeCullAOV_BCG)) # p-value = 3.361e-07
leveneTest(Age_Cull_BCG ~ Genotype * AgeGroup, data = MouseFFTotal) # p-value = 0.0002347

# Transformations don't help with meeting assumptions of 2-way ANOVA
```

```{r}
# Run one-way ANOVA to assess whether there is an age difference on each independent group (Adult and Aged)

# Adult
AgeCullAOVAdult <- aov(Age_Cull ~ Genotype, data = (subset(MouseFFTotal, AgeGroup %in% c("Adult"))))
summary(AgeCullAOVAdult) # p-value = 0.175

# Check that meets the assumptions
# Check for normality of residuals
shapiro.test(resid(AgeCullAOVAdult)) # p-value = 0.1151

# Levene's Test for Homogeneity of Variance (center = median)
leveneTest(Age_Cull ~ Genotype, data = (subset(MouseFFTotal, AgeGroup %in% c("Adult")))) # p-value = 0.6411

# No significant difference genotypes in adult group


# Aged
AgeCullAOVAged <- aov(Age_Cull ~ Genotype, data = (subset(MouseFFTotal, AgeGroup %in% c("Aged"))))
summary(AgeCullAOVAged)

# Check that meets the assumptions
# Check for normality of residuals
shapiro.test(resid(AgeCullAOVAged)) # p-value = 3.072e-07

# Levene's Test for Homogeneity of Variance (center = median)
leveneTest(Age_Cull ~ Genotype, data = (subset(MouseFFTotal, AgeGroup %in% c("Aged")))) # p-value = 0.8779


# Residuals are not normally distributed so will use a Kruskall-Wallis 
AgeCullKWAged <- kruskal.test(Age_Cull ~ Genotype, data = (subset(MouseFFTotal, AgeGroup %in% c("Aged"))))
AgeCullKWAged # p-value = 0.9483

# No significant difference between genotypes in aged group
```


### Incubation period
```{r}
# Plot incubation period data - both ages

ggplot(MouseFFTotal, aes(x = GenotypeOrder, y = Incubation, fill = GenotypeOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Incubation period") + # title of the plot
  scale_y_continuous("Incubation period length (months)",  limits = c(3.15,3.3)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  # scale_fill_brewer(palette = "Blues") + 
  scale_fill_manual(values = c("#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)

ggsave("3_Incubation_Aged.png", width = 30, height = 20, units = "cm", dpi = 300)

#"#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"
```

```{r}
# Run *two-way ANOVA* to assess whether there is an age difference - two-way this time
IncubationAOV <- aov(Incubation ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(IncubationAOV)

# Check that meets the assumptions
# Check for normality of residuals
hist(resid(IncubationAOV))
shapiro.test(resid(IncubationAOV)) # p-value = 0.00285

# Levene's Test for Homogeneity of Variance (center = median)
leveneTest(Incubation ~ Genotype * AgeGroup, data = MouseFFTotal) # p-value = 0.7879


IncubationLM <- lm(Incubation ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(IncubationLM)

resid_panel(IncubationLM) # Look okay
anova(IncubationLM)

# No significant differences

# Tukey transform the data to see if helps with assumptions (uses rcompanion)
MouseFFTotal$Incubation_tukey = 
  transformTukey(MouseFFTotal$Incubation,
                 plotit=FALSE)

IncubationAOV_tukey <- aov(Incubation_tukey ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(IncubationAOV_tukey)

shapiro.test(resid(IncubationAOV_tukey)) # p-value = 0.01449
leveneTest(Incubation_tukey ~ Genotype * AgeGroup, data = MouseFFTotal) # p-value = 0.7706

IncubationLM_tukey <- lm(Incubation_tukey ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(IncubationLM_tukey)
anova(IncubationLM_tukey)
resid_panel(IncubationLM_tukey)

# QQ-plots Look okay - so probably okay - No significant differences

# Box-Cox transformation - log likelihood method
lam.l = BoxCox.lambda(MouseFFTotal$Incubation, method = "loglik") 
MouseFFTotal$Incubation_BCL = BoxCox(MouseFFTotal$Incubation, lam.l)

IncubationAOV_BCL <- aov(Incubation_BCL ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(IncubationAOV_BCL)

shapiro.test(resid(IncubationAOV_BCL))
leveneTest(Incubation_BCL ~ Genotype * AgeGroup, data = MouseFFTotal)  


# Box-Cox transformation - guerrero method
lam.g = BoxCox.lambda(MouseFFTotal$Incubation, method = "guerrero") 
MouseFFTotal$Incubation_BCG = BoxCox(MouseFFTotal$Incubation, lam.g)

IncubationAOV_BCG <- aov(Incubation_BCG ~ Genotype * AgeGroup, data = MouseFFTotal)
summary(IncubationAOV_BCG)

shapiro.test(resid(IncubationAOV_BCG)) 
leveneTest(Incubation_BCG ~ Genotype * AgeGroup, data = MouseFFTotal) 

# Transformations don't help with meeting assumptions of 2-way ANOVA - but residual plots look okay for LMs
# Don't want to run one ways as want to directly compare the two groups? 
```

### Set reference levels
```{r}
# Genotype as APOE3 or APOEKO?
MouseFFTotal$Genotype <- relevel(MouseFFTotal$Genotype, ref = "APOE3")
MouseFFTotal$AgeGroup <- relevel(MouseFFTotal$AgeGroup, ref = "Adult")

MouseFFRegions$Genotype <- relevel(MouseFFRegions$Genotype, ref = "APOE3")
MouseFFRegions$AgeGroup <- relevel(MouseFFRegions$AgeGroup, ref = "Adult")

MouseFFRegions$Brain_Region <- relevel(MouseFFRegions$Brain_Region, ref = "ECSuper")
MouseFFECTotalRegions$Brain_Region <- relevel(MouseFFECTotalRegions$Brain_Region, ref = "ECTotal")
```


### Donor cell number (GFP+/DAPI+) - total - Average per section
```{r}
# Plot Donor cell number - Total (i.e. across all regions - average per section
ggplot(MouseFFTotal, aes(x = GenotypeOrder, y = PS_GFP_DAPI, fill = GenotypeAgeOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  scale_y_continuous("Average number of donor cells (GFP+/DAPI+) per section", limits = c(0,310)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)

ggsave("1_Donor Cell Number_All.png", width = 30, height = 20, units = "cm", dpi = 300)
```


Run linear regression - no need for mixed effect model as only one measurement per mouse
Age_Inj not included in model as groups were age matched
```{r}
LM_Donor <- lm(PS_GFP_DAPI ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_Donor)

resid_panel(LM_Donor) # Residual plots don't look great. QQplot S-shaped 

shapiro.test(resid(LM_Donor)) # p-value = 0.001145

# Tukey transform the data to see if helps with assumptions (uses rcompanion)
MouseFFTotal$PS_GFP_DAPI_tukey = 
  transformTukey(MouseFFTotal$PS_GFP_DAPI,
                 plotit=FALSE)

LM_Donor_tukey <- lm(PS_GFP_DAPI_tukey ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_Donor_tukey)

resid_panel(LM_Donor_tukey) # plots look better

shapiro.test(resid(LM_Donor_tukey)) # p-value = 0.108

# Significant effect of Age_Inj. But no other significant results

# Box-Cox transformation - log likelihood method
lam.l = BoxCox.lambda(MouseFFTotal$PS_GFP_DAPI, method = "loglik") 
MouseFFTotal$PS_GFP_DAPI_BCL = BoxCox(MouseFFTotal$PS_GFP_DAPI, lam.l)

LM_Donor_BCL <- lm(PS_GFP_DAPI_BCL ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_Donor_BCL)

resid_panel(LM_Donor_BCL) # plots look similar to tukey

shapiro.test(resid(LM_Donor_BCL)) # p-value = 0.08803


# Box-Cox transformation - guerrero method
lam.g = BoxCox.lambda(MouseFFTotal$PS_GFP_DAPI, method = "guerrero") 
MouseFFTotal$PS_GFP_DAPI_BCG = BoxCox(MouseFFTotal$PS_GFP_DAPI, lam.g)

LM_Donor_BCG <- lm(PS_GFP_DAPI_BCG ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_Donor_BCG)

resid_panel(LM_Donor_BCG) # plots look good - probably the best

shapiro.test(resid(LM_Donor_BCG)) # p-value = 0.1169

anova(LM_Donor_BCG)
```


### Donor cell (GFP+/DAPI+) number - by region (ECSuper, ECDeep, DG, CA23, CA1) - average per section
```{r}
# Plot Donor cell number - Total (i.e. across all regions - average per section
ggplot(MouseFFRegions, aes(x = GenotypeOrder, y = PS_GFP_DAPI, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.75), alpha = 0.75, size = 1.5) + 
  #ggtitle("Average number of donor cells (GFP+/DAPI+) per section") + # title of the plot
  scale_y_continuous("Average number of donor cells (GFP+/DAPI+) per section", limits = c(0,205)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)

ggsave("1_Donor Cell Number_Region.png", width = 30, height = 20, units = "cm", dpi = 300)

```



Run mixed effects model - mixed effects used because multiple measurements per mouse
```{r}
# Mixed effect model, region multiplication
ME_Donor <- lmer(PS_GFP_DAPI ~ Genotype * AgeGroup * Brain_Region + Sex + (1|UnblindedID), data = MouseFFRegions)
summary(ME_Donor)

resid_panel(ME_Donor) # Residual plots look crap

shapiro.test(resid(ME_Donor)) # p-value < 2.2e-16

# Tukey transform - didnt work as have a negative value. Arefact of producing ROIs. Replace with 0.
MouseFFRegions$PS_GFP_DAPI_tukey = 
  transformTukey(MouseFFRegions$PS_GFP_DAPI,
                 plotit=FALSE)

ME_Donor_tukey <- lmer(PS_GFP_DAPI_tukey ~ Genotype * AgeGroup * Brain_Region + Sex + (1|UnblindedID), data = MouseFFRegions)
summary(ME_Donor_tukey)

resid_panel(ME_Donor_tukey) # Residual plots look better - best one to use

shapiro.test(resid(ME_Donor_tukey)) # p-value = 0.2622

anova(ME_Donor_tukey)
# Signigficant effect of Brain_Region
# Significant interactions - Genotype:Brain_Region, AgeGroup:Brain_Region (almost significant)

# Pairwise comparisons - won't do this
# emmeans(ME_Donor_tukey, pairwise ~ Brain_Region)

# Simple comparisons - Genotype:Brain_Region
emm <- emmeans(ME_Donor_tukey, ~ Genotype * Brain_Region)
pairs(emm, simple = "each") 

# Simple comparisons - AgeGroup:Brain_Region 
emm <- emmeans(ME_Donor_tukey, ~ AgeGroup * Brain_Region)
pairs(emm, simple = "each") 

# Box-Cox transformation - log likelihood method - Error in bcloglik(x, lower, upper) : x must be positive - can't use as have multiple values of 0


```

Plot Genotype:Brain_Region and AgeGroup:Brain_Region
```{r}
# Plot Donor cell number - Total (i.e. across all regions - average per section - Genotype:Brain_Region
ggplot(MouseFFRegions, aes(x = GenotypeOrder, y = PS_GFP_DAPI, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.75), alpha = 0.75, size = 1.5) + 
  #ggtitle("Average number of donor cells (GFP+/DAPI+) per section") + # title of the plot
  scale_y_continuous("Average number of donor cells (GFP+/DAPI+) per section", limits = c(-1,205)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"))
#  facet_wrap(~AgeGroup)

ggsave("1_Donor Cell Number_Region_Genotype*Region.png", width = 30, height = 20, units = "cm", dpi = 300)



# Plot Donor cell number - Total (i.e. across all regions - average per section - AgeGroup:Brain_Region
ggplot(MouseFFRegions, aes(x = AgeGroup, y = PS_GFP_DAPI, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.75), alpha = 0.75, size = 1.5) + 
  #ggtitle("Average number of donor cells (GFP+/DAPI+) per section") + # title of the plot
  scale_y_continuous("Average number of donor cells (GFP+/DAPI+) per section", limits = c(-1,205)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Age Group")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"))
#  facet_wrap(~AgeGroup)

ggsave("1_Donor Cell Number_Region_Age*Region.png", width = 30, height = 20, units = "cm", dpi = 300)
```


### Donor cell number (GFP+/DAPI+/Tau13+) - total - Average per section
```{r}
# Plot Donor cell number - Total (i.e. across all regions - average per section
ggplot(MouseFFTotal, aes(x = GenotypeOrder, y = PS_GFP_Tau13_DAPI, fill = GenotypeAgeOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  scale_y_continuous("Average number of donor cells 
  (GFP+/DAPI+/Tau13+) per section", limits = c(0,250)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)

ggsave("2_Donor Cell Number+Tau_Total.png", width = 30, height = 20, units = "cm", dpi = 300)
```


Run linear regression - no need for mixed effect model as only one measurement per mouse
Age_Inj not included in model as groups are age matched (could alternatively include Age_Cull)
```{r}
LM_Donor2 <- lm(PS_GFP_Tau13_DAPI ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_Donor2)

resid_panel(LM_Donor2) # Residual plots don't look great. QQplot bad. 

shapiro.test(resid(LM_Donor2)) # p-value = 7.139e-06-06

# Tukey transform the data to see if helps with assumptions (uses rcompanion)
MouseFFTotal$PS_GFP_Tau13_DAPI_tukey = 
  transformTukey(MouseFFTotal$PS_GFP_Tau13_DAPI,
                 plotit=FALSE)

LM_Donor2_tukey <- lm(PS_GFP_Tau13_DAPI_tukey ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_Donor2_tukey)

resid_panel(LM_Donor2_tukey) # plots look a lot better

shapiro.test(resid(LM_Donor2_tukey)) # p-value = 0.08789

anova(LM_Donor2_tukey)



# Box-Cox transformation - guerrero method
lam.g = BoxCox.lambda(MouseFFTotal$PS_GFP_Tau13_DAPI, method = "guerrero") 
MouseFFTotal$PS_GFP_Tau13_DAPI_BCG = BoxCox(MouseFFTotal$PS_GFP_Tau13_DAPI, lam.g)

LM_Donor2_BCG <- lm(PS_GFP_Tau13_DAPI_BCG ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_Donor2_BCG)

resid_panel(LM_Donor2_BCG) # plots look baaaad

shapiro.test(resid(LM_Donor2_BCG)) # p-value < 2.2e-16

# No significant results
```


### Donor cell (GFP+/DAPI+/Tau13+) number - by region (ECSuper, ECDeep, DG, CA23, CA1) - average per section
```{r}
# Plot Donor cell number - Total (i.e. across all regions - average per section
ggplot(MouseFFRegions, aes(x = GenotypeOrder, y = PS_GFP_Tau13_DAPI, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.75), alpha = 0.75, size = 1.5) + 
  #ggtitle("Average number of donor cells (GFP+/DAPI+/Tau13+) per section") + # title of the plot
  scale_y_continuous("Average number of donor cells (GFP+/DAPI+) per section", limits = c(0,150)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)

ggsave("2_Donor Cell Number+Tau_Regions.png", width = 30, height = 20, units = "cm", dpi = 300)
```

Run mixed effects model - mixed effects used because multiple measurements per mouse
```{r}
# Mixed effect model, region multiplication
ME_Donor2 <- lmer(PS_GFP_Tau13_DAPI ~ Genotype * AgeGroup * Brain_Region + Sex + (1|UnblindedID), data = MouseFFRegions)
summary(ME_Donor)

resid_panel(ME_Donor2) # Residual plots look crap

shapiro.test(resid(ME_Donor2)) # p-value < 2.2e-16

# Tukey transform - didnt work as have a negative value. Arefact of producing ROIs. Replace with 0.
MouseFFRegions$PS_GFP_Tau13_DAPI_tukey = 
  transformTukey(MouseFFRegions$PS_GFP_Tau13_DAPI,
                 plotit=FALSE)

ME_Donor2_tukey <- lmer(PS_GFP_Tau13_DAPI_tukey ~ Genotype * AgeGroup * Brain_Region + Sex + (1|UnblindedID), data = MouseFFRegions)
summary(ME_Donor2_tukey)

resid_panel(ME_Donor2_tukey) # Residual plots look pretty good - best one to use

shapiro.test(resid(ME_Donor2_tukey)) # p-value = 0.3876

anova(ME_Donor2_tukey)
# Signigficant effect of ABrain_Region, Age_Inj
# Significant interactions - Genotype:Brain_Region

# Pairwise comparisons - won't do this
emmeans(ME_Donor2_tukey, pairwise ~ Brain_Region)


# Simple comparisons - Genotype:Brain_Region
emm <- emmeans(ME_Donor2_tukey, ~ Genotype * Brain_Region)
pairs(emm, simple = "each") 


# Box-Cox transformation - log likelihood method - Error in bcloglik(x, lower, upper) : x must be positive - can't use as have multiple values of 0

# Box-Cox transformation - guerrero method
lam.g = BoxCox.lambda(MouseFFRegions$PS_GFP_Tau13_DAPI, method = "guerrero") 
MouseFFRegions$PS_GFP_Tau13_DAPI_BCG = BoxCox(MouseFFRegions$PS_GFP_Tau13_DAPI, lam.g)

ME_Donor2_BCG <- lmer(PS_GFP_Tau13_DAPI_BCG ~ Genotype * AgeGroup * Brain_Region + Sex + (1|UnblindedID), data = MouseFFRegions)
summary(ME_Donor2_BCG)

resid_panel(ME_Donor2_BCG) # plots not great

shapiro.test(resid(ME_Donor2_BCG)) # p-value = 2.251e-10


# Signigficant effect of Brain_Region
# Significant interactions - Genotype:Brain_Region

```

Plot Genotype:Brain_Region 
```{r}

ggplot(MouseFFRegions, aes(x = Brain_RegionOrder, y = PS_GFP_Tau13_DAPI, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.75), alpha = 0.75, size = 1.5) + 
  #ggtitle("Average number of donor cells (GFP+/DAPI+/Tau13+) per section") + # title of the plot
  scale_y_continuous("Average number of donor cells (GFP+/DAPI+) per section", limits = c(0,150)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"))

ggsave("2_Donor Cell Number+Tau_Regions-Region.png", width = 30, height = 20, units = "cm", dpi = 300)


# Plot Donor cell number - Regions - Genotype:Brain_Region
ggplot(MouseFFRegions, aes(x = GenotypeOrder, y = PS_GFP_Tau13_DAPI, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.75), alpha = 0.75, size = 1.5) + 
  #ggtitle("Average number of donor cells (GFP+/DAPI+/Tau13+) per section") + # title of the plot
  scale_y_continuous("Average number of donor cells (GFP+/DAPI+/Tau13+) per section", limits = c(0,150)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"))
#  facet_wrap(~AgeGroup)

ggsave("2_Donor Cell Number+Tau_Regions-Genotype*Region.png", width = 30, height = 20, units = "cm", dpi = 300)
```


### % of donor cells (GFP+/DAPI+) colocalising with Tau13 - total 
```{r}
# Plot Donor cells colocalising with Tau13 - Total (i.e. across all regions)
ggplot(MouseFFTotal, aes(x = GenotypeOrder, y = PercentGFP_DAPI_ColocTau13, fill = GenotypeOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Percent of donor cells (GFP+/DAPI+) colocalising with Tau13 (%)") + # title of the plot
  scale_y_continuous("Percent of donor cells (GFP+/DAPI+) colocalising with Tau13 (%)", limits = c(0,100)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)

ggsave("3_Donor Cell Coloc wTau_Total.png", width = 30, height = 20, units = "cm", dpi = 300)
```


Run linear regression - no need for mixed effect model as only one measurement per mouse
Age_Inj included in model as Adult group not age matched (could alternatively include Age_Cull)
```{r}
LM_Coloc <- lm(PercentGFP_DAPI_ColocTau13 ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_Coloc)

resid_panel(LM_Coloc) # Residual plots look a bit off

shapiro.test(resid(LM_Coloc)) # p-value = 0.0007974


# Tukey transform the data to see if helps with assumptions (uses rcompanion)
MouseFFTotal$PercentGFP_DAPI_ColocTau13_tukey = 
  transformTukey(MouseFFTotal$PercentGFP_DAPI_ColocTau13,
                 plotit=FALSE)

LM_Coloc_tukey <- lm(PercentGFP_DAPI_ColocTau13_tukey ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_Coloc_tukey)

resid_panel(LM_Coloc_tukey) # plots look really good

shapiro.test(resid(LM_Coloc_tukey)) # p-value = 0.7994

anova(LM_Coloc_tukey)

# Significant effect of AgeGroup - Aged have more

```

Plot by AgeGroup - Aged mice have higher percent
```{r}
ggplot(MouseFFTotal, aes(x = AgeGroup, y = PercentGFP_DAPI_ColocTau13, fill = AgeGroup)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Percent of donor cells (GFP+/DAPI+) colocalising with Tau13 (%)") + # title of the plot
  scale_y_continuous("Percent of donor cells (GFP+/DAPI+) colocalising with Tau13 (%)", limits = c(0,100)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Age Group")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"))
#  facet_wrap(~AgeGroup)

ggsave("3_Donor Cell Coloc wTau_Total_Age.png", width = 30, height = 20, units = "cm", dpi = 300)
```



### ### % of donor cells (GFP+/DAPI+) colocalising with Tau13 - by region (ECSuper, ECDeep, DG, CA23, CA1)
```{r}
# Plot Donor cell number - Total (i.e. across all regions - average per section
ggplot(MouseFFRegions, aes(x = GenotypeOrder, y = PercentGFP_DAPI_ColocTau13, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.75), alpha = 0.75, size = 1.5) + 
  #ggtitle("Percent of donor cells (GFP+/DAPI+) colocalising with Tau13 (%)") + # title of the plot
  scale_y_continuous("Percent of donor cells (GFP+/DAPI+) colocalising with Tau13 (%)", limits = c(0,120)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)

ggsave("3_Donor Cell Coloc wTau_Regions.png", width = 30, height = 20, units = "cm", dpi = 300)
```

Run mixed effects model - mixed effects used because multiple measurements per mouse
```{r}
# Mixed effect model, region multiplication
ME_Coloc2 <- lmer(PercentGFP_DAPI_ColocTau13 ~ Genotype * AgeGroup * Brain_Region + Sex + (1|UnblindedID), data = MouseFFRegions)
summary(ME_Coloc2)

resid_panel(ME_Coloc2) # Residual plots look really good

shapiro.test(resid(ME_Coloc2)) # p-value = 0.9068

anova(ME_Coloc2)
# Signigficant effect of Brain_Region
# Significant interactions - Genotype:Brain_Region, AgeGroup:Brain_Region

# Pairwise comparisons
emmeans(ME_Coloc2, pairwise ~ Brain_Region)

# Significant effect of Age_Inj. Interaction almost significant. But no other significant results
```

```{r}
# Plot by brian region
ggplot(MouseFFRegions, aes(x = Brain_RegionOrder, y = PercentGFP_DAPI_ColocTau13, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.75), alpha = 0.75, size = 1.5) + 
  #ggtitle("Percent of donor cells (GFP+/DAPI+) colocalising with Tau13 (%)") + # title of the plot
  scale_y_continuous("Percent of donor cells (GFP+/DAPI+) colocalising with Tau13 (%)", limits = c(0,120)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"))

ggsave("3_Donor Cell Coloc wTau_Regions_Region.png", width = 30, height = 20, units = "cm", dpi = 300)
```


### Recipient cell number (GFP-/Tau13+) - total - Average per section
```{r}
# Plot recipient cell number - Total (i.e. across all regions - average per section
ggplot(MouseFFTotal, aes(x = GenotypeOrder, y = PS_Recipient, fill = GenotypeAgeOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  scale_y_continuous("Average number of recipient cells 
(GFP-/Tau13+) per section", limits = c(0,2.6)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)

ggsave("1_Recipient cell_Total.png", width = 30, height = 20, units = "cm", dpi = 300)
```


Run linear regression - no need for mixed effect model as only one measurement per mouse
Age_Inj not included in model as groups age matched (could alternatively include Age_Cull)
```{r}
LM_recipient <- lm(PS_Recipient ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_recipient)

resid_panel(LM_recipient) # Residual plots don't look great. QQplot bendy wendy 

shapiro.test(resid(LM_recipient)) # p-value = 7.56e-05


# Tukey transform the data to see if helps with assumptions (uses rcompanion)
MouseFFTotal$PS_Recipient_tukey = 
  transformTukey(MouseFFTotal$PS_Recipient,
                 plotit=FALSE)

LM_recipient_tukey <- lm(PS_Recipient_tukey ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_recipient_tukey)

resid_panel(LM_recipient_tukey) # plots look better

shapiro.test(resid(LM_recipient_tukey)) # p-value = 0.205

anova(LM_recipient_tukey)

# No significant effects.

# Box-Cox transformation - log likelihood method. Error in bcloglik(x, lower, upper) : x must be positive. can't run as have 0's in dataset



# Box-Cox transformation - guerrero method
lam.g = BoxCox.lambda(MouseFFTotal$PS_Recipient, method = "guerrero") 
MouseFFTotal$PS_Recipient_BCG = BoxCox(MouseFFTotal$PS_Recipient, lam.g)

LM_recipient_BCG <- lm(PS_Recipient_BCG ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_recipient_BCG)

resid_panel(LM_recipient_BCG) # plots don't look great

shapiro.test(resid(LM_recipient_BCG)) # p-value = 7.173e-07



```


### Recipient cell number (GFP-/Tau13+)  - by region (ECSuper, ECDeep, DG, CA23, CA1) - average per section
```{r}
ggplot(MouseFFRegions, aes(x = GenotypeOrder, y = PS_Recipient, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.75), alpha = 0.75, size = 1.5) + 
# #ggtitle("Average number of recipient cells (GFP-/Tau13+) per section") + # title of the plot
  scale_y_continuous("Average number of recipient cells (GFP-/Tau13+) per section", limits = c(0,2.2)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)

ggsave("1_Recipient cell_Regions.png", width = 30, height = 20, units = "cm", dpi = 300)
```

Run mixed effects model - mixed effects used because multiple measurements per mouse
```{r}
# Mixed effect model, region multiplication
ME_Recipient2 <- lmer(PS_Recipient ~ Genotype * AgeGroup * Brain_Region + Sex + (1|UnblindedID), data = MouseFFRegions)
summary(ME_Recipient2)

resid_panel(ME_Recipient2) # Residual plots look crap

shapiro.test(resid(ME_Recipient2)) # p-value < 2.2e-16


# Tukey transform - didnt work as have a negative value. Arefact of producing ROIs. Replace with 0.
MouseFFRegions$PS_Recipient_tukey = 
  transformTukey(MouseFFRegions$PS_Recipient,
                 plotit=FALSE)

ME_Recipient2_tukey <- lmer(PS_Recipient_tukey ~ Genotype * AgeGroup * Brain_Region + Sex + (1|UnblindedID), data = MouseFFRegions)
summary(ME_Recipient2_tukey)

resid_panel(ME_Recipient2_tukey) # Residual plots look pretty bad - QQplot is bendy woooo

shapiro.test(resid(ME_Recipient2_tukey)) # p-value = 6.491e-12


# Box-Cox transformation - log likelihood method - Error in bcloglik(x, lower, upper) : x must be positive - can't use as have multiple values of 0


# Box-Cox transformation - guerrero method
lam.g = BoxCox.lambda(MouseFFRegions$PS_Recipient, method = "guerrero") 
MouseFFRegions$PS_Recipient_BCG = BoxCox(MouseFFRegions$PS_Recipient, lam.g)

ME_Recipient2_BCG <- lmer(PS_Recipient_BCG ~ Genotype * AgeGroup * Brain_Region + Sex + (1|UnblindedID), data = MouseFFRegions)
summary(ME_Recipient2_BCG)

resid_panel(ME_Recipient2_BCG) # plots not a bit better but still not wonderful

shapiro.test(resid(ME_Recipient2_BCG)) # p-value = 1.159e-09


anova(ME_Recipient2_BCG)
# Signigficant effect Brain_Region 
# Significant interactions - AgeGroup:Brain_Region

# Pairwise comparisons
emmeans(ME_Recipient2_BCG, pairwise ~ Brain_Region)
# Significant differences between ECSuper - DG, ECSuper - CA23, ECSuper - CA1, ECSuper - ECDeep, DG - CA23, DG - CA1,  

# Simple comparisons - :Brain_Region
emm <- emmeans(ME_Recipient2_BCG, ~ AgeGroup * Brain_Region)
pairs(emm, simple = "each") 
# Aged mice have more recipient cells in ECSuper, DG, CA23, CA1, ECDeep

# Significant effect of Age_Inj. Interaction almost significant. But no other significant results
```


Plot Brain_Region 
```{r}
# Plot recipient cell number - Regions - Brain_Region
ggplot(MouseFFRegions, aes(x = Brain_Region, y = PS_Recipient, fill = Brain_Region)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Average number of recipient cells (GFP-/Tau13+) per section") + # title of the plot
  scale_y_continuous("Average number of recipient cells (GFP-/Tau13+) per section", limits = c(0,2.2)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Brain Region")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"))
#  facet_wrap(~AgeGroup)

ggsave("1_Recipient cell_Regions_BrainRegion.png", width = 30, height = 20, units = "cm", dpi = 300)
```

Plot AgeGroup:Brain_Region 
```{r}
ggplot(MouseFFRegions, aes(x = AgeGroup, y = PS_Recipient, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.75), alpha = 0.75, size = 1.5) + 
  #ggtitle("Average number of recipient cells (GFP-/Tau13+) per section") + # title of the plot
  scale_y_continuous("Average number of recipient cells (GFP-/Tau13+) per section", limits = c(0,2.2)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76"))
#  facet_wrap(~AgeGroup)

ggsave("1_Recipient cell_Regions_Age*BrainRegion.png", width = 30, height = 20, units = "cm", dpi = 300)
```


### Recipient cell number (GFP-/Tau13+) normalised to donor cells (GFP+/DAPI+) - total - expressed as %
```{r}
ggplot(MouseFFTotal, aes(x = GenotypeOrder, y = RecipientNormGFP_DAPI_100, fill = GenotypeAgeOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+) (%)") + # title of the plot
  scale_y_continuous("Number of recipient cells (GFP-/Tau13+) 
  normalised to number of donor cells (GFP+/DAPI+) (%)", limits = c(0,7.1)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)

ggsave("2_Recipient cell norm GFPDAPI_Total.png", width = 30, height = 20, units = "cm", dpi = 300)
```


Run linear regression - no need for mixed effect model as only one measurement per mouse
Age_Inj included in model as Adult group not age matched (could alternatively include Age_Cull)
```{r}
LM_RecipientNorm <- lm(RecipientNormGFP_DAPI_100 ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_RecipientNorm)

resid_panel(LM_RecipientNorm) # Residual plots don't look great. QQplot bendy wendy 

shapiro.test(resid(LM_RecipientNorm)) # p-value = 1.872e-07

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
# Tukey transform the data to see if helps with assumptions (uses rcompanion)
MouseFFTotal$RecipientNormGFP_DAPI_100_tukey = 
  transformTukey(MouseFFTotal$RecipientNormGFP_DAPI_100,
                 plotit=FALSE)

LM_RecipientNorm_tukey <- lm(RecipientNormGFP_DAPI_100_tukey ~ Genotype * AgeGroup + Sex, data = MouseFFTotal)
summary(LM_RecipientNorm_tukey)

resid_panel(LM_RecipientNorm_tukey) # plots look good

shapiro.test(resid(LM_RecipientNorm_tukey)) # p-value = 0.5725

anova(LM_RecipientNorm_tukey)


# No significant effects
#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
# Box-Cox transformation - log likelihood method. Error in bcloglik(x, lower, upper) : x must be positive. can't run as have 0's in dataset

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---


#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
# No significant effects
```


### Recipient cell number (GFP-/Tau13+) normalised to donor cells (GFP+/DAPI+) - by region (ECSuper, ECDeep, DG, CA23, CA1) - expressed as % 
```{r}
# Plot Donor cell number - Total (i.e. across all regions - average per section
ggplot(MouseFFRegions, aes(x = GenotypeOrder, y = RecipientNormGFP_DAPI_100, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.75), alpha = 0.75, size = 1.5) + 
  #ggtitle("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+) (%)") + # title of the plot
  scale_y_continuous("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+) (%)") + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup) +
  facet_zoom(ylim = c(0,6.5), show.area = TRUE, shrink = TRUE)

ggsave("2_Recipient cell norm GFPDAPI_Regions.png", width = 30, height = 20, units = "cm", dpi = 300)
```


Run mixed effects model - mixed effects used because multiple measurements per mouse
```{r}
# Mixed effect model, region multiplication
ME_RecipientNorm <- lmer(RecipientNormGFP_DAPI_100 ~ Genotype * AgeGroup * Brain_Region + Sex + Age_Inj + (1|UnblindedID), data = MouseFFRegions)
summary(ME_RecipientNorm)

resid_panel(ME_RecipientNorm) # Residual plots look crap

shapiro.test(resid(ME_RecipientNorm)) # p-value < 2.2e-16

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
# Tukey transform - didnt work as have a negative value. Arefact of producing ROIs. Replace with 0.
MouseFFRegions$RecipientNormGFP_DAPI_100_tukey = 
  transformTukey(MouseFFRegions$RecipientNormGFP_DAPI_100,
                 plotit=FALSE)

ME_RecipientNorm_tukey <- lmer(RecipientNormGFP_DAPI_100_tukey ~ Genotype * AgeGroup * Brain_Region + Sex + (1|UnblindedID), data = MouseFFRegions)
summary(ME_RecipientNorm_tukey)

resid_panel(ME_RecipientNorm_tukey) # Residual plots look better but still not great

shapiro.test(resid(ME_RecipientNorm_tukey)) # p-value = 2.355e-15

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
# Box-Cox transformation - log likelihood method - Error in bcloglik(x, lower, upper) : x must be positive - can't use as have multiple values of 0


#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
# Box-Cox transformation - guerrero method
lam.g = BoxCox.lambda(MouseFFRegions$RecipientNormGFP_DAPI_100, method = "guerrero") 
MouseFFRegions$RecipientNormGFP_DAPI_100_BCG = BoxCox(MouseFFRegions$RecipientNormGFP_DAPI_100, lam.g)

ME_RecipientNorm_BCG <- lmer(RecipientNormGFP_DAPI_100_BCG ~ Genotype * AgeGroup * Brain_Region + Sex + (1|UnblindedID), data = MouseFFRegions)
summary(ME_RecipientNorm_BCG)

resid_panel(ME_RecipientNorm_BCG) # plots not a bit better but still not wonderful - best of the bunch

shapiro.test(resid(ME_RecipientNorm_BCG)) # p-value = 4.153e-10


anova(ME_RecipientNorm_BCG)
# Signigficant effect of AgeGroup, Brain_Region, Age_Inj 
# Aged mice have more recipient cells
# Significant interactions - AgeGroup:Brain_Region very close to being significant - 0.05072

# Pairwise comparisons
emmeans(ME_RecipientNorm_BCG, pairwise ~ Brain_Region)


# Simple comparisons - AgeGroup:Brain_Region
emm <- emmeans(ME_RecipientNorm_BCG, ~ AgeGroup * Brain_Region)
pairs(emm, simple = "each") 

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
```


Plot Brain_Region 
```{r}
# Plot recipient cell number - Regions - Brain_Region
ggplot(MouseFFRegions, aes(x = Brain_RegionOrder, y = RecipientNormGFP_DAPI_100, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+) (%)") + # title of the plot
  scale_y_continuous("Number of recipient cells (GFP-/Tau13+) 
  normalised to number of donor cells (GFP+/DAPI+) (%)") + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Brain Region")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_zoom(ylim = c(0,6.5), show.area = TRUE, shrink = TRUE)
#  facet_wrap(~AgeGroup)

ggsave("2_Recipient cell norm GFPDAPI_Regions_BrainRegion.png", width = 30, height = 20, units = "cm", dpi = 300)
```

Plot AgeGroup:Brain_Region 
```{r}
ggplot(MouseFFRegions, aes(x = AgeGroup, y = RecipientNormGFP_DAPI_100, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.8), alpha = 0.75, size = 1.5) + 
  #ggtitle("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+) (%)") + # title of the plot
  scale_y_continuous("Number of recipient cells (GFP-/Tau13+) 
  normalised to number of donor cells (GFP+/DAPI+) (%)") + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_zoom(ylim = c(0,6.5), show.area = TRUE, shrink = TRUE)
#  facet_wrap(~AgeGroup)

ggsave("2_Recipient cell norm GFPDAPI_Regions_Age*BrainRegion.png", width = 30, height = 20, units = "cm", dpi = 300)
```


### Recipient cell number (GFP-/Tau13+) normalised to donor cells (GFP+/DAPI+/Tau13+) - total - expressed as %
```{r}
# Plot recipient cell number - Total (i.e. across all regions)
ggplot(MouseFFTotal, aes(x = GenotypeOrder, y = RecipientNormGFP_Tau13_DAPI_100, fill = GenotypeAgeOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.2) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+/Tau13+) (%)") + # title of the plot
  scale_y_continuous("Number of recipient cells (GFP-/Tau13+) 
  normalised to number of donor cells (GFP+/DAPI+/Tau13+) (%)", limits = c(0,11.5)) + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="none") + #removes legend since labels on x axis sufficient
  guides(fill = FALSE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)

ggsave("2_Recipient cell norm to donor cells.png", width = 30, height = 20, units = "cm", dpi = 300)
```


Run linear regression - no need for mixed effect model as only one measurement per mouse
Age_Inj included in model as Adult group not age matched (could alternatively include Age_Cull)
```{r}
LM_RecipientNorm2 <- lm(RecipientNormGFP_Tau13_DAPI_100 ~ Genotype * AgeGroup + Sex + Age_Inj, data = MouseFFTotal)
summary(LM_RecipientNorm2)

resid_panel(LM_RecipientNorm2) # Residual plots don't look great. QQplot bendy wendy 

shapiro.test(resid(LM_RecipientNorm2)) # p-value = 9.126e-10

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
# Tukey transform the data to see if helps with assumptions (uses rcompanion)
MouseFFTotal$RecipientNormGFP_Tau13_DAPI_100_tukey = 
  transformTukey(MouseFFTotal$RecipientNormGFP_Tau13_DAPI_100,
                 plotit=FALSE)

LM_RecipientNorm2_tukey <- lm(RecipientNormGFP_Tau13_DAPI_100_tukey ~ Genotype * AgeGroup + Sex + Age_Inj, data = MouseFFTotal)
summary(LM_RecipientNorm2_tukey)

resid_panel(LM_RecipientNorm2_tukey) # plots look good

shapiro.test(resid(LM_RecipientNorm2_tukey)) # p-value = 0.3278

anova(LM_RecipientNorm2_tukey)

# Pairwise comparisons
#emmeans(LM_RecipientNorm2_tukey, pairwise ~ RecipientNormGFP_Tau13_DAPI_100_tukey ~ Genotype)

# No significant effects - Genotype p = 0.08254
#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
# Box-Cox transformation - log likelihood method. Error in bcloglik(x, lower, upper) : x must be positive. can't run as have 0's in dataset

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
# Box-Cox transformation - guerrero method
lam.g = BoxCox.lambda(MouseFFTotal$RecipientNormGFP_Tau13_DAPI_100, method = "guerrero") 
MouseFFTotal$RecipientNormGFP_Tau13_DAPI_100_BCG = BoxCox(MouseFFTotal$RecipientNormGFP_Tau13_DAPI_100, lam.g)

LM_RecipientNorm2_BCG <- lm(RecipientNormGFP_Tau13_DAPI_100_BCG ~ Genotype * AgeGroup + Sex + Age_Inj, data = MouseFFTotal)
summary(LM_RecipientNorm2_BCG)

resid_panel(LM_RecipientNorm2_BCG) # plots are baaaad

shapiro.test(resid(LM_RecipientNorm2_BCG)) # p-value = 2.219e-07

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
# No significant effects
```


### Recipient cell number (GFP-/Tau13+) normalised to donor cells (GFP+/DAPI+/Tau13+) - by region (ECSuper, ECDeep, DG, CA23, CA1) - expressed as % 
```{r}
# Plot Donor cell number - Total (i.e. across all regions - average per section
ggplot(MouseFFRegions, aes(x = GenotypeOrder, y = RecipientNormGFP_Tau13_DAPI_100, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.8), alpha = 0.75, size = 1.5) + 
  #ggtitle("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+/Tau13+) (%)") + # title of the plot
  scale_y_continuous("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+/Tau13+) (%)") + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup) +
  facet_zoom(ylim = c(0,11), show.area = TRUE, shrink = TRUE)
```


Run mixed effects model - mixed effects used because multiple measurements per mouse
```{r}
# Mixed effect model, region multiplication
ME_RecipientNorm2 <- lmer(RecipientNormGFP_Tau13_DAPI_100 ~ Genotype * AgeGroup * Brain_Region + Sex + Age_Inj + (1|UnblindedID), data = MouseFFRegions)
summary(ME_RecipientNorm2)

resid_panel(ME_RecipientNorm2) # Residual plots look crap

shapiro.test(resid(ME_RecipientNorm2)) # p-value < 2.2e-16

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
# Tukey transform - didnt work as have a negative value. Arefact of producing ROIs. Replace with 0.
MouseFFRegions$RecipientNormGFP_Tau13_DAPI_100_tukey = 
  transformTukey(MouseFFRegions$RecipientNormGFP_Tau13_DAPI_100,
                 plotit=FALSE)

ME_RecipientNorm2_tukey <- lmer(RecipientNormGFP_Tau13_DAPI_100_tukey ~ Genotype * AgeGroup * Brain_Region + Sex + Age_Inj + (1|UnblindedID), data = MouseFFRegions)
summary(ME_RecipientNorm2_tukey)

resid_panel(ME_RecipientNorm2_tukey) # Residual plots look better but still not great

shapiro.test(resid(ME_RecipientNorm2_tukey)) # p-value < 2.2e-16

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
# Box-Cox transformation - log likelihood method - Error in bcloglik(x, lower, upper) : x must be positive - can't use as have multiple values of 0


#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
# Box-Cox transformation - guerrero method
lam.g = BoxCox.lambda(MouseFFRegions$RecipientNormGFP_Tau13_DAPI_100, method = "guerrero") 
MouseFFRegions$RecipientNormGFP_Tau13_DAPI_100_BCG = BoxCox(MouseFFRegions$RecipientNormGFP_Tau13_DAPI_100, lam.g)

ME_RecipientNorm2_BCG <- lmer(RecipientNormGFP_Tau13_DAPI_100_BCG ~ Genotype * AgeGroup * Brain_Region + Sex + Age_Inj + (1|UnblindedID), data = MouseFFRegions)
summary(ME_RecipientNorm2_BCG)

resid_panel(ME_RecipientNorm2_BCG) # plots not a bit better but still not the best - best of the bunch

shapiro.test(resid(ME_RecipientNorm2_BCG)) # p-value = 1.2e-11


anova(ME_RecipientNorm2_BCG)
# Signigficant effect of Brain_Region. AgeGroup v close to significance
# Aged mice have more recipient cells

# Pairwise comparisons
#emmeans(ME_RecipientNorm2_BCG, pairwise ~ RecipientNormGFP_Tau13_DAPI_100_BCG ~ Brain_Region)
# Significant differences between ECSuper - DG, ECSuper - CA23, ECSuper - CA1, ECSuper - ECDeep, DG - CA23, DG - CA1, CA1 - ECDeep
#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
```

Plot AgeGroup 
```{r}
# Plot Donor cell number - Total Regions - Age_Group
ggplot(MouseFFRegions, aes(x = AgeGroup, y = RecipientNormGFP_Tau13_DAPI_100, fill = AgeGroup)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+/Tau13+) (%)") + # title of the plot
  scale_y_continuous("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+/Tau13+) (%)") + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("AgeGroup")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_zoom(ylim = c(0,11), show.area = TRUE, shrink = TRUE)
#  facet_wrap(~AgeGroup)
```

Plot Brain_Region 
```{r}
# Plot recipient cell number - Regions - Brain_Region
ggplot(MouseFFRegions, aes(x = Brain_Region, y = RecipientNormGFP_Tau13_DAPI_100, fill = Brain_Region)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex), position = position_jitter(w = 0.1, h = 0.0), alpha = 0.6, size = 1.5) + 
  #ggtitle("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+/Tau13+) (%)") + # title of the plot
  scale_y_continuous("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+/Tau13+) (%)") + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Brain Region")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_zoom(ylim = c(0,11), show.area = TRUE, shrink = TRUE)
#  facet_wrap(~AgeGroup)
```

Plot AgeGroup:Brain_Region 
```{r}
ggplot(MouseFFRegions, aes(x = AgeGroup, y = RecipientNormGFP_DAPI_100, fill = Brain_RegionOrder)) +
  stat_boxplot(geom ='errorbar', width = 0.75) +  #  adds the little horizontal lines on top and bottom of whiskers
  geom_boxplot(outlier.size = -1, width=0.75) +  # makes it not show outliers twice (because of jitter)
  #stat_summary(fun = max, fun.max = length, geom = "text", aes(label = ..ymax..), vjust = -2)  + # adds group n's above bars
  geom_point(aes(shape=Sex, group = Brain_RegionOrder), position = position_dodge(0.8), alpha = 0.75, size = 1.5) + 
  #ggtitle("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+) (%)") + # title of the plot
  scale_y_continuous("Number of recipient cells (GFP-/Tau13+) normalised to number of donor cells (GFP+/DAPI+) (%)") + # names the y axis, sets where ticks are and sets limits of y-axis
  xlab("Genotype")  + # sets the x-axis label
  theme_classic(base_size = 20) + #classic theme makes white background without lines. 
  theme(legend.position="right", axis.text.x = element_text(angle = 45, hjust = 1)) + #removes legend since labels on x axis sufficient
#  guides(fill = TRUE) + # removes group legend
#  scale_color_manual(values = c("cornflowerblue", “deeppink”)) + # manually set the colours you want - see word doc for link to list of colours
  #scale_fill_brewer(palette = "Purples", name = "GenotypeOrder") + 
  scale_fill_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_zoom(ylim = c(0,6.5), show.area = TRUE, shrink = TRUE)
#  facet_wrap(~AgeGroup)
```


## Correlation between donor cell (DAPI+/GFP+) and recipient cell counts (GFP-/Tau13+) per section - Total counts (all regions combined)

```{r}
# Scatter plot where data points are coloured according to genotype/Age group + regression line and CI - all data together
ggplot(MouseFFTotal, aes(x=PS_GFP_DAPI, y=PS_Recipient)) + 
  geom_point(aes(color = GenotypeAgeOrder, shape = Sex), size = 2) + 
  geom_smooth(method=lm, colour="darkgrey", se=TRUE) +
  xlab("Average number of donor cells (GFP+/DAPI+) per section") +
  ylab("Average number of recipient cells (GFP-/Tau13+) per section") +
  #ggtitle("Correlation between number of donor cells and recipient cells per section") +
  theme_classic(base_size = 20) +
  scale_color_manual(values = c("#FCC5C0", "#F768A1", "#AE017E", "#7A0177", "#99D8C9", "#66C2A4", "#238B45", "#00441B")) 
ggsave("4_Correlation_donor_recipient.png", width = 30, height = 10, units = "cm", dpi = 300)


ggplot(MouseFFTotal, aes(x=PS_GFP_DAPI, y=PS_Recipient)) + 
  geom_point(aes(color = GenotypeAgeOrder, shape = Sex), size = 2) + 
  geom_smooth(method=lm, colour="darkgrey", se=TRUE) +
  xlab("Average number of donor cells (GFP+/DAPI+) per section") +
  ylab("Average number of recipient cells 
  (GFP-/Tau13+) per section") +
  #ggtitle("Correlation between number of donor cells and recipient cells per section") +
  theme_classic(base_size = 20) +
  scale_color_manual(values = c("#FDE0DD","#FCC5C0", "#FA9FB5", "#F768A1", "#CCECE6", "#99D8C9", "#66C2A4", "#41AE76")) +
  facet_wrap(~AgeGroup)
ggsave("4_Correlation_donor_recipient_byage.png", width = 30, height = 15, units = "cm", dpi = 300)

# Scatter plot with separate plot and regression lines for each genotype
ggplot(MouseFFTotal, aes(x=PS_GFP_DAPI, y=PS_Recipient)) + 
  geom_point(aes(colour=GenotypeAgeOrder, shape = Sex), alpha = 1, size = 2) + 
  geom_smooth(method=lm, colour="darkgrey", se=TRUE, fullrange = TRUE) +
  xlab("Average number of donor cells (GFP+/DAPI+) per section") +
  ylab("Average number of recipient cells (GFP-/Tau13+) per section") +
  #ggtitle("Correlation between number of donor cells and recipient cells per section") +
  scale_color_manual(values = c("#FCC5C0", "#F768A1", "#AE017E", "#7A0177", "#99D8C9", "#66C2A4", "#238B45", "#00441B")) +
  theme(legend.position="none") +
  facet_wrap(~GenotypeAgeOrder, ncol = 2)
```

Calculate correlations
First need to see if variables are normally distributed
```{r}
shapiro.test(MouseFFTotal$PS_GFP_DAPI) # p-value = 1.517e-08
shapiro.test(MouseFFTotal$PS_Recipient) # p-value = 5.933e-08

# Neither are normally distributed
# need to use non-parametric correlation test - Spearman's correlation

cor.test(MouseFFTotal$PS_GFP_DAPI, MouseFFTotal$PS_Recipient,  method = "spearman")
# Significant positive correlation - rho = 0.6569039; p-value = 6.49e-13
```

Correlations in age/genotype groups
```{r}
KOAdultTot <- subset(MouseFFTotal, GenotypeAgeOrder == "APOEKO Adult")
E2AdultTot <- subset(MouseFFTotal, GenotypeAgeOrder == "APOE2 Adult")
E3AdultTot <- subset(MouseFFTotal, GenotypeAgeOrder == "APOE3 Adult")
E4AdultTot <- subset(MouseFFTotal, GenotypeAgeOrder == "APOE4 Adult")

KOAgedTot <- subset(MouseFFTotal, GenotypeAgeOrder == "APOEKO Aged")
E2AgedTot <- subset(MouseFFTotal, GenotypeAgeOrder == "APOE2 Aged")
E3AgedTot <- subset(MouseFFTotal, GenotypeAgeOrder == "APOE3 Aged")
E4AgedTot <- subset(MouseFFTotal, GenotypeAgeOrder == "APOE4 Aged")

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
### KO Adult ###
shapiro.test(KOAdultTot$PS_GFP_DAPI) # p-value = 0.0006804
shapiro.test(KOAdultTot$PS_Recipient) # p-value = 0.00188

# Both are not normally distributed - need to use non-parametric correlation test - Spearman's correlation

cor.test(KOAdultTot$PS_GFP_DAPI, KOAdultTot$PS_Recipient,  method = "spearman")
# Significant positive correlation- rho = 0.872081; p-value = 0.0001016

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
### E2 Adult ###
shapiro.test(E2AdultTot$PS_GFP_DAPI) # p-value = 0.0002219
shapiro.test(E2AdultTot$PS_Recipient) # p-value = 0.3139

# One  not normally distributed - need to use non-parametric correlation test - Spearman's correlation

cor.test(E2AdultTot$PS_GFP_DAPI, E2AdultTot$PS_Recipient,  method = "spearman")
# Significant positive correlation- rho = 0.5592308; p-value = 0.04692

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
### E3 Adult ###
shapiro.test(E3AdultTot$PS_GFP_DAPI) # p-value = 0.01157
shapiro.test(E3AdultTot$PS_Recipient) # p-value = 0.006493

# Both not normally distributed - need to use non-parametric correlation test - Spearman's correlation

cor.test(E3AdultTot$PS_GFP_DAPI, E3AdultTot$PS_Recipient,  method = "spearman")
# Significant positive correlation- rho = 0.6344057 ; p-value = 0.01481

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
### E4 Adult ###
shapiro.test(E4AdultTot$PS_GFP_DAPI) # p-value = 0.1403
shapiro.test(E4AdultTot$PS_Recipient) # p-value = 0.3467

# Both  normally distributed - need to use parametric correlation test - Pearson's correlation

cor.test(E4AdultTot$PS_GFP_DAPI, E4AdultTot$PS_Recipient,  method = "pearson")
# No significant positive - rho = 0.347988  ; p-value = 0.2677

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
### KO Aged ###
shapiro.test(KOAgedTot$PS_GFP_DAPI) # p-value = 0.9368
shapiro.test(KOAgedTot$PS_Recipient) # p-value = 0.04105

# One not normally distributed - need to use non-parametric correlation test - Spearman's correlation

cor.test(KOAgedTot$PS_GFP_DAPI, KOAgedTot$PS_Recipient,  method = "spearman")
# No significant correlation- rho = 0.6571429 ; p-value = 0.175

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
### E2 Aged ###
shapiro.test(E2AgedTot$PS_GFP_DAPI) # p-value = 0.2549
shapiro.test(E2AgedTot$PS_Recipient) # p-value = 0.05161

# Both normally distributed - need to use parametric correlation test - Pearson's correlation

cor.test(E2AgedTot$PS_GFP_DAPI, E2AgedTot$PS_Recipient,  method = "pearson")
# No significant correlation- rho = 0.4262717  ; p-value = 0.2193

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
### E3 Aged ###
shapiro.test(E3AgedTot$PS_GFP_DAPI) # p-value = 0.0008626
shapiro.test(E3AgedTot$PS_Recipient) # p-value = 0.008215

# Both not normally distributed - need to use non-parametric correlation test - Spearman's correlation

cor.test(E3AgedTot$PS_GFP_DAPI, E3AgedTot$PS_Recipient,  method = "spearman")
# Significant positive correlation- rho = 0.7264892 ; p-value = 0.003253

#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---#---
### E4 Aged ###
shapiro.test(E4AgedTot$PS_GFP_DAPI) # p-value = 0.01228
shapiro.test(E4AgedTot$PS_Recipient) # p-value = 0.03039

# Both not normally distributed - need to use non-parametric correlation test - Spearman's correlation

cor.test(E4AgedTot$PS_GFP_DAPI, E4AgedTot$PS_Recipient,  method = "spearman")
# No significant  correlation- rho = 0.4777316  ; p-value = 0.1162
```