xeno_bulk

1 Functions

2 Data

genesets <- msigdb_download("Homo sapiens",category="H") %>% append( msigdb_download("Homo sapiens",category="C2",subcategory = "CP:KEGG"))
xeno_bulk = read.table(
  file = "./Data/osiRoxa_bulk/Oct23/gene_fpkm.xls",
  sep = "\t",
  header = TRUE
)
rownames(xeno_bulk) = make.unique(xeno_bulk[,"gene_name",drop=T])
xeno_bulk = xeno_bulk[,26:37]

cell.labels = names(xeno_bulk)
condition = str_extract(cell.labels, "osi|combo|ctrl|roxa")
metadata = data.frame(condition = condition, row.names = colnames(xeno_bulk))

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = round(xeno_bulk),
                              colData = metadata,
                              design = ~condition)

converting counts to integer mode
Warning in DESeqDataSet(se, design = design, ignoreRank) :
  some variables in design formula are characters, converting to factors

3 PCA

nrow(dds)

[1] 32780

dds1 <- dds[ rowSums(counts(dds)) >= 3, ]
nrow(dds1)

[1] 9259

vst = vst(dds, blind=FALSE)

library("ggfortify")
PCAdata <- prcomp(t(assay(vst)))
autoplot(PCAdata, data = metadata,colour = "condition",label = FALSE, main="PCA") # Show dots

4 DESeq

dds <- DESeq(dds)

estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

5 DEG

cpVSop <- results(dds,contrast = c("condition","combo","osi"))  %>% as.data.frame()
roxaVSctrl <- results(dds,contrast = c("condition","roxa","ctrl"))  %>% as.data.frame()
diff_genes = data.frame(row.names = rownames(cpVSop), cpVSop_FC = cpVSop$log2FoldChange,roxaVSctrl_FC = roxaVSctrl$log2FoldChange,  cpVSop_padj = cpVSop$padj)

library(DT)
combo_vs_ctrl = results(dds,contrast = c("condition","combo","ctrl"))  %>% as.data.frame()
roxa_vs_ctrl = results(dds,contrast = c("condition","roxa","ctrl"))  %>% as.data.frame()
osi_vs_ctrl = results(dds,contrast = c("condition","osi","ctrl"))  %>% as.data.frame()

combo_vs_roxa = results(dds,contrast = c("condition","combo","roxa"))  %>% as.data.frame()
osi_vs_roxa = results(dds,contrast = c("condition","osi","roxa"))  %>% as.data.frame()
combo_vs_osi <- results(dds,contrast = c("condition","combo","osi"))  %>% as.data.frame()

dt1 = datatable(roxa_vs_ctrl %>% filter(padj <0.05),caption = "significant roxa VS ctrl")
dt2= datatable(combo_vs_ctrl %>% filter(padj <0.05),caption = "significant combo vs ctrl")
dt3 = datatable(osi_vs_ctrl%>% filter(padj <0.05),caption = "significant osi_vs_ctrl")
dt4 =datatable(combo_vs_roxa %>% filter(padj <0.05),caption = "significant combo_vs_roxa")
dt5 =datatable(osi_vs_roxa %>% filter(padj <0.05),caption = "significant osi_vs_roxa")
dt6 =datatable(combo_vs_osi %>% filter(padj <0.05),caption = "significant combo VS osi")
print_tab(plt = dt1,title = "significant roxa VS ctrl")

##   significant roxa VS ctrl {.unnumbered }  

print_tab(plt = dt2,title = "significant combo vs ctrl")

##   significant combo vs ctrl {.unnumbered }  

print_tab(plt = dt3,title = "significant osi VS ctrl")

##   significant osi VS ctrl {.unnumbered }  

print_tab(plt = dt4,title = "significant combo VS roxa")

##   significant combo VS roxa {.unnumbered }  

print_tab(plt = dt5,title = "significant osi VS roxa")

##   significant osi VS roxa {.unnumbered }  

print_tab(plt = dt6,title = "significant combo VS osi")

##   significant combo VS osi {.unnumbered }  

NA

6 GSEA

plot_hyper <- function(genes_df,ident.1,ident.2) {
  genes_df = cpVSop[order(genes_df$log2FoldChange, genes_df$padj, decreasing = T), ] #order by FC, ties bt padj
  ranked_vec = genes_df[, "log2FoldChange"] %>% setNames(rownames(genes_df)) %>% na.omit() # make named vector
  
  hyp_obj <- hypeR_fgsea(ranked_vec, genesets, up_only = F)
  plt = hyp_dots(hyp_obj, merge = F)
  plt1 = plt$up + aes(size = nes) + ggtitle(paste("up in" ,ident.1))
  plt2 = plt$dn + aes(size = abs(nes)) + ggtitle(paste("up in" ,ident.2))
  print_tab(plt1 + plt2, title = "title")
}

