epigenetic_age.Rmd

---
title: "R Notebook"
output: html_notebook
---

Library Load
```{r}
library(ggplot2)
library(magrittr)
library(dplyr)
# install.packages("broom")
# install.packages("ggbeeswarm")
library(ggbeeswarm)
library(broom)
# install.packages("ggpubr")
library(ggpubr)
# install.packages("patchwork")
library(patchwork)
library(data.table)
```

```{r}
hugo_epi_age = read.csv("/Users/hugohudson/Desktop/epigenetic_age_added.csv")
filtered_epi <- filter(hugo_epi_age,  hugo_epi_age$Epigenetic.Age > 0)
epi_age_lm = lm(ageDeath ~ filtered_epi$Epigenetic.Age, data = filtered_epi)

epigenetic_ages <- data.frame(Epigenetic.Age = filtered_epi$Epigenetic.Age)
predicted_ages <- predict(epi_age_lm, newdata = epigenetic_ages)
model_residuals <- -residuals(epi_age_lm)

ggplot(data = filtered_epi, mapping = aes(x = ageDeath, y = predicted_ages, color = primaryDiagnosis)) +
  geom_point() + 
  geom_abline(slope = 1, intercept = 0) + 
  xlim(0, 100) + ylim(0, 100)

ggplot(data = filtered_epi, mapping = aes(x = ageDeath, y = model_residuals, color = primaryDiagnosis)) + 
  geom_point() + 
  geom_hline(yintercept = 0) +
  geom_abline(slope = 1, intercept = 0)

ggplot(data = filtered_epi, mapping = aes(x = primaryDiagnosis, y = model_residuals)) + 
  geom_quasirandom() +
  stat_compare_means()

epi_age_vs_sst_proportions <- ggplot(data = filtered_epi, mapping = aes(color = primaryDiagnosis, x = predicted_ages, y = SST)) +
  geom_point() + 
  labs( title = "Epigenetic Age vs SST Proportions", 
        x = "Epigenetic Age", 
        y = "SST Proportion") +
  stat_smooth(method = "lm") + 
  stat_cor() 

chrono_age_vs_sst_proportions <- ggplot(data = filtered_epi, mapping = aes(color = primaryDiagnosis, x = ageDeath, y = SST)) +
  geom_point() + 
  labs( title = "Chronological Age vs SST Proportions", 
        x = "Chronological Age", 
        y = "SST Proportion") +
  stat_smooth(method = "lm") + 
  stat_cor() 

chrono_age_vs_sst_proportions
epi_age_vs_sst_proportions
# VIP cells interesting

ggplot(data = filtered_epi, mapping = aes(x = model_residuals, color = primaryDiagnosis)) + 
  geom_histogram() + 
  labs( title = "Residuals of the Model", 
        x = "Residuals", 
        y = "Count")

# Now want to graph the residuals when binned by age
filtered_epi$predicted_ages = predicted_ages
filtered_epi$chrono_age_groups <- cut(x = filtered_epi$ageDeath, breaks = seq(10, 100, by = 10), labels = c("10-19", "20-29", "30-39", "40-49", "50-59", "60-69", "70-79", "80-89", "90+"))

filtered_epi$epi_age_groups <- cut(x = filtered_epi$predicted_ages, breaks = seq(10, 100, by = 10), labels = c("10-19", "20-29", "30-39", "40-49", "50-59", "60-69", "70-79", "80-89", "90+"))

ggplot(data = filtered_epi, mapping = aes(x = chrono_age_groups, y = model_residuals, color = primaryDiagnosis)) +
  geom_boxplot() + 
  labs( x = "Chronological Age Groups", y = "Residuals", title = "Residuals for Chronological Age Groups")

ggplot(data = filtered_epi, mapping = aes(x = epi_age_groups, y = model_residuals, color = primaryDiagnosis)) +
  geom_boxplot() + 
  labs( x = "Epigenetic Age Groups", y = "Residuals", title = "Residuals for Epigenetic Age Groups")

epi_age_vs_VIP_proportions <- ggplot(data = filtered_epi, mapping = aes(color = primaryDiagnosis, x = predicted_ages, y = VIP)) +
  geom_point() + 
  labs( title = "Epigenetic Age vs VIP Proportions",  
        x = "Epigenetic Age", 
        y = "VIP Proportion") +
  stat_smooth(method = "lm") + 
  stat_cor() 

chrono_age_vs_VIP_proportions <- ggplot(data = filtered_epi, mapping = aes(color = primaryDiagnosis, x = ageDeath, y = VIP)) +
  geom_point() + 
  labs( title = "Chronological Age vs VIP Proportions", 
        x = "Chronological Age", 
        y = "VIP Proportion") +
  stat_smooth(method = "lm") + 
  stat_cor() 

epi_age_vs_VIP_proportions
chrono_age_vs_VIP_proportions

epi_age_vs_VIP_proportions <- ggplot(data = filtered_epi, mapping = aes(color = primaryDiagnosis, x = predicted_ages, y = Astrocyte)) +
  geom_point() + 
  labs( title = "Epigenetic Age vs Astrocyte Proportions",  
        x = "Epigenetic Age", 
        y = "Astrocyte Proportion") +
  stat_smooth(method = "lm") + 
  stat_cor() 

chrono_age_vs_VIP_proportions <- ggplot(data = filtered_epi, mapping = aes(color = primaryDiagnosis, x = ageDeath, y = Astrocyte)) +
  geom_point() + 
  labs( title = "Chronological Age vs Astrocyte Proportions", 
        x = "Chronological Age", 
        y = "Astrocyte Proportion") +
  stat_smooth(method = "lm") + 
  stat_cor() 

epi_age_vs_VIP_proportions
chrono_age_vs_VIP_proportions

epi_age_vs_VIP_proportions <- ggplot(data = filtered_epi, mapping = aes(color = primaryDiagnosis, x = predicted_ages, y = Pericyte)) +
  geom_point() + 
  labs( title = "Epigenetic Age vs Pericyte Proportions",  
        x = "Epigenetic Age", 
        y = "Pericyte Proportion") +
  stat_smooth(method = "lm") + 
  stat_cor() 

chrono_age_vs_VIP_proportions <- ggplot(data = filtered_epi, mapping = aes(color = primaryDiagnosis, x = ageDeath, y = Pericyte)) +
  geom_point() + 
  labs( title = "Chronological Age vs Pericyte Proportions", 
        x = "Chronological Age", 
        y = "Pericyte Proportion") +
  stat_smooth(method = "lm") + 
  stat_cor() 

epi_age_vs_VIP_proportions
chrono_age_vs_VIP_proportions

epi_age_vs_VIP_proportions <- ggplot(data = filtered_epi, mapping = aes(color = primaryDiagnosis, x = predicted_ages, y = LAMP5)) +
  geom_point() + 
  labs( title = "Epigenetic Age vs LAMP5 Proportions",  
        x = "Epigenetic Age", 
        y = "LAMP5 Proportion") +
  stat_smooth(method = "lm") + 
  stat_cor() 

chrono_age_vs_VIP_proportions <- ggplot(data = filtered_epi, mapping = aes(color = primaryDiagnosis, x = ageDeath, y = LAMP5)) +
  geom_point() + 
  labs( title = "Chronological Age vs LAMP5 Proportions", 
        x = "Chronological Age", 
        y = "LAMP5 Proportion") +
  stat_smooth(method = "lm") + 
  stat_cor() 

epi_age_vs_VIP_proportions
chrono_age_vs_VIP_proportions

filtered_epi$
```

Import EpiDISH / CellDMC Required Packages
```{r}
#if (!require("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")

#BiocManager::install("EpiDISH")
library(EpiDISH)
```

Loading example data / running workflow
```{r}
data(centEpiFibIC.m)
data(DummyBeta.m)
frac.m <- epidish(beta.m = DummyBeta.m, ref.m = centEpiFibIC.m, method = "RPC")
dim(frac.m$dataREF)
dim(frac.m$estF)
pheno.v <- rep(c(0, 1), each = 5)
dim(pheno.v)
celldmc.o <- CellDMC(DummyBeta.m, pheno.v, frac.m$estF)
head(celldmc.o$coe)
class(frac.m$estF)
frac.m$estF
```

