#Installing 

# devtools::install_github(repo = "tiagochst/ELMER.data")

# devtools::install_github(repo = "tiagochst/ELMER")

if (!requireNamespace("BiocManager", quietly=TRUE))

  install.packages("BiocManager")

BiocManager::install("ELMER")

# BiocManager::install(version='devel')

# BiocManager::install("TCGAbiolinks")

# BiocManager::install("TCGAbiolinks", version = "2.5.9") 

# source("https://bioconductor.org/biocLite.R")

# biocLite("TCGAbiolinks")

library(ELMER)

library(MultiAssayExperiment)

library(ELMER.data)

library(sesameData)

library(tibble)

library(TCGAbiolinks)

#example

{

  distal.probes <- get.feature.probe(genome = "hg19",met.platform = "450K")

  head(distal.probes)

  data(LUSC_RNA_refined,package = "ELMER.data")

  data(LUSC_meth_refined,package = "ELMER.data")

  GeneExp[1:5,1:5]

  Meth[1:5,1:5]

  mae <- createMAE(exp = GeneExp, 

                   met = Meth,

                   save = TRUE,

                   linearize.exp = TRUE,

                   save.filename = "tmp.rda",

                   filter.probes = distal.probes,

                   met.platform = "450K",

                   genome = "hg19",

                   TCGA = TRUE)

  tmp = as.data.frame(colData(mae)[1:5,])

  View(tmp)

}

#data: RNA and methy

{

  setwd("E:/workplace/mywork/methy/dbgap/chf/data_chf_contr/early_chf/c1_UMN_JHU/RNA/")

  load(file="./Gene.Rdata")

  load(file="./exon.Rdata")

  exon[1:4,1:4];dim(exon)

  Gene[1:4,1:4];dim(Gene)

  exon$id <- paste("sample",exon$X,sep="_")

  exon <- column_to_rownames(exon,"id")

  exon <- exon[,-c(1:2)]

  Gene$id <- paste("sample",Gene$X,sep="_")

  Gene <- column_to_rownames(Gene,"id")

  Gene <- Gene[,-c(1:2)]

  setwd("E:\\workplace\\mywork\\methy\\dbgap\\chf\\data_chf_contr\\early_chf\\c1_UMN_JHU\\train_UMN_tset_JHU/1123_dataSummary/new_champ/")

  load("CorrectedBeta.Rdata")

  CorrectedBeta[1:4,1:4]

  dim(CorrectedBeta)

  data_cpg <- read.csv("../champ_rawdata_train.csv")

  colnames(CorrectedBeta) <- paste("sample",data_cpg$shareid,sep="_")

  id = setdiff(rownames(exon),colnames(CorrectedBeta))

  exon_data = exon[!(rownames(exon) %in% id),]

  exon_data <- t(exon_data)

  exon_data[1:4,1:4];dim(exon_data)

  Gene_data = Gene[!(rownames(Gene) %in% id),]

  Gene_data <- t(Gene_data)

  Gene_data[1:4,1:4];dim(Gene_data)

  CorrectedBeta <- subset(CorrectedBeta,select = colnames(exon_data))

  dim(CorrectedBeta);dim(exon_data);dim(Gene_data)

  CorrectedBeta[1:4,1:4]

  exon_data[1:4,1:4]

  Gene_data[1:4,1:4]

  save(CorrectedBeta,exon_data,Gene_data,file="result/MAE.Rdata")

  {

    library(stringr) 

    exon_data[1:4,1:5]

    ids=data.frame(ensembl_id=str_split(rownames(exon_data),'[.]',simplify = T)[,1],

                   median=apply(exon_data,1,median))

    head(ids)

    head(ids$ensembl_id)

    library(org.Hs.eg.db)

    g2s=unique(toTable(org.Hs.egSYMBOL))

    head(g2s)

    g2e=unique(toTable(org.Hs.egENSEMBL)) 

    head(g2e)

    s2e=merge(g2e,g2s,by='gene_id')

    head(s2e)

    table(ids$ensembl_id %in% s2e$symbol)

    # FALSE  TRUE 

    # 2866 15448 

    ids=ids[ids$ensembl_id %in% s2e$symbol,]

    ids$ENSEMBL=s2e[match(ids$ensembl_id,s2e$symbol),2]

    length(unique(ids$ENSEMBL))

    head(ids) 

    ids=ids[order(ids$ENSEMBL,ids$median,decreasing = T),]

    ids=ids[!duplicated(ids$ENSEMBL),]

    dim(ids) 

    exon_data = as.data.frame(exon_data)

    exon_data = rownames_to_column(exon_data,"ensembl_id")

    exon_data = merge(ids,exon_data,by="ensembl_id")

    exon_data = column_to_rownames(exon_data,"ENSEMBL")

    exon_data = exon_data[,-c(1,2)]

    exon_data[1:4,1:4]  

    dim(exon_data)

  }

  {

    library(stringr) 

    Gene_data[1:4,1:5]

    ids=data.frame(ensembl_id=str_split(rownames(Gene_data),'[.]',simplify = T)[,1],

                   median=apply(Gene_data,1,median))

    head(ids)

    head(ids$ensembl_id)

    library(org.Hs.eg.db)

    g2s=unique(toTable(org.Hs.egSYMBOL))

    head(g2s)

    g2e=unique(toTable(org.Hs.egENSEMBL)) 

    head(g2e)

    s2e=merge(g2e,g2s,by='gene_id')

    head(s2e)

    table(ids$ensembl_id %in% s2e$symbol)

    # FALSE  TRUE 

    # 2529 14870 

    ids=ids[ids$ensembl_id %in% s2e$symbol,]

    ids$ENSEMBL=s2e[match(ids$ensembl_id,s2e$symbol),2]

    length(unique(ids$ENSEMBL))

    head(ids) 

    ids=ids[order(ids$ENSEMBL,ids$median,decreasing = T),]

    ids=ids[!duplicated(ids$ENSEMBL),]

    dim(ids) 

    Gene_data = as.data.frame(Gene_data)

    Gene_data = rownames_to_column(Gene_data,"ensembl_id")

    Gene_data = merge(ids,Gene_data,by="ensembl_id")

    Gene_data = column_to_rownames(Gene_data,"ENSEMBL")

    Gene_data = Gene_data[,-c(1,2)]

    Gene_data[1:4,1:4]  

    dim(Gene_data)

  }

  dim(Gene_data);dim(exon_data)

  save(exon_data,Gene_data,file="result/MAE_addENSEMBL.Rdata")

}

