HCC827_JUN23

1 Functions

2 Data

HCC827Jun23 = read.table(
  file = "./Data/osiRoxa_bulk/Jun23/gene_fpkm.xls",
  sep = "\t",
  header = TRUE
)
rownames(HCC827Jun23) = make.unique(HCC827Jun23[,"gene_name",drop=T]) #set genes names
HCC827Jun23 = HCC827Jun23[,2:31]
names (HCC827Jun23) = gsub(x = names(HCC827Jun23),pattern = "op",replacement = "osiPersistors")%>% gsub(pattern = "cp",replacement = "comboPersistors") 
HCC827Jun23 = HCC827Jun23[,16:ncol(HCC827Jun23)] #take only HCC827

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

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

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

[1] 8494

vst = vst(dds1, blind=FALSE)

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

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

#DESeq

dds <- DESeq(dds)

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

dds_HCC827_JUN23 = dds

dds = dds_HCC827_JUN23

4 Top variable genes heatmap

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

mat <- HCC827Jun23[ 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)) 

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 (11.24% 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 (12.19% 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 (23.49% 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 (23.49% 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(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_GLYCOLYSIS]),title = "up genes in GLYCOLYSIS")

up genes in GLYCOLYSIS

NA

8 Intersection HCC827 and H1975 DEG

HCC827_down_genes

[1] "IFI6"      "CDKN2B"    "SELENBP1"  "RENBP"     "SLC6A9"    "LINC02582"

intersect(HCC827_up_genes,H1975_up_genes)

[1] "HMGB2"

intersect(HCC827_down_genes,H1975_down_genes)

character(0)

9 Expression heatmap

# select the 50 most differentially expressed genes 
# genes <- c("DUSP6")
# mat <- HCC827Jun23[ 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")

genes <- hif_targets
mat <- HCC827Jun23[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 <- HCC827Jun23[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

genes <- genesets$HALLMARK_HYPOXIA
mat <- HCC827Jun23[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

10 G2M heatmap

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

split = data.frame( cluster = cutree(hclust(dist(mat)), k = 3))
ha = rowAnnotation(cluster = as.character(split$cluster),col = list(cluster = structure(rainbow(n = 3), names=unique(split$cluster))))
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),right_annotation = ha)

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

HALLMARK_G2M_CHECKPOINT

NA

11 Enrichment analysis

chosen_genes = split %>% dplyr::filter(cluster == 2) %>% 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 2")
chosen_genes = split %>% dplyr::filter(cluster == 3) %>% 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 3")

wiki_p = (p1/p2)

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

wikipathways

chosen_genes = split %>% dplyr::filter(cluster == 2) %>% 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 2")
chosen_genes = split %>% dplyr::filter(cluster == 3) %>% 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 3")

wiki_p = (p1/p2)

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

biocarta

chosen_genes = split %>% dplyr::filter(cluster == 2) %>% 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 2")
chosen_genes = split %>% dplyr::filter(cluster == 3) %>% 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 3")

wiki_p = (p1/p2)

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

pid

chosen_genes = split %>% dplyr::filter(cluster == 2) %>% 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 2")
chosen_genes = split %>% dplyr::filter(cluster == 3) %>% 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 3")

wiki_p = (p1/p2)

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

REACTOME

NA

12 Distance plot

vsd <- vst(dds, blind=FALSE)

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

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)