```{r}

# Library Loading
library(ggplot2)
library(magrittr)
library(dplyr)
library(ggbeeswarm)
library(broom)
library(ggpubr)
library(patchwork)
library(data.table)
library(EpiDISH)
library(readr)
library(data.table)

# Load Cell Fractions Data
hugo_epi_age = read.csv("/Users/hugohudson/Desktop/Lab/epigenetic_age_added.csv")
filtered_epi <- filter(hugo_epi_age,  Epigenetic.Age > 0)
brnums_dan <- filtered_epi$individualID
pheno_final <- filtered_epi$primaryDiagnosis
# Re-scale the fractions
# Define Fraction Rescaling Function
softmax_row <- function(z) {
  exp_z <- exp(z)                 # Exponentiate each element in the row
  exp_z / sum(exp_z)              # Normalize to sum to 1
}
fracs <- filtered_epi[1:nrow(filtered_epi),38:ncol(filtered_epi)-1]
frac_df <- data.frame(fracs)


# Load the patient information
patient_info <- read.csv('/Users/hugohudson/Desktop/Lab/GSE74193_series_matrix.csv')
# Sort the info by fractions
brnums <- patient_info[43,2:ncol(patient_info)]
brnums <- substr(brnums, 50, nchar(brnums))
valid_positions <- na.omit(match(brnums, brnums_dan))
brnums_filtered <- brnums[valid_positions]
pheno <- patient_info[48, 2:ncol(patient_info)]
pheno <- substr(pheno, 5, nchar(pheno))
# Read the CpG file
file <- '/Users/hugohudson/Desktop/Lab/GSE74193_GEO_procData.csv'
col_names <- names(fread(file, nrows = 0))
beta_columns <- col_names[grepl("Beta", col_names)]
dataFiltered <- fread(file, select = beta_columns)
colnames(dataFiltered) <- brnums
#dataFiltered <- fread(file, select = seq(2, ncol(chunk), by = 2))
# Sort the CpG file
brnums_filtered <- brnums[brnums %in% brnums_dan]
brnums_sorted <- brnums_filtered[order(match(brnums_filtered, brnums_dan))]
pheno <- pheno[order(match(brnums_filtered, brnums_dan))]
#fracs_df <- frac_df[order(match(brnums_filtered, brnums_dan)), 1:dim(frac_df)[2]]
# Apply the softmax function to each row of the data frame
df_normalized <- apply(frac_df, 1, softmax_row)
# Convert the result back into a data frame
df_normalized <- as.matrix(df_normalized)
unique_elements <- unique(brnums_filtered)
filtered_df <- dataFiltered %>% select(any_of(unique_elements))
dim(df_normalized)
select(dataFiltered, any_of(unique_elements))
filtered_df <- as.matrix(filtered_df)
result <- CellDMC(beta.m = filtered_df, pheno.v = pheno_final, frac.m = t(df_normalized), 
                     adjPMethod = 'fdr', adjPThresh = 0.05)

# Write the coefficients to a csv
coefficient_dataframe <- data.frame()
cell_names <- colnames(df_normalized)
for (i in length(result$coe)) {
  coefficient_dataframe$cell_names[i] = result$coe[i]$Estimate
}
write.csv(coefficient_dataframe, '/Users/hugohudson/Desktop/Lab/CellCoefficients.csv')
coefficients_list <- celldmc.o$coe

# Combine coefficients into a single data frame
# Assuming coefficients_list is a list of data frames for each cell type
result$coe[1]
write.csv(result$coe[1], '/Users/hugohudson/Desktop/Lab/AstrocyteCoefficients.csv')
for (i in length(result$coe)) {
  write.csv(result$coe[i], paste0('/Users/hugohudson/Desktop/Lab/',cell_names[i],'Coefficients.csv'))
}
cell_is <- 'yellow'
paste0('hello',cell_is)
paste0('/Users/hugohudson/Desktop/Lab/',cell_names[1],'Coefficients.csv')
```

```{r}
the_ages <- read.csv('/Users/hugohudson/Desktop/Lab/AstrocyteResults.csv')
length(the_ages$X0)

hugo_epi_age = read.csv("/Users/hugohudson/Desktop/Lab/epigenetic_age_added.csv")
hugo_epi_age <- hugo_epi_age[1:1350, 1:ncol(hugo_epi_age)]
filtered_epi <- filter(hugo_epi_age,  hugo_epi_age$Epigenetic.Age > 0)
hugo_epi_age$CELLAGE = the_ages$X0
epi_age_lm = lm(ageDeath ~ CELLAGE, data = hugo_epi_age)
new_df$CELLAGE <- hugo_epi_age$CELLAGE
new_df
predicted_ages <- predict(epi_age_lm, newdata = new_df)
model_residuals <- -residuals(epi_age_lm)
ggplot(data = hugo_epi_age, mapping = aes(x = ageDeath, y = predicted_ages, color = primaryDiagnosis)) +
  geom_point() + 
  geom_abline(slope = 1, intercept = 0) + 
  xlim(0, 100) + ylim(0, 100)
```

Reworked Workflow
```{r}
result$dmct
# Library Loading
library(ggplot2)
library(magrittr)
library(dplyr)
library(ggbeeswarm)
library(broom)
library(ggpubr)
library(patchwork)
library(data.table)
library(EpiDISH)
library(readr)
library(data.table)

# Load Cell Fractions Data
hugo_epi_age = read.csv("/Users/hugohudson/Desktop/Lab/epigenetic_age_added.csv")
filtered_epi <- filter(hugo_epi_age,  Epigenetic.Age > 0)
brnums_dan <- filtered_epi$individualID
pheno <- filtered_epi$primaryDiagnosis
# Re-scale the fractions
# Define Fraction Rescaling Function
softmax_row <- function(z) {
  exp_z <- exp(z)                 # Exponentiate each element in the row
  exp_z / sum(exp_z)              # Normalize to sum to 1
}
fracs <- filtered_epi[1:nrow(filtered_epi),38:ncol(filtered_epi)-1]
frac_df <- data.frame(fracs)
# Apply the softmax function to each row of the data frame
df_normalized <- apply(frac_df, 1, softmax_row)
# Convert the result back into a data frame
df_normalized <- as.matrix(df_normalized)

# Load the patient information
patient_info <- read.csv('/Users/hugohudson/Desktop/Lab/GSE74193_series_matrix.csv')
# Sort the info by fractions
brnums <- patient_info[43,2:ncol(patient_info)]
brnums <- substr(brnums, 50, nchar(brnums))
valid_positions <- na.omit(match(brnums, brnums_dan))
brnums_filtered <- brnums[valid_positions]

# Load the CpG data
# Initialize an empty list to store chunks
chunk_list <- list()

# Set chunk size (number of rows to read per chunk)
chunk_size <- 10000  # You can adjust this based on memory availability

# Read the file in chunks
file <- '/Users/hugohudson/Desktop/Lab/GSE74193_GEO_procData.csv'
data_filtered <- fread(file, select = seq(2, ncol(chunk), by = 2))
head(data_filtered)
# Keep track of the total number of rows read
total_rows_read <- 0

# Read and process in chunks
repeat {
  # Read the next chunk
  chunk <- read_csv(file, skip = total_rows_read, n_max = chunk_size, col_types = cols())
  
  # Break the loop if the chunk is empty (end of file)
  if (nrow(chunk) == 0) break
  
  # Perform your operations on the chunk (e.g., filtering or transformations)
  # chunk <- chunk %>% filter(some_condition)
  # Filter for only the beta rows
  #beta_cols <- chunk[, grepl("Beta", names(chunk))]
  beta_cols <- chunk[, seq(1, ncol(chunk), by = 2)]
  betas_filtered <- beta_cols[, valid_positions]
  # Add the processed chunk to the list
  chunk_list[[length(chunk_list) + 1]] <- betas_filtered
  
  # Update the total number of rows read
  total_rows_read <- total_rows_read + nrow(beta_cols)
}
# Run CellDMC
results_list <- list()
for (i in seq(1, length(brnums_filtered), by = 5)) {
  # Gather Chunks Information - Use min to handle the last group if it has less than 5 elements
  # Need fractions
  chunk_fractions <- df_normalized[i:min(i + 4, length(filtered_brnums)), 1:ncol(df_normalized)]
  # Need CpGs
  #chunk_CpGs <- chunk_list[, 1:nrow(chunk_list[1]), i:min(i + 4, length(filtered_brnums))]
  #chunk_CpGs <- chunk_list[, 1:nrow(chunk_list[[1]]), drop = FALSE][i:min(i + 4, length(filtered_brnums))]
  chunk_CpGs <- lapply(chunk_list, function(chunk) {
  # Calculate the number of rows in the current chunk
  n_rows <- nrow(chunk)
  
  # Determine the end index for rows to extract
  end_index <- min(i + n_rows - 1, n_rows)
  
  # Extract the desired rows and all columns
  chunk[i:end_index, , drop = FALSE]
})
  # Need Group Info
  chunk_groups <- pheno[i:min(i + 4, length(filtered_brnums))]
  # Call CellDMC
  result <- CellDMC(beta.m = chunk_CpGs, pheno.v = chunk_groups, frac.m = chunk_fractions, 
                     adjPMethod = 'fdr', adjPThresh = 0.05)
  # Add the result to the list
  results_list[[length(results_list) + 1]] <- result
}
# Initialize variables to store the sum of column averages and the count of data frames
sum_averages <- NULL
count <- 0

# Step 1: Calculate the column averages incrementally
for (df in results_list) {
  # Calculate the column averages for the current data frame
  column_averages <- colMeans(df)
  
  # Initialize sum_averages if it's NULL
  if (is.null(sum_averages)) {
    sum_averages <- column_averages
  } else {
    sum_averages <- sum_averages + column_averages
  }
  
  # Increment the count of data frames processed
  count <- count + 1
}

# Step 2: Calculate the overall average of the column sums
overall_averages <- sum_averages / count

# Define the probe names


# Convert to DF
coefficient_dataframe <- as.data.frame(probe_names)
coefficient_dataframe$coefficients <- overall_averages

# Write the coefficients to a csv
write.csv(result$dmct, '/Users/hugohudson/Desktop/Lab/DMCT.csv')

class(result$dmct)

```