#creat mea

{

  load("result/MAE_addENSEMBL.Rdata") #exon_data,Gene_data

  distal.probes <- get.feature.probe(genome = "hg19",met.platform = "450K")

  assay <- c(rep("DNA methylation", ncol(CorrectedBeta)),

             rep("Gene expression", ncol(exon_data)))

  primary <- c(colnames(CorrectedBeta),colnames(exon_data))

  data_cpg <- read.csv("../champ_rawdata_train.csv")

  data_cpg$shareid <- paste("sample",data_cpg$shareid,sep="_")

  colname <- data_cpg[colnames(CorrectedBeta) %in% data_cpg$shareid ,"Sample_Group"]

  colname <- c(colnames(CorrectedBeta),colnames(exon_data))

  sampleMap <- data.frame(assay,primary,colname)

  #distal.probes <- get.feature.probe(genome = "hg19",met.platform = "450K")

  colData <- data.frame(sample = colnames(CorrectedBeta))

  tmp = data_cpg[data_cpg$shareid %in% colData$sample ,]

  colData <- merge(colData,tmp[,c(1,11)],by.x="sample",by.y="shareid")

  colData = unique(colData)

  rownames(colData) <- colData$sample

  mae <- createMAE(exp = exon_data, 

                   met = CorrectedBeta,

                   save = TRUE,

                   filter.probes = distal.probes,

                   colData = colData,

                   #sampleMap = sampleMap,

                   linearize.exp = TRUE,

                   save.filename = "mae.rda",

                   met.platform = "450K",

                   genome = "hg19",

                   TCGA = FALSE)

  save(mae,"result/mae.rda")

}

#DMP