plot_hyper(combo_vs_ctrl,ident.1 = "combo",ident.2 = "ctrl")

Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (15.25% of the list). The order of those tied genes will be arbitrary, which may produce unexpected results. ## title {.unnumbered }

plot_hyper(osi_vs_ctrl,ident.1 = "osi",ident.2 = "ctrl")

Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (15.25% of the list). The order of those tied genes will be arbitrary, which may produce unexpected results. ## title {.unnumbered }

plot_hyper(combo_vs_roxa,ident.1 = "combo",ident.2 = "roxa")

Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (15.25% of the list). The order of those tied genes will be arbitrary, which may produce unexpected results. ## title {.unnumbered }

plot_hyper(osi_vs_roxa,ident.1 = "osi",ident.2 = "roxa")

Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (15.25% of the list). The order of those tied genes will be arbitrary, which may produce unexpected results. ## title {.unnumbered }

NA

HALLMARK_HYPOXIA genes in combo_vs_ctrl_genes

genesets$HALLMARK_HYPOXIA[genesets$HALLMARK_HYPOXIA %in% combo_vs_ctrl_genes] 

[1] "BGN"    "FOSL2"  "GPC1"   "KLHL24" "NDST1"  "S100A4"

7 DEG shrinked FC

dds$condition = relevel(dds$condition, ref = "osi")
dds <- nbinomWaldTest(dds)
cpVSop <- lfcShrink(dds,coef = "condition_combo_vs_osi")  %>% as.data.frame()

dds$condition = relevel(dds$condition, ref = "ctrl")
dds <- nbinomWaldTest(dds)
roxaVSctrl <- lfcShrink(dds,coef  = "condition_roxa_vs_ctrl")  %>% as.data.frame()


diff_genes = data.frame(row.names = rownames(cpVSop), cpVSop_FC = cpVSop$log2FoldChange,roxaVSctrl_FC = roxaVSctrl$log2FoldChange,  cpVSop_padj = cpVSop$padj)

ranked_vec = diff_genes[, 1] %>% setNames(rownames(diff_genes)) %>% sort(decreasing = TRUE)
hyp_obj <- hypeR_fgsea(ranked_vec, genesets, up_only = F)
plt = hyp_dots(hyp_obj,merge = F)
plt1 = plt$up+ aes(size=nes)+ggtitle("up in comboPersistor") + theme(  axis.text.y = element_text(size=10))
plt2 = plt$dn+ aes(size=abs(nes))+ggtitle("up in osiPersistors") + theme(axis.text.y = element_text(size=10))
print_tab(plt1+plt2,title = "cpVSop")

ranked_vec = diff_genes[, 2] %>% setNames(rownames(diff_genes)) %>% sort(decreasing = TRUE)
hyp_obj <- hypeR_fgsea(ranked_vec, genesets, up_only = F)
plt = hyp_dots(hyp_obj,merge = F)
plt1 = plt$up+ aes(size=nes)+ggtitle("up in roxa")
plt2 = plt$dn+ aes(size=abs(nes))+ggtitle("up in ctrl")
print_tab(plt1+plt2,title = "cpVSop")

8 DEG in comboVSosi but not in roxaVSctrl

cpVSop <- results(dds,contrast = c("condition","combo","osi"))  %>% as.data.frame()
roxaVSctrl <- results(dds,contrast = c("condition","roxa","ctrl"))  %>% as.data.frame()
diff_genes = data.frame(row.names = rownames(cpVSop), cpVSop_FC = 2**cpVSop$log2FoldChange,roxaVSctrl_FC = 2**roxaVSctrl$log2FoldChange,  cpVSop_padj = cpVSop$padj)
up_genes_df =  diff_genes %>% filter(cpVSop_FC > 2 & roxaVSctrl_FC<1.2 & cpVSop_padj<0.05) 
down_genes_df = diff_genes %>% filter(cpVSop_FC < 0.5 & roxaVSctrl_FC>0.8 & cpVSop_padj<0.05)
up_genes = diff_genes %>% filter(cpVSop_FC > 2 & roxaVSctrl_FC<1.2 & cpVSop_padj<0.05) %>% rownames()
down_genes = diff_genes %>% filter(cpVSop_FC < 0.5 & roxaVSctrl_FC>0.8 & cpVSop_padj<0.1)%>% rownames()

print_tab(up_genes_df,title = "up")
print_tab(down_genes_df,title = "down")

cpVSop <- results(dds,contrast = c("condition","roxa","osi"))  %>% as.data.frame() %>% filter(padj<0.05)