```{r}
library(readr)
library(dplyr)

# Define Fraction Rescaling Function
softmax_row <- function(z) {
  exp_z <- exp(z)                 # Exponentiate each element in the row
  exp_z / sum(exp_z)              # Normalize to sum to 1
}

# Initialize an empty list to store chunks
chunk_list <- list()

# Set chunk size (number of rows to read per chunk)
chunk_size <- 10000  # You can adjust this based on memory availability

# Read the file in chunks
file <- '/Users/hugohudson/Desktop/Lab/GSE74193_GEO_procData.csv'

# Keep track of the total number of rows read
total_rows_read <- 0

# Read and process in chunks
repeat {
  # Read the next chunk
  chunk <- read_csv(file, skip = total_rows_read, n_max = chunk_size, col_types = cols())
  
  # Break the loop if the chunk is empty (end of file)
  if (nrow(chunk) == 0) break
  
  # Perform your operations on the chunk (e.g., filtering or transformations)
  # chunk <- chunk %>% filter(some_condition)
  # Filter for only the beta rows
  #beta_cols <- chunk[, grepl("Beta", names(chunk))]
  beta_cols <- chunk[, seq(1, ncol(chunk), by = 2)]
  # Add the processed chunk to the list
  chunk_list[[length(chunk_list) + 1]] <- beta_cols
  
  # Update the total number of rows read
  total_rows_read <- total_rows_read + nrow(beta_cols)
}

# Combine the list of chunks into a single data frame
#final_data <- bind_rows(chunk_list)

# DATA SETUP / ORGANIZATION NOW

# Load the patient information
patient_info <- read.csv('/Users/hugohudson/Desktop/Lab/GSE74193_series_matrix.csv')
# Extract phenotype information
patient_info <- data.frame(patient_info)
# Load the cell fractions 
hugo_epi_age = read.csv("/Users/hugohudson/Desktop/Lab/epigenetic_age_added.csv")
filtered_epi <- filter(hugo_epi_age,  Epigenetic.Age > 0)
# Check that the brnum samples are in alignment
brnums <- patient_info[43,2:ncol(patient_info)]
brnums
brnums <- substr(brnums, 50, nchar(brnums))
brnums
brnums_dan <- filtered_epi$individualID
filtered_brnums <- brnums_dan[brnums_dan %in% brnums]
filtered_data <- filter(filtered_epi, individualID %in% brnums)
# Sort CpG Data
# Might be able to form into one df at this point
chunk_list
# final_data <- bind_rows(chunk_list)
setnames(final_data, new = brnums)
sorted_df <- final_data[, match(filtered_brnums, colnames(final_data))]

# Now collect other data required for cell DMC
# Now in same order as Dan's data, so can use filtered_epi
pheno <- filtered_epi$primaryDiagnosis
fracs <- filtered_epi[1:nrow(filtered_epi),38:ncol(filtered_epi)-1]
frac_df <- data.frame(fracs)
softmax_row <- function(z) {
  exp_z <- exp(z)                 # Exponentiate each element in the row
  exp_z / sum(exp_z)              # Normalize to sum to 1
}
# Apply the softmax function to each row of the data frame
df_normalized <- apply(frac_df, 1, softmax_row)
# Convert the result back into a data frame
df_normalized <- as.matrix(df_normalized)
sorted_df <- as.matrix(sorted_df)
# Now can run cell DMC
#sorted_df
#celldmc.o <- CellDMC(beta.m = sorted_df, pheno.v = pheno, frac.m = df_normalized, adjPMethod = 'fdr', adjPThresh = 0.05)
#head(celldmc.o$coe)

# OTHER APPROACH 
# Placeholder for storing results
#chunk_list_func <- function(x, chunk_size) {
#    split(x, ceiling(seq_along(x) / chunk_size))
#}

# Chunk the list into groups of 10,000
#chunk_df <- function(df, chunk_size) {
#    split(df, ceiling(seq_along(df) / chunk_size))
#}

# Chunk the data frame into groups of 10,000 rows
#chunked_df <- chunk_df(fracs, 10000)
#chunked_list$'0.018'
results_list <- list()
library(data.table)
sorted_df
# Loop through each chunk (assuming feature_chunks is your list of 40 chunks)
for (i in seq_along(chunk_list)) {
    # Extract the current chunk
    current_chunk <- as.data.frame(chunk_list[i])  # Each chunk is (10000 x 180)
    setnames(current_chunk, new = c('probes', brnums))
    print(current_chunk)
    valid_positions <- na.omit(match(filtered_brnums, colnames(current_chunk)))
    reordered_chunk <- current_chunk[, valid_positions]
    current_chunk <- as.data.frame(current_chunk, )
    # Fit the model using CellDMC
    result <- CellDMC(beta.m = current_chunk, pheno.v = pheno, frac.m = sorted_df, 
                     adjPMethod = 'fdr', adjPThresh = 0.05)
    
    # Store the coefficients or relevant metrics from the result
    results_list[i] <- result$coe  # Assuming you want the coefficients
}

# Combine results into a single matrix or data frame
combined_results <- do.call(rbind, results_list)

# Calculate the mean coefficients across chunks if needed
mean_coefficients <- rowMeans(combined_results)

# Inspect the mean coefficients
mean_coefficients
```

```{r}
# Load the probe data
#beta_matrix <- read.csv('/Users/hugohudson/Desktop/GSE74193_GEO_procData.csv')
# Load the patient information
patient_info <- read.csv('/Users/hugohudson/Desktop/GSE74193_series_matrix.csv')
# Extract phenotype information
patient_info <- data.frame(patient_info)
# Load the cell fractions 
hugo_epi_age = read.csv("/Users/hugohudson/Desktop/epigenetic_age_added.csv")
filtered_epi <- filter(hugo_epi_age,  Epigenetic.Age > 0)
dim(filtered_epi)
filtered_epi$individualID
# Check that the brnum samples are in alignment
brnums <- patient_info[43,2:ncol(patient_info)]
brnums
brnums <- substr(brnums, 50, nchar(brnums))
brnums
brnums_dan <- filtered_epi$individualID
filtered_brnums <- brnums_dan[brnums_dan %in% brnums]
filtered_data <- filter(filtered_epi, individualID %in% brnums)
# Create a dataframe with the brnums and cpg info
final_data <- final_data[1:nrow(final_data), 1:length(brnums)]
brnums <- unlist(brnums)
setnames(final_data, new = brnums)
sorted_df <- final_data[, match(filtered_brnums, colnames(final_data))]
# Now in same order as Dan's data, so can use filtered_epi
pheno <- filtered_epi$primaryDiagnosis
fracs <- filtered_epi[1:nrow(filtered_epi),38:ncol(filtered_epi)-1]
frac_df <- data.frame(fracs)
softmax_row <- function(z) {
  exp_z <- exp(z)                 # Exponentiate each element in the row
  exp_z / sum(exp_z)              # Normalize to sum to 1
}
# Apply the softmax function to each row of the data frame
df_normalized <- t(apply(frac_df, 1, softmax_row))
# Convert the result back into a data frame
df_normalized <- as.matrix(df_normalized)
sorted_df <- as.matrix(sorted_df)
```


#frac_df[] <- lapply(frac_df, as.numeric)
frac_matrix <- as.matrix(df_normalized)
sorted_df <- as.matrix(sorted_df)
```

Running CellDMC workflow on our data
```{r}
celldmc.o <- CellDMC(beta.m = sorted_df, pheno.v = pheno, frac.m = frac_matrix, adjPMethod = 'fdr', adjPThresh = 0.05)
head(celldmc.o$coe)
sum(beta_matrix$Br1764 * celldmc.o$coe$Astrocyte$Estimate)
```

Load horvath clock
```{r}
# Load the package / functions
#install.packages("WGCNA")
#install.packages("sqldf")
#install.packages("BiocManager")
#BiocManager::install("GO.db")1
# BiocManager::install("impute")
#install the Bioconductor installer
#install.packages("BiocInstaller",repos="http://www.bioconductor.org/packages/2.13/bioc")
#install “impute” from Bioconductor 
#biocLite("impute")
#BiocManager::install("preprocessCore")
library(WGCNA)
library(sqldf)
#BiocManager::install("BMIQ")
#library(BMIQ)
# Load the betas
#horvath_betas <- read.csv('/Users/hugohudson/Downloads/gb-2013-14-10-r115-S3.csv')
# Load the example data
library(RPMM)
```

Test horvath clock on horvath example data
```{r}
# Age transofrmation and probe annotation functions
trafo= function(x,adult.age=20) { x=(x+1)/(1+adult.age); y=ifelse(x<=1, log( x),x-1);y }
anti.trafo= function(x,adult.age=20) { ifelse(x<0, (1+adult.age)*exp(x)-1, (1+adult.age)*x+adult.age) }
probeAnnotation21kdatMethUsed=read.csv("/Users/hugohudson/Downloads/gb-2013-14-10-r115-S22.csv")
probeAnnotation27k=read.csv("/Users/hugohudson/Downloads/gb-2013-14-10-r115-S21.csv")
datClock=read.csv("/Users/hugohudson/Downloads/gb-2013-14-10-r115-S23.csv")

