source("hu_dataprep.R")  ## read table and prepare normalized abundances (NA per row < 24); change table to read either with adj.pval or normal p.val
source("heatmap,pca,dist_functions.R") ## get di(), hm(), pca() functions
library(openxlsx)
library(tidyverse)
library(ggrepel)
x = med_nabu[,!grepl("_SP",colnames(med_nabu))]
z = log_nabu[,!grepl("_SP",colnames(log_nabu))]
y <-"./figures/hu_all_overview.pdf"
a_all<-data.frame(group=c( rep("NEUR_M", 5),
rep("RMG_M", 5),
rep("NEUR_P", 5),
rep("RMG_P", 5)),
row.names = colnames(x))
acol_all <- list(group=c(RMG_M="#6A94D4",
RMG_P="#FFCB73",
NEUR_M = "#224C8A",
NEUR_P = "#A67625"
))
di(z)
pdf(y, height = 10/2.54, width = 10/2.54)
hm(x, dis, dis2)
pca(z)
plot(dis)
dev.off()
library(limma)
l = log2(RMG)
colnames(l) = c(rep("RMG_M",5),rep("RMG_P",5))
fac = factor(colnames(l))
design <- model.matrix(~0+fac)
colnames(design) = c("RMG_M","RMG_P")
fit = lmFit(l, design = design)
contrast.matrix = makeContrasts(RMG_M-RMG_P, levels = design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
deProt = topTable(fit2,number = Inf, adjust="BH")
deProt = deProt[deProt$adj.P.Val<0.05 & !is.na(deProt$adj.P.Val),]
deProt = deProt[order(deProt$logFC, decreasing = T),]
log_limma = data.matrix(log2(RMG[rownames(deProt),]))
med_limma = f(a, log_limma)
log_limma[is.na(log_limma)] = 0
rownames(deProt) = all[rownames(deProt),1]
## glia specific genes: adjs.p.val < 0.05 & detected in at least 3 glia samples
RMG_M = all[(all$adj.p.val_RvN_M < 0.05 | is.na(all$adj.p.val_RvN_M))  & rowSums(!is.na(all[,7:11]))>=3,] # RMG vs. Neuron specific in Macula
RMG_P = all[(all$adj.p.val_RvN_P < 0.05 | is.na(all$adj.p.val_RvN_P)) & rowSums(!is.na(all[,17:21]))>=3,] # RMG vs. Neuron specific in periphery
RMG_SP = all[(all$adj.p.val_RvN_SP < 0.05 | is.na(all$adj.p.val_RvN_SP)) & rowSums(!is.na(all[,26:30]))>=3,] #RMG vs. Neuron specific in superperiphery
RMG_M[is.na(RMG_M)] = 0
RMG_P[is.na(RMG_P)] = 0
RMG_SP[is.na(RMG_SP)] = 0
# glia vs. neuron ratio >1
RMG_M = rownames(RMG_M[rowMeans(RMG_M[,7:11], na.rm = T) > rowMeans(RMG_M[,2:6], na.rm = T),])
RMG_P = rownames(RMG_P[rowMeans(RMG_P[,17:21], na.rm = T) > rowMeans(RMG_P[,12:16], na.rm = T) ,])
RMG_SP = rownames(RMG_SP[rowMeans(RMG_SP[,26:30], na.rm = T) > rowMeans(RMG_SP[,22:25], na.rm = T) ,])
### glia-spec, only glia, only macula, periphery
RMG = all[unique(c(RMG_M,RMG_P)),2:30]
RMG = RMG[,grep("RMG", colnames(RMG))]
RMG = RMG[,1:10]
log_RMG = data.matrix(log2(RMG))
med_RMG = f(a, log_RMG)
log_RMG[is.na(log_RMG)] = 0
#write.xlsx(data.frame(RMG,all[rownames(RMG),1]), "rmg.xlsx")
library(limma)
l = log2(RMG)
colnames(l) = c(rep("RMG_M",5),rep("RMG_P",5))
fac = factor(colnames(l))
design <- model.matrix(~0+fac)
colnames(design) = c("RMG_M","RMG_P")
fit = lmFit(l, design = design)
contrast.matrix = makeContrasts(RMG_M-RMG_P, levels = design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
deProt = topTable(fit2,number = Inf, adjust="BH")
deProt = deProt[deProt$adj.P.Val<0.05 & !is.na(deProt$adj.P.Val),]
deProt = deProt[order(deProt$logFC, decreasing = T),]
log_limma = data.matrix(log2(RMG[rownames(deProt),]))
med_limma = f(a, log_limma)
log_limma[is.na(log_limma)] = 0
rownames(deProt) = all[rownames(deProt),1]
source("hu_dataprep.R")  ## read table and prepare normalized abundances (NA per row < 24); change table to read either with adj.pval or normal p.val
source("heatmap,pca,dist_functions.R") ## get di(), hm(), pca() functions
library(openxlsx)
library(tidyverse)
library(ggrepel)
### general overview of all, filtered normalized abundances
x = med_nabu[,!grepl("_SP",colnames(med_nabu))]
z = log_nabu[,!grepl("_SP",colnames(log_nabu))]
y <-"./figures/hu_all_overview.pdf"
a_all<-data.frame(group=c( rep("NEUR_M", 5),
rep("RMG_M", 5),
rep("NEUR_P", 5),
rep("RMG_P", 5)),
row.names = colnames(x))
acol_all <- list(group=c(RMG_M="#6A94D4",
RMG_P="#FFCB73",
NEUR_M = "#224C8A",
NEUR_P = "#A67625"
))
di(z)
pdf(y, height = 10/2.54, width = 10/2.54)
hm(x, dis, dis2)
pca(z)
plot(dis)
dev.off()
### glia specific genes: adjs.p.val < 0.05 & detected in at least 3 glia samples
RMG_M = all[(all$adj.p.val_RvN_M < 0.05 | is.na(all$adj.p.val_RvN_M))  & rowSums(!is.na(all[,7:11]))>=3,] # RMG vs. Neuron specific in Macula
RMG_P = all[(all$adj.p.val_RvN_P < 0.05 | is.na(all$adj.p.val_RvN_P)) & rowSums(!is.na(all[,17:21]))>=3,] # RMG vs. Neuron specific in periphery
RMG_SP = all[(all$adj.p.val_RvN_SP < 0.05 | is.na(all$adj.p.val_RvN_SP)) & rowSums(!is.na(all[,26:30]))>=3,] #RMG vs. Neuron specific in superperiphery
RMG_M[is.na(RMG_M)] = 0
RMG_P[is.na(RMG_P)] = 0
RMG_SP[is.na(RMG_SP)] = 0
# glia vs. neuron ratio >1
RMG_M = rownames(RMG_M[rowMeans(RMG_M[,7:11], na.rm = T) > rowMeans(RMG_M[,2:6], na.rm = T),])
RMG_P = rownames(RMG_P[rowMeans(RMG_P[,17:21], na.rm = T) > rowMeans(RMG_P[,12:16], na.rm = T) ,])
RMG_SP = rownames(RMG_SP[rowMeans(RMG_SP[,26:30], na.rm = T) > rowMeans(RMG_SP[,22:25], na.rm = T) ,])
### glia-spec, only glia, only macula, periphery
RMG = all[unique(c(RMG_M,RMG_P)),2:30]
RMG = RMG[,grep("RMG", colnames(RMG))]
RMG = RMG[,1:10]
log_RMG = data.matrix(log2(RMG))
med_RMG = f(a, log_RMG)
log_RMG[is.na(log_RMG)] = 0
#write.xlsx(data.frame(RMG,all[rownames(RMG),1]), "rmg.xlsx")
##### explore with limma
library(limma)
l = log2(RMG)
colnames(l) = c(rep("RMG_M",5),rep("RMG_P",5))
fac = factor(colnames(l))
design <- model.matrix(~0+fac)
colnames(design) = c("RMG_M","RMG_P")
fit = lmFit(l, design = design)
contrast.matrix = makeContrasts(RMG_M-RMG_P, levels = design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
deProt = topTable(fit2,number = Inf, adjust="BH")
deProt = deProt[deProt$adj.P.Val<0.05 & !is.na(deProt$adj.P.Val),]
deProt = deProt[order(deProt$logFC, decreasing = T),]
log_limma = data.matrix(log2(RMG[rownames(deProt),]))
med_limma = f(a, log_limma)
log_limma[is.na(log_limma)] = 0
rownames(deProt) = all[rownames(deProt),1]
## write diff proteins
write.xlsx(deProt, file = "./tables/hu_diff_limma_adj.pval.xlsx", row.names =T)
## heatmap
x = med_limma
z = log_limma
y <-"./figures/hu_glia_spec_RMG_only_diff_limma.pdf"
a_all<-data.frame(group=c( rep("RMG_M", 5),
rep("RMG_P", 5)),
row.names = colnames(x))
acol_all <- list(group=c(RMG_M="#6A94D4",
RMG_P="#FFCB73"
))
pdf(y, height = 15, width = 5)
hm(x, F, F)
dev.off()
tb_deProt = topTable(fit2,number = Inf, adjust="BH")
tb_deProt = as_tibble(tb_deProt) %>% mutate(prots = rownames(tb_deProt))
tb_deProt = tb_deProt %>% mutate(pval = ifelse(adj.P.Val<0.05, T, F),
FC = ifelse(abs(logFC)>1, T,F),
label = ifelse(adj.P.Val<0.01 & abs(logFC)>3, all[prots,1],NA)) %>%
mutate(reg = ifelse(pval & FC & logFC>0, "M", "NA")) %>%
mutate(reg = ifelse(pval & FC & logFC<0, "P", reg))
vlc = tb_deProt %>%
ggplot(aes(x = logFC, y = -log10(adj.P.Val), color = reg, label = label )) +
geom_hline(yintercept = -log10(0.05), color = "darkgrey") +
geom_vline(xintercept = c(-1,1), color = "darkgrey") +
geom_point() +
geom_text_repel()+
scale_color_manual(values = c("#FFCB73","#6A94D4", "black"), limits = c("P", "M","NA"))+
theme_light() +
theme(panel.grid= element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
axis.line = element_line(color = "black"),
text =  element_text(size = 12),
legend.position = "none") +
labs(x = "logFC" , y = "-log10(adj.P.Val.)")
pdf(file = "./figures/hu_vulcano_limma.pdf", width = 10/2.54,height = 10/2.54)
vlc
dev.off()
source("ms_dataprep.R")  ## read table and prepare normalized abundances
source("heatmap,pca,dist_functions.R") ## get di(), hm(), pca() functions
library(openxlsx)
library(tidyverse)
library(ggrepel)
cnames = data.frame(cnames) ##cnames defined in dataprep.R
samples = paste(cnames[1:40, "celltype"], cnames[1:40, "gt"], sep="_")
