Skip to content

Example Script

Ewoud Ewing edited this page Jul 17, 2020 · 2 revisions

require(devtools)
install_github("https://github.com/TranslationalBioinformaticsUnit/GeneSetCluster/")
require(GeneSetCluster)

Canonical.files <- c(system.file("extdata", "MM10.IPA.KO.uGvsMac.Canonical_pathways.xls", package = "GeneSetCluster"),
                     system.file("extdata", "MM10.IPA.WT.uGvsMac.Canonical_pathways.xls", package = "GeneSetCluster"))

require(readxl)
MM10.IPA.KO.uGvsMac.Canonical <- read_excel(path = Canonical.files[1], skip = 1, sheet = 1)
MM10.IPA.WT.uGvsMac.Canonical <- read_excel(path = Canonical.files[2], skip = 1, sheet = 1)

MM10.IPA.KO.uGvsMac.Canonical <- as.data.frame(MM10.IPA.KO.uGvsMac.Canonical)
MM10.IPA.WT.uGvsMac.Canonical <- as.data.frame(MM10.IPA.WT.uGvsMac.Canonical)

require(limma)
MM10.IPA.KO.uGvsMac.Canonical$MoleculesCount <- NA

for(can.i in 1:nrow(MM10.IPA.KO.uGvsMac.Canonical))
{
  mol <- as.vector(strsplit2(MM10.IPA.KO.uGvsMac.Canonical[can.i,"Molecules"], split=","))
  MM10.IPA.KO.uGvsMac.Canonical[can.i,"MoleculesCount"]  <- length(mol)
}

head(MM10.IPA.KO.uGvsMac.Canonical)

MM10.IPA.WT.uGvsMac.Canonical$MoleculesCount <- NA

for(can.i in 1:nrow(MM10.IPA.WT.uGvsMac.Canonical))
{
  mol <- as.vector(strsplit2(MM10.IPA.WT.uGvsMac.Canonical[can.i,"Molecules"], split=","))
  MM10.IPA.WT.uGvsMac.Canonical[can.i,"MoleculesCount"]  <- length(mol)
}

####################

MM10.IPA.KO.uGvsMac.Canonical.sig <- MM10.IPA.KO.uGvsMac.Canonical[MM10.IPA.KO.uGvsMac.Canonical$`-log(p-value)` > 1.31 & MM10.IPA.KO.uGvsMac.Canonical$MoleculesCount > 5,]
MM10.IPA.WT.uGvsMac.Canonical.sig <- MM10.IPA.WT.uGvsMac.Canonical[MM10.IPA.WT.uGvsMac.Canonical$`-log(p-value)` > 1.31 & MM10.IPA.WT.uGvsMac.Canonical$MoleculesCount > 5,]

####################

IPA.KOvsWT <- ObjectCreator(Pathways = c(MM10.IPA.KO.uGvsMac.Canonical.sig$`Ingenuity Canonical Pathways`, 
                                         MM10.IPA.WT.uGvsMac.Canonical.sig$`Ingenuity Canonical Pathways`), 
                            Molecules = c(MM10.IPA.KO.uGvsMac.Canonical.sig$Molecules, 
                                          MM10.IPA.WT.uGvsMac.Canonical.sig$Molecules),
                            Groups = c(rep("KO", times=nrow(MM10.IPA.KO.uGvsMac.Canonical.sig)),
                                       rep("WT", times=nrow(MM10.IPA.WT.uGvsMac.Canonical.sig))), 
                            Source = "IPA", Type = "Canonical", structure = "SYMBOL", sep = ",", organism = "org.Mm.eg.db")

IPA.KOvsWT <- CombineGeneSets(IPA.KOvsWT)

require(ggplot2)
OptimalGeneSets(object = IPA.KOvsWT, method = "gap", max_cluster = 15, cluster_method = "kmeans", main = "IPA.KOvsWT")

IPA.KOvsWT <- ClusterGeneSets(Object = IPA.KOvsWT, clusters = 6,method = "kmeans", order = "cluster" )

PlotGeneNetworks(IPA.KOvsWT, labels = F, RRmin = 5)
PlotGeneSets(Object = IPA.KOvsWT, annotation.mol = F, main = "IPA.KOvsWT with 10 kmeans clusters")
PlotGeneSets(Object = IPA.KOvsWT, annotation.mol = F, main = "IPA.KOvsWT with 10 kmeans clusters", RR.max = 90, cluster.order = c(1,8,3,4,5,6,7,9,10,2))

########################
require(WebGestaltR)
IPA.KOvsWT.ORA.all <- ORAGeneSets(Object = IPA.KOvsWT, ORA.returned = 10, unique.per.cluster = F)

View(IPA.KOvsWT.ORA.all)



########################
IPA.KOvsWT.String.all <- GetSTRINGdbPerGeneSets(Object = IPA.KOvsWT, unique.per.cluster = F)

par(mfrow=c(2,3))

plotSTRINGdbPerGeneSets(StringObject = IPA.KOvsWT.String.all, plot.cluster = 1)
plotSTRINGdbPerGeneSets(StringObject = IPA.KOvsWT.String.all, plot.cluster = 2)
plotSTRINGdbPerGeneSets(StringObject = IPA.KOvsWT.String.all, plot.cluster = 3)
plotSTRINGdbPerGeneSets(StringObject = IPA.KOvsWT.String.all, plot.cluster = 4)
plotSTRINGdbPerGeneSets(StringObject = IPA.KOvsWT.String.all, plot.cluster = 5)
plotSTRINGdbPerGeneSets(StringObject = IPA.KOvsWT.String.all, plot.cluster = 6)


Clone this wiki locally