# Read in DNA methylation data
# For a small file, e.g. measured on the 27k platform you could just use read.csv. 
# But for large files, e.g. those measured on the 450K platform, I recommend you use read.csv.sql.
dat0=read.csv.sql("/Users/hugohudson/Downloads/gb-2013-14-10-r115-S26.csv") ;
nSamples=dim(dat0)[[2]]-1
nProbes= dim(dat0)[[1]]
# the following command may not be needed. But it is sometimes useful when you use read.csv.sql
dat0[,1]= gsub(x=dat0 [,1],pattern="\"",replacement="") 
#Create a log file which will be output into your directory
# The code looks a bit complicated because it serves to create a log file (for error checks etc).
# It will automatically create a log file.
file.remove("LogFile.txt")
file.create("LogFile.txt")
DoNotProceed=FALSE
cat(paste( "The methylation data set contains", nSamples, "samples (e.g. arrays) and ", nProbes, " probes."),file="LogFile.txt")
if (nSamples==0) {DoNotProceed=TRUE; cat(paste( "\n ERROR: There must be a data input error since there seem to be no samples.\n Make sure that you input a comma delimited file (.csv file)\n that can be read using the R command read.csv.sql . Samples correspond to columns in that file  ."), file="LogFile.txt",append=TRUE) } 
if (nProbes==0) {DoNotProceed=TRUE; cat(paste( "\n ERROR: There must be a data input error since there seem to be zero probes.\n Make sure that you input a comma delimited file (.csv file)\n that can be read using the R command read.csv.sql  CpGs correspond to rows.")   , file="LogFile.txt",append=TRUE) } 
if (  nSamples > nProbes  ) { cat(paste( "\n MAJOR WARNING: It worries me a lot that there are more samples than CpG probes.\n Make sure that probes correspond to rows and samples to columns.\n I wonder whether you want to first transpose the data and then resubmit them? In any event, I will proceed with the analysis."),file="LogFile.txt",append=TRUE) }
if (  is.numeric(dat0[,1]) ) { DoNotProceed=TRUE; cat(paste( "\n Error: The first column does not seem to contain probe identifiers (cg numbers from Illumina) since these entries are numeric values. Make sure that the first column of the file contains probe identifiers such as cg00000292. Instead it contains ", dat0[1:3,1]  ),file="LogFile.txt",append=TRUE)  } 
if (  !is.character(dat0[,1]) ) {  cat(paste( "\n Major Warning: The first column does not seem to contain probe identifiers (cg numbers from Illumina) since these entries are numeric values. Make sure that the first column of the file contains CpG probe identifiers such as cg00000292. Instead it contains ", dat0[1:3,1]  ),file="LogFile.txt",append=TRUE)  } 
datout=data.frame(Error=c("Input error. Please check the log file for details","Please read the instructions carefully."), Comment=c("", "email Steve Horvath."))
if ( ! DoNotProceed ) {
nonNumericColumn=rep(FALSE, dim(dat0)[[2]]-1)
for (i in 2:dim(dat0)[[2]] ){ nonNumericColumn[i-1]=! is.numeric(dat0[,i]) }
if (  sum(nonNumericColumn) >0 ) { cat(paste( "\n MAJOR WARNING: Possible input error. The following samples contain non-numeric beta values: ", colnames(dat0)[-1][ nonNumericColumn], "\n Hint: Maybe you use the wrong symbols for missing data. Make sure to code missing values as NA in the Excel file. To proceed, I will force the entries into numeric values but make sure this makes sense.\n" ),file="LogFile.txt",append=TRUE)  } 
XchromosomalCpGs=as.character(probeAnnotation27k$Name[probeAnnotation27k$Chr=="X"])
selectXchromosome=is.element(dat0[,1], XchromosomalCpGs )
selectXchromosome[is.na(selectXchromosome)]=FALSE
meanXchromosome=rep(NA, dim(dat0)[[2]]-1)
if (   sum(selectXchromosome) >=500 )  {
meanXchromosome= as.numeric(apply( as.matrix(dat0[selectXchromosome,-1]),2,mean,na.rm=TRUE)) }
if (  sum(is.na(meanXchromosome)) >0 ) { cat(paste( "\n \n Comment: There are lots of missing values for X chromosomal probes for some of the samples. This is not a problem when it comes to estimating age but I cannot predict the gender of these samples.\n " ),file="LogFile.txt",append=TRUE)  } 

match1=match(probeAnnotation21kdatMethUsed$Name , dat0[,1])
if  ( sum( is.na(match1))>0 ) { 
missingProbes= probeAnnotation21kdatMethUsed$Name[!is.element( probeAnnotation21kdatMethUsed$Name , dat0[,1])]    
DoNotProceed=TRUE; cat(paste( "\n \n Input error: You forgot to include the following ", length(missingProbes), " CpG probes (or probe names):\n ", paste( missingProbes, sep="",collapse=", ")),file="LogFile.txt",append=TRUE)  } 
#STEP 2: Restrict the data to 21k probes and ensure they are numeric
match1=match(probeAnnotation21kdatMethUsed$Name , dat0[,1])
if  ( sum( is.na(match1))>0 ) stop(paste(sum( is.na(match1)), "CpG probes cannot be matched"))
dat1= dat0[match1,]
asnumeric1=function(x) {as.numeric(as.character(x))}
dat1[,-1]=apply(as.matrix(dat1[,-1]),2,asnumeric1)

#STEP 3: Create the output file called datout
set.seed(1)
# Do you want to normalize the data (recommended)?
normalizeData=TRUE
#source("/Users/hugohudson/Downloads/gb-2013-14-10-r115-S24.txt")





### ORIGINAL AUTHOR: Andrew Teschendorff
# The original BMIQ function from Teschendorff 2013 adjusts for the type-2 bias in 
# Illumina Infinium 450k data.
# Later functions and edits were provided by yours truly, Steve Horvath.
# I changed the code so that one can calibrate methylation data to a gold standard.
# Specifically, I took version v_1.2 by Teschendorff  and fixed minor issues. 
# Also I made the code more robust e.g. by changing the optimization algorithm.
# Toward this end, I used the method="Nelder-Mead" in optim()

### Later functions and edits by Steve Horvath 
### # Steve Horvath took version v_1.2 by Teschendorff 
# and fixed minor errors. Also he made the code more robust.
# Importantly, SH changed the optimization algorithm to make it #more robust. 
# SH used method="Nelder-Mead" in optim() since the other #optimization method sometimes gets stuck.
#Toward this end, the function blc was replaced by blc2.



require(RPMM);



betaEst2=function (y, w, weights) 
{
    yobs = !is.na(y)
    if (sum(yobs) <= 1) 
        return(c(1, 1))
    y = y[yobs]
    w = w[yobs]
    weights = weights[yobs]
    N = sum(weights * w)
    p = sum(weights * w * y)/N
    v = sum(weights * w * y * y)/N - p * p
    logab = log(c(p, 1 - p)) + log(pmax(1e-06, p * (1 - p)/v - 
        1))
    if (sum(yobs) == 2) 
        return(exp(logab))
    opt = try(optim(logab, betaObjf, ydata = y, wdata = w, weights = weights, 
        method = "Nelder-Mead",control=list(maxit=50) ), silent = TRUE)
    if (inherits(opt, "try-error")) 
        return(c(1, 1))
    exp(opt$par)
} # end of function betaEst



blc2=function (Y, w, maxiter = 25, tol = 1e-06, weights = NULL, verbose = TRUE) 
{
    Ymn = min(Y[Y > 0], na.rm = TRUE)
    Ymx = max(Y[Y < 1], na.rm = TRUE)
    Y = pmax(Y, Ymn/2)
    Y = pmin(Y, 1 - (1 - Ymx)/2)
    Yobs = !is.na(Y)
    J = dim(Y)[2]
    K = dim(w)[2]
    n = dim(w)[1]
    if (n != dim(Y)[1]) 
        stop("Dimensions of w and Y do not agree")
    if (is.null(weights)) 
        weights = rep(1, n)
    mu = a = b = matrix(Inf, K, J)
    crit = Inf
    for (i in 1:maxiter) {
        warn0 = options()$warn
        options(warn = -1)
        eta = apply(weights * w, 2, sum)/sum(weights)
        mu0 = mu
        for (k in 1:K) {
            for (j in 1:J) {
                ab = betaEst2(Y[, j], w[, k], weights)
                a[k, j] = ab[1]
                b[k, j] = ab[2]
                mu[k, j] = ab[1]/sum(ab)
            }
        }
        ww = array(0, dim = c(n, J, K))
        for (k in 1:K) {
            for (j in 1:J) {
                ww[Yobs[, j], j, k] = dbeta(Y[Yobs[, j], j], 
                  a[k, j], b[k, j], log = TRUE)
            }
        }
        options(warn = warn0)
        w = apply(ww, c(1, 3), sum, na.rm = TRUE)
        wmax = apply(w, 1, max)
        for (k in 1:K) w[, k] = w[, k] - wmax
        w = t(eta * t(exp(w)))
        like = apply(w, 1, sum)
        w = (1/like) * w
        llike = weights * (log(like) + wmax)
        crit = max(abs(mu - mu0))
        if (verbose) 
            print(crit)
        if (crit < tol) 
            break
    }
    return(list(a = a, b = b, eta = eta, mu = mu, w = w, llike = sum(llike)))
}











# The function BMIQcalibration was created by Steve Horvath by heavily recycling code
# from A. Teschendorff's BMIQ function.
# BMIQ stands for beta mixture quantile normalization.
# Explanation: datM is a data frame with Illumina beta values (rows are samples, colums are CpGs.
# goldstandard is a numeric vector with beta values that is used as gold standard for calibrating the columns of datM.
# The length of goldstandard has to equal the number of columns of datM.
# Example code: First we impute missing values.
# library(WGCNA); dimnames1=dimnames(datMeth)
# datMeth= data.frame(t(impute.knn(as.matrix(t(datMeth)))$data))
# dimnames(datMeth)=dimnames1
# gold.mean=as.numeric(apply(datMeth,2,mean,na.rm=TRUE))
#datMethCalibrated=BMIQcalibration(datM=datMeth,goldstandard.beta=gold.mean)

BMIQcalibration=function(datM,goldstandard.beta,nL=3,doH=TRUE,nfit=20000,th1.v=c(0.2,0.75),th2.v=NULL,niter=5,tol=0.001,plots=FALSE,calibrateUnitInterval=TRUE){
if (length(goldstandard.beta) !=dim(datM)[[2]] ) {stop("Error in function arguments length(goldstandard.beta) !=dim(datM)[[2]]. Consider transposing datM.")}
if (plots ) {par(mfrow=c(2,2))}
beta1.v = goldstandard.beta

if (calibrateUnitInterval ) {datM=CalibrateUnitInterval(datM)}

### estimate initial weight matrix from type1 distribution
w0.m = matrix(0,nrow=length(beta1.v),ncol=nL);
w0.m[which(beta1.v <= th1.v[1]),1] = 1;
w0.m[intersect(which(beta1.v > th1.v[1]),which(beta1.v <= th1.v[2])),2] = 1;
w0.m[which(beta1.v > th1.v[2]),3] = 1;
### fit type1
print("Fitting EM beta mixture to goldstandard probes");
set.seed(1)
rand.idx = sample(1:length(beta1.v),min(c(nfit, length(beta1.v))  )   ,replace=FALSE)
em1.o = blc(matrix(beta1.v[rand.idx],ncol=1),w=w0.m[rand.idx,],maxiter=niter,tol=tol);
subsetclass1.v = apply(em1.o$w,1,which.max);
subsetth1.v = c(mean(max(beta1.v[rand.idx[subsetclass1.v==1]]),min(beta1.v[rand.idx[subsetclass1.v==2]])),mean(max(beta1.v[rand.idx[subsetclass1.v==2]]),min(beta1.v[rand.idx[subsetclass1.v==3]])));
class1.v = rep(2,length(beta1.v));
class1.v[which(beta1.v < subsetth1.v[1])] = 1;
class1.v[which(beta1.v > subsetth1.v[2])] = 3;
nth1.v = subsetth1.v;
print("Done");

### generate plot from estimated mixture
if(plots){
print("Check");
tmpL.v = as.vector(rmultinom(1:nL,length(beta1.v),prob=em1.o$eta));
tmpB.v = vector();
for(l in 1:nL){
  tmpB.v = c(tmpB.v,rbeta(tmpL.v[l],em1.o$a[l,1],em1.o$b[l,1]));
}
plot(density(beta1.v),main= paste("Type1fit-", sep=""));
d.o = density(tmpB.v);
points(d.o$x,d.o$y,col="green",type="l")
legend(x=0.5,y=3,legend=c("obs","fit"),fill=c("black","green"),bty="n");
}

### Estimate Modes 
if (  sum(class1.v==1)==1 ){ mod1U= beta1.v[class1.v==1]}
if (  sum(class1.v==3)==1 ){ mod1M= beta1.v[class1.v==3]}
if (  sum(class1.v==1) >1){ 
d1U.o = density(beta1.v[class1.v==1])
mod1U = d1U.o$x[which.max(d1U.o$y)]
}
if (  sum(class1.v==3)>1 ){ 
d1M.o = density(beta1.v[class1.v==3])
mod1M = d1M.o$x[which.max(d1M.o$y)]
}

### BETA 2
for (ii in 1:dim(datM)[[1]] ){
printFlush(paste("ii=",ii))
sampleID=ii
beta2.v = as.numeric(datM[ii,])

d2U.o = density(beta2.v[which(beta2.v<0.4)]);
d2M.o = density(beta2.v[which(beta2.v>0.6)]);
mod2U = d2U.o$x[which.max(d2U.o$y)]
mod2M = d2M.o$x[which.max(d2M.o$y)]

### now deal with type2 fit
th2.v = vector();
th2.v[1] = nth1.v[1] + (mod2U-mod1U);
th2.v[2] = nth1.v[2] + (mod2M-mod1M);

### estimate initial weight matrix 
w0.m = matrix(0,nrow=length(beta2.v),ncol=nL);
w0.m[which(beta2.v <= th2.v[1]),1] = 1;
w0.m[intersect(which(beta2.v > th2.v[1]),which(beta2.v <= th2.v[2])),2] = 1;
w0.m[which(beta2.v > th2.v[2]),3] = 1;

print("Fitting EM beta mixture to input probes");
# I fixed an error in the following line (replaced beta1 by beta2)
set.seed(1)
rand.idx = sample(1:length(beta2.v),min(c(nfit, length(beta2.v)),na.rm=TRUE)   ,replace=FALSE)
em2.o = blc2(Y=matrix(beta2.v[rand.idx],ncol=1),w=w0.m[rand.idx,],maxiter=niter,tol=tol,verbose=TRUE);
print("Done");

### for type II probes assign to state (unmethylated, hemi or full methylation)
subsetclass2.v = apply(em2.o$w,1,which.max);


if (sum(subsetclass2.v==2)>0 ){
subsetth2.v = c(mean(max(beta2.v[rand.idx[subsetclass2.v==1]]),min(beta2.v[rand.idx[subsetclass2.v==2]])),
mean(max(beta2.v[rand.idx[subsetclass2.v==2]]),min(beta2.v[rand.idx[subsetclass2.v==3]])));
}
if (sum(subsetclass2.v==2)==0 ){
subsetth2.v = c(1/2*max(beta2.v[rand.idx[subsetclass2.v==1]])+ 1/2*mean(beta2.v[rand.idx[subsetclass2.v==3]]), 1/3*max(beta2.v[rand.idx[subsetclass2.v==1]])+ 2/3*mean(beta2.v[rand.idx[subsetclass2.v==3]]));
}



class2.v = rep(2,length(beta2.v));
class2.v[which(beta2.v <= subsetth2.v[1])] = 1;
class2.v[which(beta2.v >= subsetth2.v[2])] = 3;

### generate plot
if(plots){
tmpL.v = as.vector(rmultinom(1:nL,length(beta2.v),prob=em2.o$eta));
tmpB.v = vector();
for(lt in 1:nL){
  tmpB.v = c(tmpB.v,rbeta(tmpL.v[lt],em2.o$a[lt,1],em2.o$b[lt,1]));
}
plot(density(beta2.v),  main= paste("Type2fit-",sampleID,sep="")  );
d.o = density(tmpB.v);
points(d.o$x,d.o$y,col="green",type="l")
legend(x=0.5,y=3,legend=c("obs","fit"),fill=c("black","green"),bty="n");
}

classAV1.v = vector();classAV2.v = vector();
for(l in 1:nL){
  classAV1.v[l] =  em1.o$mu[l,1];
  classAV2.v[l] =  em2.o$mu[l,1];
}

### start normalising input probes
print("Start normalising input probes");
nbeta2.v = beta2.v;
### select U probes
lt = 1;
selU.idx = which(class2.v==lt);
selUR.idx = selU.idx[which(beta2.v[selU.idx] > classAV2.v[lt])];
selUL.idx = selU.idx[which(beta2.v[selU.idx] < classAV2.v[lt])];
### find prob according to typeII distribution
p.v = pbeta(beta2.v[selUR.idx],em2.o$a[lt,1],em2.o$b[lt,1],lower.tail=FALSE);
### find corresponding quantile in type I distribution
q.v = qbeta(p.v,em1.o$a[lt,1],em1.o$b[lt,1],lower.tail=FALSE);
nbeta2.v[selUR.idx] = q.v;
p.v = pbeta(beta2.v[selUL.idx],em2.o$a[lt,1],em2.o$b[lt,1],lower.tail=TRUE);
### find corresponding quantile in type I distribution
q.v = qbeta(p.v,em1.o$a[lt,1],em1.o$b[lt,1],lower.tail=TRUE);
nbeta2.v[selUL.idx] = q.v;

### select M probes
lt = 3;
selM.idx = which(class2.v==lt);
selMR.idx = selM.idx[which(beta2.v[selM.idx] > classAV2.v[lt])];
selML.idx = selM.idx[which(beta2.v[selM.idx] < classAV2.v[lt])];
### find prob according to typeII distribution
p.v = pbeta(beta2.v[selMR.idx],em2.o$a[lt,1],em2.o$b[lt,1],lower.tail=FALSE);
### find corresponding quantile in type I distribution
q.v = qbeta(p.v,em1.o$a[lt,1],em1.o$b[lt,1],lower.tail=FALSE);
nbeta2.v[selMR.idx] = q.v;


if(doH){ ### if TRUE also correct type2 hemimethylated probes
### select H probes and include ML probes (left ML tail is not well described by a beta-distribution).
lt = 2;
selH.idx = c(which(class2.v==lt),selML.idx);
minH = min(beta2.v[selH.idx],na.rm=TRUE)
maxH = max(beta2.v[selH.idx],na.rm=TRUE)
deltaH = maxH - minH;
#### need to do some patching
deltaUH = -max(beta2.v[selU.idx],na.rm=TRUE) + min(beta2.v[selH.idx],na.rm=TRUE)
deltaHM = -max(beta2.v[selH.idx],na.rm=TRUE) + min(beta2.v[selMR.idx],na.rm=TRUE)

## new maximum of H probes should be
nmaxH = min(nbeta2.v[selMR.idx],na.rm=TRUE) - deltaHM;
## new minimum of H probes should be
nminH = max(nbeta2.v[selU.idx],na.rm=TRUE) + deltaUH;
ndeltaH = nmaxH - nminH;

### perform conformal transformation (shift+dilation)
## new_beta_H(i) = a + hf*(beta_H(i)-minH);
hf = ndeltaH/deltaH ;
### fix lower point first
nbeta2.v[selH.idx] = nminH + hf*(beta2.v[selH.idx]-minH);

}


### generate final plot to check normalisation
if(plots){
 print("Generating final plot");
 d1.o = density(beta1.v);
 d2.o = density(beta2.v);
 d2n.o = density(nbeta2.v);
 ymax = max(d2.o$y,d1.o$y,d2n.o$y);
plot(density(beta2.v),type="l",ylim=c(0,ymax),xlim=c(0,1), main=paste("CheckBMIQ-",sampleID,sep="") );
 points(d1.o$x,d1.o$y,col="red",type="l");
 points(d2n.o$x,d2n.o$y,col="blue",type="l");
 legend(x=0.5,y=ymax,legend=c("type1","type2","type2-BMIQ"),bty="n",fill=c("red","black","blue"));
}

datM[ii,]= nbeta2.v ;
} # end of for (ii=1 loop
datM
} # end of function BMIQcalibration





BMIQ = function(beta.v,design.v,nL=3,doH=TRUE,nfit=50000,th1.v=c(0.2,0.75),th2.v=NULL,niter=5,tol=0.001,plots=TRUE,sampleID=1,calibrateUnitInterval=TRUE){

if (calibrateUnitInterval) {
rangeBySample=range(beta.v,na.rm=TRUE)
minBySample=rangeBySample[1]
maxBySample=rangeBySample[2]
if ( (minBySample<0 | maxBySample>1) & !is.na(minBySample) & !is.na(maxBySample) ) {
y1=c(0.001,.999) 
x1=c(minBySample,maxBySample)
lm1=lm( y1 ~ x1 )
intercept1=coef(lm1)[[1]]
slope1=coef(lm1)[[2]]
beta.v=intercept1+slope1*beta.v
} # end of if
} # end of if (calibrateUnitInterval


type1.idx = which(design.v==1);
type2.idx = which(design.v==2);

beta1.v = beta.v[type1.idx];
beta2.v = beta.v[type2.idx];


### estimate initial weight matrix from type1 distribution
w0.m = matrix(0,nrow=length(beta1.v),ncol=nL);
w0.m[which(beta1.v <= th1.v[1]),1] = 1;
w0.m[intersect(which(beta1.v > th1.v[1]),which(beta1.v <= th1.v[2])),2] = 1;
w0.m[which(beta1.v > th1.v[2]),3] = 1;

### fit type1
print("Fitting EM beta mixture to goldstandard probes");
set.seed(1)
rand.idx = sample(1:length(beta1.v),min(c(nfit, length(beta1.v))  )   ,replace=FALSE)
em1.o = blc2(Y=matrix(beta1.v[rand.idx],ncol=1),w=w0.m[rand.idx,],maxiter=niter,tol=tol);
subsetclass1.v = apply(em1.o$w,1,which.max);
subsetth1.v = c(mean(max(beta1.v[rand.idx[subsetclass1.v==1]]),min(beta1.v[rand.idx[subsetclass1.v==2]])),mean(max(beta1.v[rand.idx[subsetclass1.v==2]]),min(beta1.v[rand.idx[subsetclass1.v==3]],na.rm=TRUE)));
class1.v = rep(2,length(beta1.v));
class1.v[which(beta1.v < subsetth1.v[1])] = 1;
class1.v[which(beta1.v > subsetth1.v[2])] = 3;
nth1.v = subsetth1.v;
print("Done");

### generate plot from estimated mixture
if(plots){
print("Check");
tmpL.v = as.vector(rmultinom(1:nL,length(beta1.v),prob=em1.o$eta));
tmpB.v = vector();
for(l in 1:nL){
  tmpB.v = c(tmpB.v,rbeta(tmpL.v[l],em1.o$a[l,1],em1.o$b[l,1]));
}

pdf(paste("Type1fit-",sampleID,".pdf",sep=""),width=6,height=4);
plot(density(beta1.v));
d.o = density(tmpB.v);
points(d.o$x,d.o$y,col="green",type="l")
legend(x=0.5,y=3,legend=c("obs","fit"),fill=c("black","green"),bty="n");
dev.off();
}



### Estimate Modes 
if (  sum(class1.v==1)==1 ){ mod1U= beta1.v[class1.v==1]}
if (  sum(class1.v==3)==1 ){ mod1M= beta1.v[class1.v==3]}
if (  sum(class1.v==1) >1){ 
d1U.o = density(beta1.v[class1.v==1])
mod1U = d1U.o$x[which.max(d1U.o$y)]
}
if (  sum(class1.v==3)>1 ){ 
d1M.o = density(beta1.v[class1.v==3])
mod1M = d1M.o$x[which.max(d1M.o$y)]
}


d2U.o = density(beta2.v[which(beta2.v<0.4)]);
d2M.o = density(beta2.v[which(beta2.v>0.6)]);
mod2U = d2U.o$x[which.max(d2U.o$y)]
mod2M = d2M.o$x[which.max(d2M.o$y)]


### now deal with type2 fit
th2.v = vector();
th2.v[1] = nth1.v[1] + (mod2U-mod1U);
th2.v[2] = nth1.v[2] + (mod2M-mod1M);

### estimate initial weight matrix 
w0.m = matrix(0,nrow=length(beta2.v),ncol=nL);
w0.m[which(beta2.v <= th2.v[1]),1] = 1;
w0.m[intersect(which(beta2.v > th2.v[1]),which(beta2.v <= th2.v[2])),2] = 1;
w0.m[which(beta2.v > th2.v[2]),3] = 1;

print("Fitting EM beta mixture to input probes");
set.seed(1)
rand.idx = sample(1:length(beta2.v),min(c(nfit, length(beta2.v)),na.rm=TRUE)   ,replace=FALSE)
em2.o = blc2(Y=matrix(beta2.v[rand.idx],ncol=1),w=w0.m[rand.idx,],maxiter=niter,tol=tol);
print("Done");

### for type II probes assign to state (unmethylated, hemi or full methylation)
subsetclass2.v = apply(em2.o$w,1,which.max);



if (sum(subsetclass2.v==2)>0 ){
subsetth2.v = c(mean(max(beta2.v[rand.idx[subsetclass2.v==1]]),min(beta2.v[rand.idx[subsetclass2.v==2]])),
mean(max(beta2.v[rand.idx[subsetclass2.v==2]]),min(beta2.v[rand.idx[subsetclass2.v==3]])));
}
if (sum(subsetclass2.v==2)==0 ){
subsetth2.v = c(1/2*max(beta2.v[rand.idx[subsetclass2.v==1]])+ 1/2*mean(beta2.v[rand.idx[subsetclass2.v==3]]), 1/3*max(beta2.v[rand.idx[subsetclass2.v==1]])+ 2/3*mean(beta2.v[rand.idx[subsetclass2.v==3]]));
}


class2.v = rep(2,length(beta2.v));
class2.v[which(beta2.v <= subsetth2.v[1])] = 1;
class2.v[which(beta2.v >= subsetth2.v[2])] = 3;


### generate plot
if(plots){
tmpL.v = as.vector(rmultinom(1:nL,length(beta2.v),prob=em2.o$eta));
tmpB.v = vector();
for(lt in 1:nL){
  tmpB.v = c(tmpB.v,rbeta(tmpL.v[lt],em2.o$a[lt,1],em2.o$b[lt,1]));
}
pdf(paste("Type2fit-",sampleID,".pdf",sep=""),width=6,height=4);
plot(density(beta2.v));
d.o = density(tmpB.v);
points(d.o$x,d.o$y,col="green",type="l")
legend(x=0.5,y=3,legend=c("obs","fit"),fill=c("black","green"),bty="n");
dev.off();
}

classAV1.v = vector();classAV2.v = vector();
for(l in 1:nL){
  classAV1.v[l] =  em1.o$mu[l,1];
  classAV2.v[l] =  em2.o$mu[l,1];
}

### start normalising input probes
print("Start normalising input probes");
nbeta2.v = beta2.v;
### select U probes
lt = 1;
selU.idx = which(class2.v==lt);
selUR.idx = selU.idx[which(beta2.v[selU.idx] > classAV2.v[lt])];
selUL.idx = selU.idx[which(beta2.v[selU.idx] < classAV2.v[lt])];
### find prob according to typeII distribution
p.v = pbeta(beta2.v[selUR.idx],em2.o$a[lt,1],em2.o$b[lt,1],lower.tail=FALSE);
### find corresponding quantile in type I distribution
q.v = qbeta(p.v,em1.o$a[lt,1],em1.o$b[lt,1],lower.tail=FALSE);
nbeta2.v[selUR.idx] = q.v;
p.v = pbeta(beta2.v[selUL.idx],em2.o$a[lt,1],em2.o$b[lt,1],lower.tail=TRUE);
### find corresponding quantile in type I distribution
q.v = qbeta(p.v,em1.o$a[lt,1],em1.o$b[lt,1],lower.tail=TRUE);
nbeta2.v[selUL.idx] = q.v;

### select M probes
lt = 3;
selM.idx = which(class2.v==lt);
selMR.idx = selM.idx[which(beta2.v[selM.idx] > classAV2.v[lt])];
selML.idx = selM.idx[which(beta2.v[selM.idx] < classAV2.v[lt])];
### find prob according to typeII distribution
p.v = pbeta(beta2.v[selMR.idx],em2.o$a[lt,1],em2.o$b[lt,1],lower.tail=FALSE);
### find corresponding quantile in type I distribution
q.v = qbeta(p.v,em1.o$a[lt,1],em1.o$b[lt,1],lower.tail=FALSE);
nbeta2.v[selMR.idx] = q.v;


if(doH){ ### if TRUE also correct type2 hemimethylated probes
### select H probes and include ML probes (left ML tail is not well described by a beta-distribution).
lt = 2;
selH.idx = c(which(class2.v==lt),selML.idx);
minH = min(beta2.v[selH.idx],na.rm=TRUE)
maxH = max(beta2.v[selH.idx],na.rm=TRUE)
deltaH = maxH - minH;
#### need to do some patching
deltaUH = -max(beta2.v[selU.idx],na.rm=TRUE) + min(beta2.v[selH.idx],na.rm=TRUE)
deltaHM = -max(beta2.v[selH.idx],na.rm=TRUE) + min(beta2.v[selMR.idx],na.rm=TRUE)

## new maximum of H probes should be
nmaxH = min(nbeta2.v[selMR.idx],na.rm=TRUE) - deltaHM;
## new minimum of H probes should be
nminH = max(nbeta2.v[selU.idx],na.rm=TRUE) + deltaUH;
ndeltaH = nmaxH - nminH;

### perform conformal transformation (shift+dilation)
## new_beta_H(i) = a + hf*(beta_H(i)-minH);
hf = ndeltaH/deltaH ;
### fix lower point first
nbeta2.v[selH.idx] = nminH + hf*(beta2.v[selH.idx]-minH);

}

pnbeta.v = beta.v;
pnbeta.v[type1.idx] = beta1.v;
pnbeta.v[type2.idx] = nbeta2.v;

### generate final plot to check normalisation
if(plots){
 print("Generating final plot");
 d1.o = density(beta1.v);
 d2.o = density(beta2.v);
 d2n.o = density(nbeta2.v);
 ymax = max(d2.o$y,d1.o$y,d2n.o$y);
 pdf(paste("CheckBMIQ-",sampleID,".pdf",sep=""),width=6,height=4)
 plot(density(beta2.v),type="l",ylim=c(0,ymax),xlim=c(0,1));
 points(d1.o$x,d1.o$y,col="red",type="l");
 points(d2n.o$x,d2n.o$y,col="blue",type="l");
 legend(x=0.5,y=ymax,legend=c("type1","type2","type2-BMIQ"),bty="n",fill=c("red","black","blue"));
 dev.off();
}

print(paste("Finished for sample ",sampleID,sep=""));

return(list(nbeta=pnbeta.v,class1=class1.v,class2=class2.v,av1=classAV1.v,av2=classAV2.v,hf=hf,th1=nth1.v,th2=th2.v));

}



CheckBMIQ = function(beta.v,design.v,pnbeta.v){### pnbeta is BMIQ normalised profile

type1.idx = which(design.v==1);
type2.idx = which(design.v==2);

beta1.v = beta.v[type1.idx];
beta2.v = beta.v[type2.idx];
pnbeta2.v = pnbeta.v[type2.idx];
  
} # end of function CheckBMIQ















CalibrateUnitInterval=function(datM,onlyIfOutside=TRUE){

rangeBySample=data.frame(lapply(data.frame(t(datM)),range,na.rm=TRUE))
minBySample=as.numeric(rangeBySample[1,])
maxBySample=as.numeric(rangeBySample[2,])
if (onlyIfOutside) { indexSamples=which((minBySample<0 | maxBySample>1) & !is.na(minBySample) & !is.na(maxBySample))
   }
if (!onlyIfOutside) { indexSamples=1:length(minBySample)}
if ( length(indexSamples)>=1 ){
for ( i in indexSamples) {
y1=c(0.001,0.999) 
x1=c(minBySample[i],maxBySample[i])
lm1=lm( y1 ~ x1 )
intercept1=coef(lm1)[[1]]
slope1=coef(lm1)[[2]]
datM[i,]=intercept1+slope1*datM[i,]
} # end of for loop
}
datM
} #end of function for calibrating to [0,1]





# STEP 4: Output the results 
if (  sum(  datout$Comment  != "" )   ==0 ) { cat(paste( "\n The individual samples appear to be fine. "),file="LogFile.txt",append=TRUE)  } 
if (  sum(  datout$Comment != "" )   >0 ) { cat(paste( "\n Warnings were generated for the following samples.\n", datout[,1][datout$Comment != ""], "\n Hint: Check the output file for more details."),file="LogFile.txt",append=TRUE)  } 
} 
# output the results into the directory
write.table(datout,"Output.csv", row.names=F, sep="," )
datout
```

Use horvath to validate our previously calculated epigenetic ages
```{r}
horvath_predictions <-read.csv('/Users/hugohudson/CAMH_projects/horvath_age_example.csv')
real_ages <- c(60,39,28,39,8,22,4,28,5,2,30,1,60,22, 8, 30)
horvath_df <- data.frame(horvath_predictions, real_ages)
epi_age_lm <- lm(horvath_df$real_ages ~ horvath_df$X0, data = horvath_df)
info <- read.csv('/Users/hugohudson/Downloads/gb-2013-14-10-r115-S27.csv')
horvath_df$diagnosis <- substr(info$title, 1, 7)

