GSEA algorithms for the function getPatternHallmarks about cogaps

[dimalvovs]

@LiuCanidk thanks for your questions. Please note that getPatternHallmarks is discontinued in the current version (see #92) and replaced with getPatternGeneSet method to aid stability, eliminate dependency on possibly outdated packages, and give users more freedom regarding the choice of hallmarks. Do you want to check that out?

[LiuCanidk]

@dimalvovs of course! Should I remove the old package and install the developmental version on the github? Then I would check out the new function and provide feedback.

[dimalvovs]

yes, that's the way

[LiuCanidk]

@dimalvovs , I try to use the new function. But first I need an order ranked gene list. So I plan to obtain this from the result of patternMarkers. But I found strange results that the rank and the score of genes has a positive correlation. As stated in the tutorial, lower rank represent higher association with specific patterns and higher score represent higher association. They should be negatively correlated.

code

#already with a cogaps result
pm <- patternMarkers(cogapsresult, threshold="cut")
save(pm, file='05.NMF_expression_program/NMF_CoGAPS_marker.rda')
#check marker rank
head(pm$PatternMarkers[[1]])
head(pm$PatternMarkerRanks[pm$PatternMarkers[[1]],])
#check marker score
head(pm$PatternMarkerScores[pm$PatternMarkers[[1]],])
#check the correlation between rank and score
cor(pm$PatternMarkerRanks[,1], pm$PatternMarkerScores[,1])

result:

Also, I did find that the rank of marker genes in pm$PatternMarkerRanks were not corresponding to the rank in the pm$PatternMarkers one by one. As you can see in the screenshot, gene PNMA8B is not NO.2 but NO.3.

So this made me confusing about how to choose the gene rank when perform the enrichement test in the getPatternGeneSet function.

[LiuCanidk]

@dimalvovs sorry, I may misunderstand the use of getPatternGeneSet, the geneset parameter should be some pathway gene list, such as Hallmark, KEGG GO etc. This was a little bit misleading.... And also, I recommend to provide additional choice for gmt file input, which can be easily downloaded from MsigDB.

[LiuCanidk]

@dimalvovs the new function is great. Only note that avoid the misleading use of geneset parameter as I mentioned above. Maybe providing the parameter choice of KEGG GO or Hallmark pathway list is more straightforward and user-friendly. In addition, for GSEA, the GSEA curve plot is necessary for specific pathway visualization (see below). Considering the usage of clusterProfiler function gseaplot and database R package msigdbr may help.

GSEA curve plot

Without provided list, my code was as follows:

#Enrichment analysis
library(clusterProfiler)
#readin the hallmark pathway gmt file
hallmark=read.gmt('F:/MSigDB/h.all.v2023.2.Hs.symbols.gmt')
#convert this dataframe to a gene list with names corresponding to each pathway
hallmark.list=list()
for(i in unique(hallmark$term)){
  term=tolower(gsub('HALLMARK_','',i))
  hallmark.list[[term]]=hallmark$gene[hallmark$term==i]
}
#perform the GSEA
ORA.hallmark=getPatternGeneSet(cogapsresult,
                      hallmark.list,
                      'enrichment')
#draw plot based on padj & NES, top10 pathways
for(i in 1:7){
  ORA.hallmark[[i]]=ORA.hallmark[[i]][order(ORA.hallmark[[i]]$padj),][1:10,]
  ORA.hallmark[[i]]$pathway=factor(ORA.hallmark[[i]]$pathway,
                                   levels = ORA.hallmark[[i]]$pathway)
  assign(paste0('p',i),
         ggplot(ORA.hallmark[[i]], aes(-log10(padj),pathway,fill=NES))+
           geom_col()+theme_bw(base_size = 12)+
           xlab(expression(-log[10])~padj)+ylab(NULL)+
           geom_vline(xintercept = -log10(0.0005), lty=2)+
           scale_x_continuous(expand = c(0,0))+
           scale_fill_continuous(low='lightblue',high='indianred')+
           geom_text(data = ORA.hallmark[[i]],
                     aes(x = 0.05, 
                         y = pathway,
                         label = pathway),
                     size = 4,
                     hjust = 0)+
           theme(axis.text.y = element_blank()))
}
#arrange and save PDF
library(ggpubr)
ggsave('05.NMF_expression_program/GSEA_hallmark.pdf',
       ggarrange(p1,p2,p3,p4,p5,p6,p7,
                    nrow=3, ncol=3, common.legend = T),
       width = 32, height = 24, units = 'cm')

plot result

plotPatternGeneSet result

Obviously, maybe providing a function that has better visualization than plotPatternGeneSet is highly recommended.