diff --git a/NAMESPACE b/NAMESPACE index f834807..7fa691c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -36,7 +36,7 @@ export(loadMIDAS) export(loadMIDAS2) export(loadMIDAS2Legacy) export(loadmetaSNV) -export(obtain_positions) +export(obtainPositions) export(plotCirclize) export(plotCoverage) export(plotGenes) @@ -46,7 +46,7 @@ export(plotMAF) export(plotPCA) export(plotTree) export(returnGenes) -export(reverse_annot) +export(reverseAnnot) export(setAnnotation) export(setGroup) export(setMetadata) diff --git a/R/checkEGGNOG.R b/R/checkEGGNOG.R index 47df030..25981e9 100644 --- a/R/checkEGGNOG.R +++ b/R/checkEGGNOG.R @@ -222,5 +222,6 @@ summariseAbundance <- function(stana, sp, anno, how=mean, verbose=FALSE) { return(NULL) } }) + names(merged) <- annoval do.call(rbind, merged) } \ No newline at end of file diff --git a/R/doGSEA.R b/R/doGSEA.R index 89d1f4c..60ed01a 100644 --- a/R/doGSEA.R +++ b/R/doGSEA.R @@ -1,6 +1,7 @@ #' doGSEA #' @export -doGSEA <- function(stana, candSp=NULL, cl=NULL, eps=1e-2, how=mean) { +doGSEA <- function(stana, candSp=NULL, cl=NULL, eps=1e-2, how=mean, + zeroPerc=0) { if (is.null(candSp)) {candSp <- stana@ids[1]} if (is.null(cl)) {cl <- stana@cl} if (length(cl)!=2) {stop("Only the two group is supported")} @@ -22,6 +23,11 @@ doGSEA <- function(stana, candSp=NULL, cl=NULL, eps=1e-2, how=mean) { ko_df_filt <- stana@kos[[candSp]] } + ## If filter based on number of zero per KOs + ko_df_filt <- data.frame(ko_df_filt[!rowSums(ko_df_filt==0)>dim(ko_df_filt)[2] * zeroPerc, ]) |> + `colnames<-`(colnames(ko_df_filt)) + + ## eps values ko_df_filt <- ko_df_filt + 1e-2 ko_sum <- log2(apply(ko_df_filt[,intersect(inSample, aa)],1,mean) / @@ -32,7 +38,7 @@ doGSEA <- function(stana, candSp=NULL, cl=NULL, eps=1e-2, how=mean) { ## KO to PATHWAY mapping url <- "https://rest.kegg.jp/link/ko/pathway" kopgsea <- data.frame(data.table::fread(url, header = FALSE, sep = "\t")) - kopgsea <- kopgsea |> filter(startsWith(V1, "path:ko")) + kopgsea <- kopgsea |> dplyr::filter(startsWith(V1, "path:ko")) kopgsea$V1 <- kopgsea$V1 |> strsplit(":") |> vapply("[",2,FUN.VALUE="a") enr <- clusterProfiler::GSEA(ko_sum, @@ -44,6 +50,7 @@ doGSEA <- function(stana, candSp=NULL, cl=NULL, eps=1e-2, how=mean) { #' @export calcKO <- function(stana, candSp=NULL, how=mean) { if (is.null(candSp)) {candSp <- stana@ids[1]} + if (is.null(stana@eggNOG[[candSp]])) {stop("Please provide list of path to annotation file by `setAnnotation` function.")} ko_df_filt <- summariseAbundance(stana, sp = candSp, checkEGGNOG(annot_file=stana@eggNOG[[candSp]], "KEGG_ko"), how=how) @@ -53,7 +60,7 @@ calcKO <- function(stana, candSp=NULL, how=mean) { #' reverse_annot #' @export -reverse_annot <- function(stana, candSp, candidate, col="KEGG_ko") { +reverseAnnot <- function(stana, candSp, candidate, col="KEGG_ko") { annot <- checkEGGNOG(annot_file=stana@eggNOG[[candSp]], col) annot %>% filter(.data[["value"]] %in% candidate) %>% pull(ID) %>% unique() @@ -62,7 +69,7 @@ reverse_annot <- function(stana, candSp, candidate, col="KEGG_ko") { #' obtain_positions #' @export -obtain_positions <- function(stana, candSp, geneID) { +obtainPositions <- function(stana, candSp, geneID) { tmp <- stana@snpsInfo[[candSp]] tmp[tmp$gene_id %in% geneID, ] %>% row.names() } \ No newline at end of file diff --git a/R/loadMIDAS2Refine.R b/R/loadMIDAS2Refine.R index e7973ac..2664d95 100644 --- a/R/loadMIDAS2Refine.R +++ b/R/loadMIDAS2Refine.R @@ -67,39 +67,10 @@ loadMIDAS2 <- function(midas_merge_dir, stana@cl <- cl } - grnm <- NULL - for (sn in snpsSummary$sample_name) { - tmpgr <- NULL - for (i in names(cl)) { - if (sn %in% cl[[i]]) { - tmpgr <- c(tmpgr, i) - } - } - if (length(tmpgr)>1) { - stop("Sample in multiple groups") - } else { - grnm <- c(grnm, tmpgr) - } - } - - snpsSummary$group <- grnm - - grnm <- NULL - for (sn in genesSummary$sample_name) { - tmpgr <- NULL - for (i in names(cl)) { - if (sn %in% cl[[i]]) { - tmpgr <- c(tmpgr, i) - } - } - if (length(tmpgr)>1) { - stop("Sample in multiple groups") - } else { - grnm <- c(grnm, tmpgr) - } - } - - genesSummary$group <- grnm + changer <- listToNV(cl) + ## error in duplicate + snpsSummary$group <- changer[snpsSummary$sample_name] + genesSummary$group <- changer[genesSummary$sample_name] snpRet <- snpsSummary |> dplyr::group_by(species_id) |> dplyr::count(group) geneRet <- genesSummary |> dplyr::group_by(species_id) |> diff --git a/R/utils.R b/R/utils.R index 1908a17..e04d9c1 100644 --- a/R/utils.R +++ b/R/utils.R @@ -48,3 +48,15 @@ returnGenes <- function(stana, species, snvs) { stop("No genes mapped") } } + + +listToNV <- function(l) { + nm <- NULL + val <- NULL + for (i in names(l)) { + val <- c(val, l[[i]]) + nm <- c(nm, rep(i, length(l[[i]]))) + } + names(nm) <- val + return(nm) +} diff --git a/man/doGSEA.Rd b/man/doGSEA.Rd index 15fe46d..147b769 100644 --- a/man/doGSEA.Rd +++ b/man/doGSEA.Rd @@ -4,7 +4,7 @@ \alias{doGSEA} \title{doGSEA} \usage{ -doGSEA(stana, candSp = NULL, cl = NULL, eps = 0.01, how = mean) +doGSEA(stana, candSp = NULL, cl = NULL, eps = 0.01, how = mean, zeroPerc = 0) } \description{ doGSEA diff --git a/man/obtain_positions.Rd b/man/obtainPositions.Rd similarity index 63% rename from man/obtain_positions.Rd rename to man/obtainPositions.Rd index eba7885..6e9427c 100644 --- a/man/obtain_positions.Rd +++ b/man/obtainPositions.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/doGSEA.R -\name{obtain_positions} -\alias{obtain_positions} +\name{obtainPositions} +\alias{obtainPositions} \title{obtain_positions} \usage{ -obtain_positions(stana, candSp, geneID) +obtainPositions(stana, candSp, geneID) } \description{ obtain_positions diff --git a/man/reverse_annot.Rd b/man/reverseAnnot.Rd similarity index 59% rename from man/reverse_annot.Rd rename to man/reverseAnnot.Rd index 7045603..878cfa1 100644 --- a/man/reverse_annot.Rd +++ b/man/reverseAnnot.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/doGSEA.R -\name{reverse_annot} -\alias{reverse_annot} +\name{reverseAnnot} +\alias{reverseAnnot} \title{reverse_annot} \usage{ -reverse_annot(stana, candSp, candidate, col = "KEGG_ko") +reverseAnnot(stana, candSp, candidate, col = "KEGG_ko") } \description{ reverse_annot