1. 简介
AAV与免疫细胞关系密切，在该疾病发生时，各类免疫细胞会被异常激活，探究AAV发生发展过程中免疫细胞的变化情况及其它们之间的相互关系具有重要意义。
2. 数据信息
- 关键基因的表达量数据
- 样本分组信息
3. 思路
- 使用反卷积算法对转录组训练集中所有样本的免疫细胞浸润情况进行评估。
- 利用 Wilcoxon 秩和检验 AAV 组和 Control 组样本之间免疫细胞的差异情况，将两组间具有显著性差异的免疫细胞记为差异免疫细胞，并可视化。
- 利用R包进行Spearman相关性分析差异免疫细胞间的相关性，并使用R包绘制热图展示差异免疫细胞的相关性。
- 利用R包进行Spearman相关性分析转录组训练集所有样本中关键基因与差异免疫细胞之间的相关性，并使用R包绘制热图展示结果。
4. 代码
library(IOBR)
library(ggplot2)
library(reshape2)
library(ggpubr)
library(dplyr)
library(psych)
exp_train <- read.csv("../00_DataPre/exp_train.csv", row.names = 1)
train_group <- read.csv("../00_DataPre/train_group.csv", row.names = 1)
colnames(train_group) <- c("ID","group")
ciber_res <- deconvo_tme(eset = exp_train, method = "cibersort",tumor = FALSE,arrays = TRUE, perm = 100)
save(ciber_res,file = "CIBERSORT.RData")
load("CIBERSORT.RData")
ciber_res <- as.data.frame(ciber_res)
colnames(ciber_res) <- c("ID","Naive B cells","Memory B cells","Plasma cells","CD8 T cells",
"Naive CD4 T cells","Resting memory CD4 T cells","Activated memory CD4 T cells",
"Follicular helper T cells","Regulatory(Tregs) T cells","Gamma delta T cells",
"Resting NK cells","Activated NK cells","Monocytes","M0 Macrophages",
"M1 Macrophages","M2 Macrophages","Resting dendritic cells","Activated dendritic cells",
"Resting mast cells","Activated mast cells","Eosinophils","Neutrophils",
"P-value","Correlation","RMSE")
ciber_res <- ciber_res[,-c(24:26)]
rownames(ciber_res) <- ciber_res$ID
ciber_res$ID <- NULL
ciber_res$group <- train_group$group[match(rownames(ciber_res),train_group$ID)]
ciber_data <- ciber_res %>% tibble::rownames_to_column('sample') %>% tidyr::gather(.,key = 'cell',value = 'score',c(-group,-sample))
ggplot()+
geom_col(data = ciber_data,mapping = aes(x=sample,y=score,fill=cell),position = 'fill',alpha=0.7)+
labs(x='',y='Estimated proportion',fill='Cell type')+
guides(fill=guide_legend(ncol = 1))+
theme(axis.text.x = element_blank(),
axis.text.y = element_text(size = 15,color = 'black'),
axis.title = element_text(size = 18,color = 'black'),
axis.ticks.x = element_blank(),
panel.background = element_rect(fill = NA,color = NA),
legend.title = element_text(size = 18,colour = 'black'),
legend.text = element_text(size = 12,colour = 'black'),
legend.key.spacing.y = unit(0.1,'mm'),
legend.key.size = unit(0.1,'mm'),
strip.text = element_text(size = 15),
strip.background = element_rect(fill = 'gray',colour = 'black'))+
facet_grid(. ~ group, scales = "free")
h <- 7;w <- 9
ggsave("09.barplot_of_cell_ratio.pdf",height=h,width=w)
ggsave("09.barplot_of_cell_ratio.png",height=h,width=w)
#ciber_res$ID <- train_group$group
#pdf("11_CIBERSORT/01.barplot_of_cell_ratio.pdf",w = 10,h = 50,family = 'Times')
#ciber_plot <- cell_bar_plot(input = ciber_res, feature = colnames(ciber_res)[2:23],
# title = "CIBERSORT Cell Fraction",id = NULL)
#dev.off()
#png("11_CIBERSORT/01.barplot_of_cell_ratio.png",w = 10,h = 50,units = "in",res = 300, family = 'Times')
#ciber_plot <- cell_bar_plot(input = ciber_res, feature = colnames(ciber_res)[2:23],
# title = "CIBERSORT Cell Fraction")
#dev.off()
load("CIBERSORT.RData")
colnames(ciber_res) <- c("ID","Naive B cells","Memory B cells","Plasma cells","CD8 T cells",
"Naive CD4 T cells","Resting memory CD4 T cells","Activated memory CD4 T cells",
"Follicular helper T cells","Regulatory(Tregs) T cells","Gamma delta T cells",
"Resting NK cells","Activated NK cells","Monocytes","M0 Macrophages",
"M1 Macrophages","M2 Macrophages","Resting dendritic cells","Activated dendritic cells",
"Resting mast cells","Activated mast cells","Eosinophils","Neutrophils",
"P-value","Correlation","RMSE")
ciber_res <- ciber_res[,-c(24:26)]
ciber_res$ID <- train_group$group
ciber = melt(ciber_res)
colnames(ciber) <- c("group","celltype","composition")
ggplot(ciber, aes(x = celltype, y = composition)) +
labs(y = "Cell composition",y = NULL,title = "Immune cell composition") +
geom_boxplot(aes(fill = group),position = position_dodge(0.5),width = 0.5,outlier.alpha = 0) +
scale_fill_manual(values = c("#1CB4B8", "#EB7369"))+
theme_classic(base_size = 14) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
stat_compare_means(aes(group = group),
label = "p.signif",
method = "wilcox.test",
hide.ns = F)
ggsave("09.boxplot.png",w = 15,h = 8,dpi = 300)
ggsave("09.boxplot.pdf",w = 15,h = 8,dpi = 300)
sig_immune <- c("Memory B cells","CD8 T cells","Regulatory(Tregs) T cells",
"M1 Macrophages","M2 Macrophages","Resting dendritic cells",
"Resting mast cells","Activated mast cells","Neutrophils")
corr_result <- cor(ciber_res[,sig_immune],method = "spearman")
corr_result[abs(corr_result)<0.3] <- NA
corr_result_p <- WGCNA::corPvalueStudent(corr_result,nrow(ciber_res))
corr_result_p[corr_result_p >= 0.05] <- NA
corr_result[is.na(corr_result_p)] <- NA
library(corrplot)
col2 <- colorRampPalette(c("#053061","#2166AC", "#4393C3","#92C5DE","#D1E5F0","#FDDBC7","#F4A582","#D6604D","#B2182B","#67001F"))(50)
corr_plot <- function(){
corrplot(corr_result,
na.label = "no",
diag=F,
method="number",
tl.col="black",
tl.srt=45,
tl.cex=1,
cl.length = 8,
cl.cex=0.6,
col = col2,
number.cex=1
)
corrplot(corr_result,
p.mat = corr_result_p,
insig = "label_sig",
add=TRUE,
col = col2,
na.label = "ns",
tl.pos="n",
cl.pos="n",
sig.level = c(.0001,.001, .01, .05),
pch.cex = 1,
number.cex=1,
pch.col = "black"
)
}
pdf("09.sig_immune_cells_heatmap.pdf",w = 8,h = 6, family = 'Times')
corr_plot()
dev.off()
png("09.sig_immune_cells_heatmap.png",w = 8,h = 6,units = "in",res = 300, family = 'Times')
corr_plot()
dev.off()
##########
## biomarker: "VCAN","CD74"
geneset <- c("VCAN","CD74")
exp_train1 <- as.data.frame(t(exp_train))
exp_train2 <- exp_train1[,geneset]
load("CIBERSORT.RData")
colnames(ciber_res) <- c("ID","Naive B cells","Memory B cells","Plasma cells","CD8 T cells",
"Naive CD4 T cells","Resting memory CD4 T cells","Activated memory CD4 T cells",
"Follicular helper T cells","Regulatory(Tregs) T cells","Gamma delta T cells",
"Resting NK cells","Activated NK cells","Monocytes","M0 Macrophages",
"M1 Macrophages","M2 Macrophages","Resting dendritic cells","Activated dendritic cells",
"Resting mast cells","Activated mast cells","Eosinophils","Neutrophils",
"P-value","Correlation","RMSE","score")
sample_name <- ciber_res$ID
ciber_res$ID <- NULL
rownames(ciber_res) <- sample_name
ct1<- corr.test(exp_train2,ciber_res[,sig_immune],method = "spearman")
r1<- ct1$r
corr_result <- cor(exp_train2,ciber_res[,sig_immune],method = "spearman")
corr_result[abs(corr_result)<0.3] <- NA
corr_result_p <- WGCNA::corPvalueStudent(corr_result,nrow(ciber_res))
corr_result_p[corr_result_p >= 0.05] <- NA
corr_result[is.na(corr_result_p)] <- NA
col2 <- colorRampPalette(c("#053061","#2166AC", "#4393C3","#92C5DE","#D1E5F0","#FDDBC7","#F4A582","#D6604D","#B2182B","#67001F"))(50)
corr_plot <- function(){
corrplot(corr_result,
na.label = "ns",#NA的标签用什么表示
diag=F,#是否展示对角线数据
method="number",#展示图类型
# type = 'upper',
# tl.pos="lt",#文本标签的位置
tl.col="black",#文本标签的颜色
tl.srt=45,#文本标签的旋转度数
tl.cex=1,#文本标签字体大小
cl.length = 8,#颜色图例中的刻度数
cl.cex=0.6,#颜色图例字体大小
col = col2,#图的颜色
number.cex=1
)
corrplot(corr_result,
p.mat = corr_result_p,
# method="square",#展示图类型
# type = 'lower',
insig = "label_sig",
add=TRUE,#是否将该图添加到已存在的图
col = col2,
na.label = "ns",
# diag=F,
tl.pos="n",
cl.pos="n",
sig.level = c(.0001,.001, .01, .05),
pch.cex = 1,
number.cex=1,
pch.col = "black"
)
}
pdf("09.biomarker_sig_immune_heatmap.pdf",w = 6,h = 4, family = 'Times')
corr_plot()
dev.off()
png("09.biomarker_sig_immune_heatmap.png",w = 6,h = 4,units = "in",res = 300, family = 'Times')
corr_plot()
dev.off()
5. 结果展示
09.barplot_of_cell_ratio.png
09.boxplot.png
09.sig_immune_cells_heatmap.png
09.biomarker_sig_immune_heatmap.png
