Commit 3ee7e4c1 authored by Marek Ostaszewski's avatar Marek Ostaszewski
Browse files

configuration cap on retrieved maps and pathways

parent a9500e72
type resource
map https://pdmap.uni.lu/minerva/api/
pathway WikiPathways_2019_Human
pathway Panther_2016
type resource value parameter max_areas_per_map 5 parameter max_areas_per_pathway_db 5 map https://pdmap.uni.lu/minerva/api/ map https://progeria.uni.lu/minerva/api/ pathway_db WikiPathways_2019_Human pathway_db Panther_2016
\ No newline at end of file
......
......@@ -64,11 +64,11 @@ for(map in config[config$type == "map","resource"]) {
message(paste0("Going through individual models..."))
for(i in 1:length(model_elements)) {
me <- model_elements[[i]]
pVals <- c()
pnames <- unique(me[me$type == "Pathway","name"])
if(length(unique(me[me$type == "Pathway","name"])) == 0) { next }
path_pVals <- data.frame(pname = unique(me[me$type == "Pathway","name"]), pVal = 1)
### Calculating the inclusion of the elements per bioentity set
for(up in pnames) {
res <- apply(me[me$name == up,"bounds"], 1,
for(up in 1:nrow(path_pVals)) {
res <- apply(me[me$name == path_pVals$pname[up],"bounds"], 1,
function(x) me$bounds$x >= x["x"] & me$bounds$x <= (x["x"] + x["width"]) &
me$bounds$y >= x["y"] & me$bounds$y <= (x["y"] + x["height"]))
......@@ -78,22 +78,28 @@ for(map in config[config$type == "map","resource"]) {
k <- length(intersect(input,hgncs))
if (k > 0) {
pVal <- (((m * (n - k)) / (k * (N - m))) ** k) * (((n * (N - m)) / (N * (n - k))) ** n)
pVals <- c(pVals, pVal)
} else {
pVals <- c(pVals, 1)
path_pVals$pVal[up] <- (((m * (n - k)) / (k * (N - m))) ** k) * (((n * (N - m)) / (N * (n - k))) ** n)
}
}
### Bonferroni correction
# pVals <- pVals*length(pVals)
if(length(pnames[pVals < 0.1])) {
pbounds <- cbind(name = me[me$name %in% pnames[pVals < 0.1],]$name,
me[me$name %in% pnames[pVals < 0.1],]$bounds)
# path_pVals$pVal <- path_pVals$pVal*nrow(path_pVals)
path_pVals <- path_pVals[path_pVals$pVal < 0.1, ]
if(nrow(path_pVals) > 0) {
path_pVals <- path_pVals[order(path_pVals$pVal), ]
if(!is.null(config[config$resource == "max_areas_per_map",]$value)) {
path_pVals <- path_pVals[1:config[config$resource == "max_areas_per_map",]$value, ]
}
pbounds <- cbind(name = me[me$name %in% path_pVals$pname,]$name,
me[me$name %in% path_pVals$pname,]$bounds)
pathways_to_pull <- cbind(minerva_instance = map,
project_id = project_id,
model_id = models$idObject[i],
pbounds, stringsAsFactors = F)
write.table(pathways_to_pull, file = "enriched_disease_maps.txt",
sep = "\t", quote = F, col.names = F, row.names = F, append = T)
message(paste0("Output file updated..."))
......@@ -110,14 +116,23 @@ cat(file = "enriched_pathways.txt",
"pathway_db\tTerm\tAdjusted.P.value\n",
append = F)
pathway_dbs <- config[config$type == "pathway","resource"]
pathway_dbs <- config[config$type == "pathway_db","resource"]
enriched_pathways <- enrichr(input, databases = pathway_dbs)
adj_pathways <- lapply(enriched_pathways, function(x) x[x$Adjusted.P.value < 0.1, ])
adj_pathways <- lapply(enriched_pathways,
function(x) {
selected <- which(x$Adjusted.P.value < 0.1)
### Cap the amount of pathways to the predefined max in config
cap <- min(length(selected), config[config$resource == "max_areas_per_pathway_db",]$value)
selected <- head(selected, cap)
x[selected, ]
})
for(aep in 1:length(adj_pathways)) {
write.table(cbind(pathway_db = pathway_dbs[aep], adj_pathways[[aep]][,c("Term","Adjusted.P.value")], stringsAsFactors = F),
write.table(cbind(pathway_db = pathway_dbs[aep],
adj_pathways[[aep]][,c("Term","Adjusted.P.value")],
stringsAsFactors = F),
file = "enriched_pathways.txt", sep = "\t", quote = F, col.names = F, row.names = F, append = T)
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment