load(file = "./DNAme_Retina.RData")
library(limma)
library(RColorBrewer)
library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
library(IlluminaHumanMethylation450kmanifest)
pal.type <- brewer.pal(8, "Paired")
ann450k <- getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)

########DMP#############
mVals <- log2(bVals/(1-bVals))
##DMP between normal and RB
cellType <- factor(target$type)
# use the above to create a design matrix
design <- model.matrix(~0+cellType, data=target)
colnames(design) <- c(levels(cellType))
# fit the linear model 
fit <- lmFit(mVals, design)
# create a contrast matrix for specific comparisons
contMatrix <- makeContrasts(RB-retina,
                            #RB-DBT,
                            #RB-FVM,
                            #RB-PVR,
                            levels=design)
contMatrix
fit2 <- contrasts.fit(fit, contMatrix)
fit2 <- eBayes(fit2)
# look at the numbers of DM CpGs at FDR < 0.05
summary(decideTests(fit2))
ann450kSub <- ann450k[match(rownames(mVals),ann450k$Name),
                      c(1:4,12:19,24:ncol(ann450k))]
DMPs <- topTable(fit2, num=Inf, coef=1, genelist=ann450kSub)
head(DMPs)
DMPs <- DMPs[DMPs$adj.P.Val< 0.01 & DMPs$logFC > 2 ,] # only hypermethylated
DMP_RB.cg <- rownames(DMPs)
DMP_RB.cg <- DMP_RB.cg[complete.cases(DMP_RB.cg)]# remove NAs 

table(target$type)
bVals <- bVals[DMP_RB.cg,]
idx.RB <- which(target$type == "RB")
idx.retina <- which(target$type == "retina")
idx.DBT <- which(target$type == "DBT")
idx.FVM <- which(target$type == "FVM")

DMP_keep <- which(rowMeans(bVals[,idx.RB]) > 0.5 & rowMeans(bVals[,idx.DBT]) <0.2 & rowMeans(bVals[,idx.FVM]) <0.2)
DMP_keep <- rownames(bVals[DMP_keep,])

library(pheatmap)
col.type <- pal.type
annotation_col = data.frame(
  Type = target$type
)
rownames(annotation_col) = colnames(bVals)
names(col.type) <- levels(as.factor(target$type))
breaksList = seq(0, 1, by = 0.01)

pdf("./plots/101Heatmap_DMP_normal.pdf", width = 5, height = 5)
pheatmap(bVals[DMP_keep,grep("Embryo|PVR", target$type, invert = T)], 
         scale = "none",
         clustering_distance_cols = "euclidean",
         clustering_method = "complete",
         show_rownames = T, show_colnames = F,
         fontsize_row = 10, fontsize_col = 1,
         cluster_cols = T,
         cluster_rows = T,
         annotation_legend = T,
         color = colorRampPalette(c("darkblue", "yellow"))(length(breaksList)), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
         border_color = NA,
         annotation_col = annotation_col,
         annotation_colors  = list(
           #Time = c("white", "firebrick"),
           # subtype = c(subtype1 = "#A6CEE3", subtype2 = "#1F78B4",
           #            subtype3 = "#B2DF8A", subtype4 = "#33A02C"),
           Type = col.type[c(1,3:4,6:8)]
         ),
         main = "Methylation Levels \nBetween Types"
)
dev.off()


########boxplot between types#########
t(bVals[DMP_keep,]) #select promoter cg
RB.marker <- data.frame(accession = target$geo_accession, 
                        tissue = target$type, 
                        bVal = t(bVals[DMP_keep,])[,1], 
                        bVal = t(bVals[DMP_keep,])[,2], 
                        stemness = target$stemness,
                        stringsAsFactors = F)
ann450kSub[DMP_keep,]
#chr6:10391237-10391412
colnames(RB.marker)[3] <- "cg17754510"
colnames(RB.marker)[4] <- "cg21995304"
RB.marker <- RB.marker[grep("Embryo|PVR", RB.marker$tissue, invert = T),]

pdf("./plots/101Boxplot_RB_marker.pdf", width = 5, height = 5)
boxplot(cg17754510 ~ tissue , data = RB.marker , col = pal.type, 
        ylim = c(0,1), notch = F,  las = 2,
        main = "cg17754510")
boxplot(cg21995304 ~ tissue , data = RB.marker , col = pal.type, 
        ylim = c(0,1), notch = F,  las = 2,
        main = "cg21995304")
dev.off()

cor(RB.marker$cg17754510, RB.marker$stemness)
cor.test(RB.marker$cg17754510, RB.marker$stemness)


##########ROC###########
library(pROC)

pdf("./plots/101ROC_RBmarker.pdf", width = 4, height = 4)
RB.marker$outcome <- RB.marker$tissue
RB.marker$outcome <- NA
RB.marker$outcome[RB.marker$tissue == "RB"] <- 1
RB.marker$outcome[RB.marker$tissue == "retina"] <- 0
modelroc <- roc(RB.marker$outcome, RB.marker$cg17754510)
plot(modelroc, print.auc=TRUE, auc.polygon=F, grid=c(0.1, 0.2),
     #grid.col=c("green", "red"),
     max.auc.polygon=F,
     #auc.polygon.col="skyblue", 
     print.thres=F,
     main = "ROC for cg17754510 (RB vs retina)")

modelroc <- roc(RB.marker$outcome, RB.marker$cg21995304)
plot(modelroc, print.auc=TRUE, auc.polygon=F, grid=c(0.1, 0.2),
     #grid.col=c("green", "red"),
     max.auc.polygon=F,
     #auc.polygon.col="skyblue", 
     print.thres=F,
     main = "ROC for cg21995304 (RB vs retina)")


RB.marker$outcome <- RB.marker$tissue
RB.marker$outcome <- NA
RB.marker$outcome[RB.marker$tissue == "RB"] <- 1
RB.marker$outcome[RB.marker$tissue == "DBT"] <- 0
modelroc <- roc(RB.marker$outcome, RB.marker$cg17754510)
plot(modelroc, print.auc=TRUE, auc.polygon=F, grid=c(0.1, 0.2),
     #grid.col=c("green", "red"),
     max.auc.polygon=F,
     #auc.polygon.col="skyblue", 
     print.thres=F,
     main = "ROC for cg17754510 (RB vs DBT)")

modelroc <- roc(RB.marker$outcome, RB.marker$cg21995304)
plot(modelroc, print.auc=TRUE, auc.polygon=F, grid=c(0.1, 0.2),
     #grid.col=c("green", "red"),
     max.auc.polygon=F,
     #auc.polygon.col="skyblue", 
     print.thres=F,
     main = "ROC for cg21995304 (RB vs DBT)")
dev.off()


###genomic location #####
ann450kSub[DMP_keep,]
#chr6:10391237-10391412