signature_compare

1 Functions

source_from_github(repositoy = "DEG_functions",version = "0.2.54")

ℹ SHA-1 hash of file is ee1432006d2446b41e3bbd7e217de81ee047abc7
Welcome to enrichR
Checking connection ... 
Enrichr ... Connection is Live!
FlyEnrichr ... Connection is available!
WormEnrichr ... Connection is available!
YeastEnrichr ... Connection is available!
FishEnrichr ... Connection is available!

2 signature

# find hpv signature genes
hmsc_deg = read.table(file = "./Data_out/HPV_analysis/hpv_deg_df.csv",row.names = 1)
opscc_deg = readRDS(file = "./Data_out/opscc_deg_5Kvargenes.rds")
scc_deg = readRDS(file = "./Data_out/scc_deg_5Kvargenes.rds")


DT::datatable(hmsc_deg, caption = "hmsc_deg")

DT::datatable(opscc_deg, caption = "opscc_deg")

DT::datatable(scc_deg, caption = "scc_deg")

avg_diff_cutoff = 0.1
fdr_cutoff = 0.1

hmsc_top = hmsc_deg %>% filter(avg_diff > avg_diff_cutoff &
                                 fdr < fdr_cutoff) %>% rownames()
opscc_top = opscc_deg %>% filter(avg_diff > avg_diff_cutoff &
                                   fdr < fdr_cutoff) %>% rownames()
scc_top = scc_deg %>% filter(avg_diff > avg_diff_cutoff &
                               fdr < fdr_cutoff) %>% rownames()

print ("intersect genes:")

[1] "intersect genes:"

intersect(intersect(hmsc_top,opscc_top),scc_top)

 [1] "MT2A"   "TUBA1C" "TUBA1B" "UBE2T"  "LGALS1" "NUCB2"  "KRT5"   "MYC"    "KRT18"  "BTG2"  

all_deg = list(hmsc_top = hmsc_top, opscc_top = opscc_top, scc_top = scc_top)

3 Venn

library(ggVennDiagram)
ggVennDiagram(all_deg)

4 HMSC GSEA

5 OPSCC GSEA

ranked_vec = opscc_deg[,"avg_diff"]%>% setNames(rownames(opscc_deg)) %>% na.omit() # make named vector
hyp_obj <-hypeR_fgsea(signature = ranked_vec,genesets = geneIds(genesets_h),up_only = F)

plt = hyp_dots(hyp_obj,merge = F,fdr = 0.2)
plt[[1]] + aes(size=nes)+ 
  guides(
    size = guide_legend(title="NES",reverse=T))

plt[[2]] + aes(size=nes) +scale_size(trans = 'reverse')+
  guides(
    size = guide_legend(title="NES",reverse=F))

Scale for 'size' is already present. Adding another scale for 'size', which will replace the existing scale.

6 SCC GSEA

ranked_vec = scc_deg[,"avg_diff"]%>% setNames(rownames(scc_deg)) %>% na.omit() # make named vector
hyp_obj <-hypeR_fgsea(signature = ranked_vec,genesets = geneIds(genesets_h),up_only = F)

plt = hyp_dots(hyp_obj,merge = F,fdr = 0.2)
plt[[1]] + aes(size=nes)+ 
  guides(
    size = guide_legend(title="NES",reverse=T))

plt[[2]] + aes(size=nes) +scale_size(trans = 'reverse')+
  guides(
    size = guide_legend(title="NES",reverse=F))

Scale for 'size' is already present. Adding another scale for 'size', which will replace the existing scale.

7 OPSCC boxplot

top_genes = hmsc_deg %>% arrange(desc(.$avg_diff)) %>% head(5) %>% rownames()

top_genes_expr = FetchData(object = opscc,vars = c("hpv",top_genes))
df = reshape2::melt(top_genes_expr,value.name = "Expression") %>% dplyr::rename(gene = variable)

library(rstatix)
stat.test <- df %>%
    group_by(gene) %>%
  wilcox_test(Expression ~ hpv) %>%
  mutate(y.position = 5)

stat.test

stat.test <- stat.test %>% 
  add_xy_position(x = "gene", dodge = 0.8)

ggboxplot(
  df,
  x = "gene",
  y = "Expression",
  color = "hpv",
  palette = "jco",
  add = c("mean","boxplot"),trim = T
)+ stat_pvalue_manual(stat.test, label = "p = {p}",remove.bracket = T)

8 SCC boxplot

top_genes = hmsc_deg %>% arrange(desc(.$avg_diff)) %>% head(5) %>% rownames()

top_genes_expr = FetchData(object = scc_myb_patients,vars = c("hpv_positive",top_genes))
df = reshape2::melt(top_genes_expr,value.name = "Expression") %>% dplyr::rename(gene = variable)

library(rstatix)
stat.test <- df %>%
    group_by(gene) %>%
  wilcox_test(Expression ~ hpv_positive) %>%
  mutate(y.position = 5)

stat.test

stat.test <- stat.test %>% 
  add_xy_position(x = "gene", dodge = 0.8)

ggboxplot(
  df,
  x = "gene",
  y = "Expression",
  color = "hpv_positive",
  palette = "jco",
  add = c("mean","boxplot"),trim = T
)+ stat_pvalue_manual(stat.test, label = "p = {p}",remove.bracket = T)