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

dds_xeno = dds

dds = dds_xeno

5 Top variable genes heatmap

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

mat <- xeno_bulk[ 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

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

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

significant combo vs ctrl

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

significant osi VS ctrl

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

significant combo VS roxa

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

significant osi VS roxa

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

significant combo VS osi

roxa_vs_ctrl["C4orf3",]

7 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))
 return(plt1+plt2)
}

plot_hyper(combo_vs_ctrl,ident.1 = "combo",ident.2 = "ctrl") %>%  print_tab(title = "combo VS ctrl")

combo VS ctrl

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.

plot_hyper(osi_vs_ctrl,ident.1 = "osi",ident.2 = "ctrl")%>%  print_tab(title = "osi VS ctrl")

osi VS ctrl

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.

plot_hyper(combo_vs_roxa,ident.1 = "combo",ident.2 = "roxa") %>%  print_tab(title = "combo VS roxa")

combo VS roxa

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.

plot_hyper(osi_vs_roxa,ident.1 = "osi",ident.2 = "roxa") %>%  print_tab(title = "osi VS roxa")

osi VS roxa

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.

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"

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

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

10 Expression heatmap

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

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

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

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

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)