H1975_JUN23

1 Functions

2 Data

names (osiRoxaJun23)

 [1] "H1975ctrl1"             "H1975ctrl2"             "H1975ctrl3"             "H1975roxa1"             "H1975roxa2"            
 [6] "H1975roxa3"             "H1975osi1"              "H1975osi2"              "H1975osi3"              "H1975osiPersistors1"   
[11] "H1975osiPersistors2"    "H1975osiPersistors3"    "H1975comboPersistors1"  "H1975comboPersistors2"  "H1975comboPersistors3" 
[16] "HCC827ctrl1"            "HCC827ctrl2"            "HCC827ctrl3"            "HCC827roxa1"            "HCC827roxa2"           
[21] "HCC827roxa3"            "HCC827osi1"             "HCC827osi2"             "HCC827osi3"             "HCC827osiPersistors1"  
[26] "HCC827osiPersistors2"   "HCC827osiPersistors3"   "HCC827comboPersistors1" "HCC827comboPersistors2" "HCC827comboPersistors3"

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

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

3 PCA

nrow(dds)

[1] 27827

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

[1] 8340

vst = vst(dds1, 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)

5 DEG FC

cpVSop <- results(dds,contrast = c("condition","comboPersistors","osiPersistors"))  %>% 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)

cpVSop = cpVSop[order(cpVSop$log2FoldChange, cpVSop$padj,decreasing = T),] #order by FC, ties bt padj
ranked_vec = cpVSop[,"log2FoldChange"]%>% setNames(rownames(cpVSop)) %>% na.omit() # make named vector

hyp_obj <- hypeR_fgsea(ranked_vec, genesets, up_only = F)

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

plt = hyp_dots(hyp_obj,merge = F)
plt1 = plt$up+ aes(size=nes)+ggtitle("up in comboPersistor")
plt2 = plt$dn+ aes(size=abs(nes))+ggtitle("up in osiPersistors")
print_tab(plt1+plt2,title = "cpVSop")

cpVSop

roxaVSctrl = roxaVSctrl[order(roxaVSctrl$log2FoldChange, roxaVSctrl$padj,decreasing = T),] #order by FC, ties bt padj
ranked_vec = roxaVSctrl[,"log2FoldChange"]%>% setNames(rownames(roxaVSctrl)) %>% na.omit()  # make named vector

hyp_obj <- hypeR_fgsea(ranked_vec, genesets, up_only = F)

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

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")

cpVSop

NA

6 DEG shrinked FC

dds$condition = relevel(dds$condition, ref = "osiPersistors")
dds <- nbinomWaldTest(dds)
cpVSop <- lfcShrink(dds,coef = "condition_comboPersistors_vs_osiPersistors")  %>% 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)

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

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")

cpVSop

ranked_vec = diff_genes[, 2] %>% setNames(rownames(diff_genes)) %>% sort(decreasing = TRUE)
hyp_obj <- hypeR_fgsea(ranked_vec, genesets, up_only = F)

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

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")

cpVSop

NA

7 DEG in comboVSosi but not in roxaVSctrl

cpVSop <- results(dds,contrast = c("condition","comboPersistors","osiPersistors"))  %>% 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")

up

print_tab(down_genes_df,title = "down")

down

NA

hyp_obj <- hypeR(up_genes, genesets, test = "hypergeometric", fdr=1, plotting=F,background = rownames(osiRoxaJun23))
plt1 = hyp_dots(hyp_obj,title = "up in comboVSosi but not in roxaVSctrl")
 

hyp_obj <- hypeR(down_genes, genesets, test = "hypergeometric", fdr=1, plotting=F,background = rownames(osiRoxaJun23))
plt2 = hyp_dots(hyp_obj,title = "down in comboVSosi but not in roxaVSctrl")

plt1 + plt2