{

  mae <- get(load("result/mae.rda"))

  sig.diff <- get.diff.meth(data = mae, 

                            group.col = "Sample_Group",

                            group1 =  "chf",

                            group2 = "nochf",

                            # if supervised mode set to 1

                            mode = "supervised",

                            minSubgroupFrac = 1, 

                            sig.dif = 0,

                            diff.dir = "both", # “hypo” and “hyper” 

                            cores = 1, 

                            dir.out ="result", 

                            save =  FALSE ,

                            pvalue = 0.05#0.05#1

                            )

  head(sig.diff);

  dim(sig.diff)#pvalue = 0.05,19;0.5,8537

  write.table(sig.diff,"result/getMethdiff.both.probes.csv")#pvalue = 0.5,8537

  write.table(sig.diff,"result/getMethdiff.both.probes.significant.csv")#pvalue = 0.05,19

  ##"getMethdiff.hypo.probes.csv"            

  ##"getMethdiff.hypo.probes.significant.csv"

  save(sig.diff,file="result/sig_diff.Rdata")#pvalue = 0.5,8537

}

#Identifying putative probe-gene pairs

{

  setwd("E:\\workplace\\mywork\\methy\\dbgap\\chf\\data_chf_contr\\early_chf\\c1_UMN_JHU\\train_UMN_tset_JHU/1123_dataSummary/new_champ/")

  mae <- get(load("result/mae.rda"))

  #sig.diff <- read.csv("result/getMethdiff.both.probes.significant.csv",sep = " ")

  nearGenes <- GetNearGenes(data = mae, 

                            probes = sig.diff$probe, 

                            numFlankingGenes = 20) # 10 upstream and 10 dowstream genes

  save(nearGenes,file="result/nearGenes.Rdata")

  Hypo.pair <- get.pair(data = mae,

                        group.col = "Sample_Group",

                        group1 =  "chf",

                        group2 = "nochf",

                        nearGenes = nearGenes,

                        #minSubgroupFrac = 1,

                        mode = "unsupervised",#supervised or unsupervised

                        permu.dir = "result/permu",

                        #Please set to 100000 to get significant results

                        #permu.size = 100000,

                        raw.pvalue = 0.05,#1,

                        #Please set to 0.001 to get significant results

                        Pe = 0.05,#1, 

                        cores = 1,

                        #See preAssociationProbeFiltering function

                        filter.probes = FALSE, 

                        # filter.percentage = 0.01,

                        # filter.portion = 0.3,

                        dir.out = "result",

                        label = "both",

                        save = FALSE)

  dim(Hypo.pair)

  write.table(Hypo.pair,"result/getPair.both.all.pairs.statistic5.csv")

  write.table(Hypo.pair,"result/getPair.both.all.pairs.statistic.significant.csv")

  #"getPair.hypo.all.pairs.statistic.csv"                  

  #"getPair.hypo.pairs.significant.csv"                    

  #"getPair.hypo.pairs.statistic.with.empirical.pvalue.csv"

}

#motif

{

  # Load results from previous sections

  mae <- get(load("result/mae.rda"))

  Hypo.pair1 = read.csv("result/getPair.both.all.pairs.statistic1.csv",sep=" ")

  Hypo.pair2 = read.csv("result/getPair.both.all.pairs.statistic2.csv",sep=" ")

  Hypo.pair3 = read.csv("result/getPair.both.all.pairs.statistic3.csv",sep=" ")

  Hypo.pair4 = read.csv("result/getPair.both.all.pairs.statistic4.csv",sep=" ")

  Hypo.pair5 = read.csv("result/getPair.both.all.pairs.statistic5.csv",sep=" ")

  pair <- rbind(rbind(rbind(rbind(Hypo.pair1,Hypo.pair2),Hypo.pair3),Hypo.pair4),Hypo.pair5)

  head(pair) # significantly hypomethylated probes with putative target genes

  pair = Hypo.pair

  # Identify enriched motif for significantly hypomethylated probes which 

  # have putative target genes.

  enriched.motif <- get.enriched.motif(data = mae,

                                       #probes = pair$Probe, 

                                       probes = sig.diff$probe, 

                                       dir.out = "result", 

                                       label = "both",

                                       min.incidence = 0,

                                       lower.OR = 1,

                                       save = FALSE)

  names(enriched.motif)

  head(enriched.motif[names(enriched.motif)[1]]) ## probes in the given set that have the first motif.

  write.table(enriched.motif$FOSL2_HUMAN.H11MO.0.A,"result/getMotif.both.motif.enrichment.txt")

  save(enriched.motif,file="result/getMotif.both.enriched.motifs.rda")

  #"getMotif.hypo.enriched.motifs.rda"  "getMotif.both.motif.enrichment.csv" "motif.enrichment.pdf") 

}

