xeno_bulk_TPM_noMouseReads

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 = "/sci/labs/yotamd/lab_share/lung_sc/bulk/RoxaOsi/Oct23/star/mm39unmapped_noMTGLKI_tpm.txt",
  sep = "\t",
  header = TRUE
)
rownames(xeno_bulk) = make.unique(xeno_bulk[,"gene_name",drop=T])
xeno_bulk = xeno_bulk[,8:19]
names (xeno_bulk) = gsub(pattern = "_mm.*$",replacement = "",x = names (xeno_bulk))
xeno_bulk = xeno_bulk[,c("M2_ctrl", "M21_ctrl", "M36_ctrl", "M14_osi", "M33_osi", "M14_2_osi", 
"M20_roxa", "M30_roxa", "M31_roxa", "M3_combo", "M26_combo", 
"M32_combo")]

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] 60580

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

[1] 21163

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

6 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 = "TPM 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 = "TPM Z-score",column_title = "HIF targets",row_names_gp = gpar(fontsize = 8)) 

genes <- genesets$HALLMARK_G2M_CHECKPOINT
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 = "TPM Z-score",column_title = "G2M targets",row_names_gp = gpar(fontsize = 0)) 

genes <- genesets$HALLMARK_HYPOXIA
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 = "TPM Z-score",column_title = "Hallmark Hypoxia",row_names_gp = gpar(fontsize = 0)) 

NA
NA

7 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)
roxa_vs_ctrl = results(dds,contrast = c("condition","roxa","ctrl"))  %>% as.data.frame()
combo_vs_ctrl = results(dds,contrast = c("condition","combo","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

NA

library(ggvenn)
deg = list(combo_vs_ctrl = combo_vs_ctrl%>% filter(padj <0.05) %>% rownames(),osi_vs_ctrl = osi_vs_ctrl%>% filter(padj <0.05) %>% rownames())
ggvenn(deg, c("combo_vs_ctrl", "osi_vs_ctrl"))            # draw two-set venn

deg = list(combo_vs_roxa = combo_vs_roxa%>% filter(padj <0.05) %>% rownames(),osi_vs_roxa = osi_vs_roxa%>% filter(padj <0.05) %>% rownames())
ggvenn(deg, c("combo_vs_roxa", "osi_vs_roxa"))            # draw two-set venn

deg = list(osi_vs_ctrl = osi_vs_ctrl%>% filter(padj <0.05) %>% rownames(),osi_vs_roxa = osi_vs_roxa%>% filter(padj <0.05) %>% rownames())
ggvenn(deg, c("osi_vs_ctrl", "osi_vs_roxa"))            # draw two-set venn