H1975_OCT23

1 Functions

2 Data

genesets_cp <-  msigdb_download("Homo sapiens",category="C2",subcategory = "CP:BIOCARTA ") %>% append( msigdb_download("Homo sapiens",category="C2",subcategory = "CP:KEGG"))

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'append': unknown subcategory

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

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = round(H1975Oct23),
                              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] 9065

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)
dds_H1975_OCT23 = dds

dds = dds_H1975_OCT23

5 Top variable genes heatmap

genes <- head(order(rowVars(assay(dds)), decreasing = TRUE), 1000)

mat <- H1975Oct23[ genes, ]
mat <- t(scale(t(mat)))
anno <- as.data.frame(mat)

library(ComplexHeatmap) 
library(ggplot2) 
Heatmap(mat, cluster_rows = T, cluster_columns = F, column_labels = colnames(anno), name = "fpkm Z-score",row_names_gp = gpar(fontsize = 0)) 

6 DEG FC

cpVSop <- results(dds,contrast = c("condition","comboPersistors","osiPersistors"))  %>% as.data.frame()
roxaVSctrl <- results(dds,contrast = c("condition","roxa","ctrl"))  %>% as.data.frame()
osiVSctrl <- results(dds,contrast = c("condition","osi","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 (10.89% 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 (19.69% 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 = "roxa VS ctrl")

roxa VS ctrl

osiVSctrl = osiVSctrl[order(osiVSctrl$log2FoldChange, osiVSctrl$padj,decreasing = T),] #order by FC, ties bt padj
ranked_vec = osiVSctrl[,"log2FoldChange"]%>% setNames(rownames(osiVSctrl)) %>% 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 (12.4% 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 osi")
plt2 = plt$dn+ aes(size=abs(nes))+ggtitle("up in ctrl")
print_tab(plt1+plt2,title = "osi VS ctrl")

osi VS ctrl

NA

7 CC upregulted in cp VS op

osiVSctrl_genes <- results(dds,contrast = c("condition","osi","ctrl"))  %>% as.data.frame() %>% filter(log2FoldChange>0 & padj < 0.05) %>% rownames()
cpVSctrl_genes <- results(dds,contrast = c("condition","comboPersistors","ctrl"))  %>% as.data.frame() %>% filter(log2FoldChange>0 & padj < 0.05) %>% rownames()
roxaVSctrl_genes = roxaVSctrl %>% filter(log2FoldChange>0 & padj < 0.05)%>% rownames()

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

9 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

H1975_up_genes = up_genes
H1975_down_genes = down_genes

hyp_obj <- hypeR(up_genes, genesets, test = "hypergeometric", fdr=1, plotting=F,background = rownames(H1975Oct23))
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(H1975Oct23))
plt2 = hyp_dots(hyp_obj,title = "down in comboVSosi but not in roxaVSctrl")

plt1 + plt2

print_tab(data.frame(up_genes[up_genes %in% genesets$HALLMARK_E2F_TARGETS]),title = "up genes in E2F")

up genes in E2F

print_tab(data.frame(up_genes[up_genes %in% genesets$HALLMARK_HYPOXIA]),title = "up genes in Hypoxia")

up genes in Hypoxia

NA

10 Expression heatmap

# select the 50 most differentially expressed genes 
genes <- c("DUSP6","MKI67")
mat <- H1975Oct23[ genes, ]
mat <- t(scale(t(mat)))
anno <- as.data.frame(mat)

library(ComplexHeatmap) 
library(ggplot2) 
p = Heatmap(mat, cluster_rows = F, cluster_columns = F, column_labels = colnames(anno), name = "fpkm Z-score") 
print_tab(plt = p,title = "markers")

markers

genes <- hif_targets
mat <- H1975Oct23[genes, ] %>% filter(rowSums(across(where(is.numeric)))!=0)
mat <- t(scale(t(mat)))
anno <- as.data.frame(mat)

 
p = Heatmap(mat, cluster_rows = T, cluster_columns = F, column_labels = colnames(anno), name = "fpkm Z-score",column_title = "HIF targets",row_names_gp = gpar(fontsize = 8))

print_tab(plt = p,title = "HIF targets")

HIF targets

genes <- genesets$HALLMARK_HYPOXIA
mat <- H1975Oct23[genes, ] %>% filter(rowSums(across(where(is.numeric)))!=0)
mat <- t(scale(t(mat)))
anno <- as.data.frame(mat)


p = Heatmap(mat, cluster_rows = T, cluster_columns = F, column_labels = colnames(anno), name = "fpkm Z-score",column_title = "HALLMARK_HYPOXIA",row_names_gp =gpar(fontsize = 0)) 

print_tab(plt = p,title = "HALLMARK_HYPOXIA")

HALLMARK_HYPOXIA

NA

11 G2M heatmap

print_tab(plt = p,title = "HALLMARK_G2M_CHECKPOINT")

HALLMARK_G2M_CHECKPOINT

NA

12 Enrichment analysis

chosen_genes = split %>% dplyr::filter(cluster == 1) %>% rownames() #take relevant genes
hyp_obj <- hypeR(chosen_genes, genesets_wikipathways, test = "hypergeometric", fdr=1, plotting=F,background = rownames(H1975Oct23))
p1 = hyp_dots(hyp_obj,size_by = "none",title = "cluster 1")
chosen_genes = split %>% dplyr::filter(cluster == 4) %>% rownames() #take relevant genes
hyp_obj <- hypeR(chosen_genes, genesets_wikipathways, test = "hypergeometric", fdr=1, plotting=F,background = rownames(H1975Oct23))
p2 = hyp_dots(hyp_obj,size_by = "none",title = "cluster 4")

wiki_p = (p1/p2)

print_tab(plt = wiki_p,title = "wikipathways")

wikipathways

chosen_genes = split %>% dplyr::filter(cluster == 1) %>% rownames() #take relevant genes
hyp_obj <- hypeR(chosen_genes, genesets_biocarta, test = "hypergeometric", fdr=1, plotting=F,background = rownames(H1975Oct23))
p1 = hyp_dots(hyp_obj,size_by = "none",title = "cluster 1")
chosen_genes = split %>% dplyr::filter(cluster == 4) %>% rownames() #take relevant genes
hyp_obj <- hypeR(chosen_genes, genesets_biocarta, test = "hypergeometric", fdr=1, plotting=F,background = rownames(H1975Oct23))
p2 = hyp_dots(hyp_obj,size_by = "none",title = "cluster 4")

wiki_p = (p1/p2)

print_tab(plt = wiki_p,title = "biocarta")

biocarta

chosen_genes = split %>% dplyr::filter(cluster == 1) %>% rownames() #take relevant genes
hyp_obj <- hypeR(chosen_genes, genesets_pid, test = "hypergeometric", fdr=1, plotting=F,background = rownames(H1975Oct23))
p1 = hyp_dots(hyp_obj,size_by = "none",title = "cluster 1")
chosen_genes = split %>% dplyr::filter(cluster == 4) %>% rownames() #take relevant genes
hyp_obj <- hypeR(chosen_genes, genesets_pid, test = "hypergeometric", fdr=1, plotting=F,background = rownames(H1975Oct23))
p2 = hyp_dots(hyp_obj,size_by = "none",title = "cluster 4")

wiki_p = (p1/p2)

print_tab(plt = wiki_p,title = "pid")

pid

chosen_genes = split %>% dplyr::filter(cluster == 1) %>% rownames() #take relevant genes
hyp_obj <- hypeR(chosen_genes, genesets_REACTOME, test = "hypergeometric", fdr=1, plotting=F,background = rownames(H1975Oct23))
p1 = hyp_dots(hyp_obj,size_by = "none",title = "cluster 1")
chosen_genes = split %>% dplyr::filter(cluster == 4) %>% rownames() #take relevant genes
hyp_obj <- hypeR(chosen_genes, genesets_REACTOME, test = "hypergeometric", fdr=1, plotting=F,background = rownames(H1975Oct23))
p2 = hyp_dots(hyp_obj,size_by = "none",title = "cluster 4")

wiki_p = (p1/p2)

print_tab(plt = wiki_p,title = "REACTOME")

REACTOME

NA

for (chosen_clusters in 1:length(unique(split$cluster))) {
chosen_genes = split %>% dplyr::filter(cluster == chosen_clusters) %>% rownames() #take relevant genes
hyp_obj <- hypeR(chosen_genes, genesets_wikipathways, test = "hypergeometric", fdr=1, plotting=F,background = rownames(H1975Oct23))
print(hyp_dots(hyp_obj,size_by = "none"))
  
  }

13 Distance plot

vsd <- vst(dds, blind=FALSE)
sampleDists <- dist(t(assay(vsd)))
library("RColorBrewer")
sampleDistMatrix <- as.matrix(sampleDists)
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows=sampleDists,
         clustering_distance_cols=sampleDists,
         col=colors)