H1975_OCT23

1 Functions

2 Data

genesets <- msigdb_download("Homo sapiens",category="H") %>% append( msigdb_download("Homo sapiens",category="C2",subcategory = "CP:KEGG"))
a = msigdb_download("Homo sapiens",category="C2",subcategory = "CP")
H1975Oct23 = read.table(
  file = "./Data/osiRoxa_bulk/Oct23/gene_fpkm.xls",
  sep = "\t",
  header = TRUE
)
rownames(H1975Oct23) = make.unique(H1975Oct23[,"gene_name",drop=T])
H1975Oct23 = H1975Oct23[,2:16]
names (H1975Oct23) = gsub(x = names(H1975Oct23),pattern = "_C",replacement = "_ctrl")%>% gsub(pattern = "p_OR",replacement = "_comboPersistors") %>% gsub(pattern = "p_O",replacement = "_osiPersistors") %>% gsub(pattern = "_R",replacement = "_roxa")%>% gsub(,pattern = "_O",replacement = "_osi")

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


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

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)

Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (16.5% 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 (16.41% 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

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

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

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

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_G2M_CHECKPOINT
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_G2M_CHECKPOINT",row_names_gp =gpar(fontsize = 0)) 

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

HALLMARK_G2M_CHECKPOINT

NA

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