#TF

{

  # Load results from previous sections

  mae <- get(load("mae.rda"))

  load("result/getMotif.both.enriched.motifs12.rda")

  TF <- get.TFs(data = mae, 

                group.col = "Sample_Group",

                group1 =  "chf",

                group2 = "nochf",

                mode = "unsupervised",

                enriched.motif = enriched.motif,

                dir.out = "result", 

                cores = 1, 

                label = "both",

                save =  FALSE )

  ##"getTF.hypo.TFs.with.motif.pvalue.rda"             

  ##"getTF.hypo.significant.TFs.with.motif.summary.csv"

  write.table(TF,"result/getTF.both.significant.TFs.with.motif.summary.csv")

  save(TF,file="result/getTF.both.TFs.with.motif.pvalue.rda")

}

#Scatter plots

{

  mae <- get(load("mae.rda"))

  load("result/getMotif.hypo.enriched.motifs.rda")

  scatter.plot(data = mae,

               byProbe = list(probe = c("cg27401945"), numFlankingGenes = 20), 

               category = "Sample_Group", 

               lm = TRUE, # Draw linear regression curve

               save = TRUE) 

  scatter.plot(data = mae,

               byPair = list(probe = c("cg27401945"), gene = c("ENSG00000148704")), 

               category = "Sample_Group", save = TRUE, lm_line = TRUE) 

  load("result/getMotif.hypo.enriched.motifs.rda")

  names(enriched.motif)[1]

  scatter.plot(data = mae,

               byTF = list(TF = c("VAX1","DIP2C"),

                           probe = enriched.motif[[names(enriched.motif)[2]]]), 

               category = "Sample_Group",

               save = TRUE, 

               lm_line = TRUE)

}

#XX 

{

  # Load results from previous sections

  mae <- get(load("mae.rda"))

  #pair <- read.csv("result/getPair.hypo.pairs.significant.csv")

  schematic.plot(pair = pair, 

                 data = mae,

                 group.col = "Sample_Group",

                 byProbe = pair$Probe[1],

                 save = FALSE)

  schematic.plot(pair = pair, 

                 data = mae,   

                 group.col = "Sample_Group", 

                 byGene = pair$GeneID[1],

                 save = FALSE)

}

#XX motif

{

  motif.enrichment.plot(motif.enrichment = enriched.motif,

                        significant = list(OR = 1.5,lowerOR = 1.3),

                        label = "both",

                        save = FALSE)

  motif.enrichment.plot(motif.enrichment = "result/getMotif.hypo.motif.enrichment.csv",

                        significant = list(OR = 1.5,lowerOR = 1.3),

                        label = "hypo",

                        summary = TRUE,

                        save = FALSE)  

}

#TF

{

  load("result/getTF.hypo.TFs.with.motif.pvalue.rda")

  motif <- colnames(TF.meth.cor)[1]

  TF.rank.plot(motif.pvalue = TF.meth.cor,

               motif = motif,

               save = FALSE)

}

#XX heatmap

{

  # Load results from previous sections

  mae <- get(load("mae.rda"))

  pair <- read.csv("result/getPair.hypo.pairs.significant.csv")

  heatmapPairs(data = mae, 

               group.col = "Sample_Group",

               group1 = "Chf", 

               #annotation.col = c("years_smoked","gender"),

               group2 = "Nochf",

               pairs = pair,

               filename =  NULL)

}

{

  library(plyr)

  A<-strsplit(as.character(names(enriched.motif)), "_")

  tmp2 <- mapply( cbind,  A )

  df0 <- ldply (tmp2[1,], data.frame)

  #median level of methylation estimated from all distal probes

  length(enriched.motif$FOSL2_HUMAN.H11MO.0.A)

  FOSL2 = enriched.motif$FOSL2_HUMAN.H11MO.0.A

  FOSL1 = enriched.motif$FOSL1_HUMAN.H11MO.0.A

}