epigenetic_ages <- data.frame(filtered_epi$Epigenetic.Age)
predicted_ages <- predict(epi_age_lm, newdata = epigenetic_ages)
horvath_df$predicted_ages <- predicted_ages
horvath_calculated_ages <- c(60.00,43.00,28.00,38.00,8.20,20.00,4.80,38.00,6.80,3.60,31.00,0.98,62.00,24.00,8.00,43.00)
horvath_df$horvath_calculated_ages <- horvath_calculated_ages
his_calculations <- lm(horvath_df$real_ages ~ horvath_df$horvath_calculated_ages, data = horvath_df)
my_calculations <- lm(horvath_df$real_ages ~ horvath_df$predicted_ages, data = horvath_df)
his_plot <- ggplot(data = horvath_df, mapping = aes(x = real_ages, y = residuals(his_calculations), color = diagnosis)) + 
  geom_point() + 
  geom_hline(yintercept = 0) +
  labs(title = 'Horvath Clock Residuals', x='Chronological Age', y= 'Residuals')

my_plot <- ggplot(data = horvath_df, mapping = aes(x = real_ages, y = residuals(my_calculations), color = diagnosis)) + 
  geom_point() + 
  geom_hline(yintercept = 0) +
  labs(title = 'Custom Implementation Horvath Clock Residuals', x='Chronological Age', y= 'Residuals')

his_plot | my_plot
horvath_df
ggplot(horvath_df, aes(x = real_ages)) + 
  geom_smooth(horvath_df, mapping = aes(y = horvath_calculated_ages, color = 'horvath_calculated_ages')) +   # Line for y1
  geom_smooth(horvath_df, mapping = aes(y = predicted_ages, color = 'predicted_ages')) + 
  geom_point(horvath_df, mapping = aes(y = real_ages, color = 'real_ages')) +
  labs(title = "Clock Performance Comparison for example Cancer Data", x = "Chronological Ages (years)", y= "Epigenetic Age (years)")

```


Use horvath clock on CellDMC output
```{r}
read_file <- "/Users/hugohudson/Desktop/Lab/DMCT.csv"
dmctCSV <- read.csv(read_file)


```