书接上文，gsva后续的处理很简单，如果熟悉差异基因分析包limma的话，更是简单。之前我们得到的gsva分数的矩阵就类似于基因表达矩阵，按照这个思路继续往下即可：
从通路的表达矩阵开始，我们进行差异分析，首先将之前保存的文件读入：
T_gsva <- read.csv("T_gsva.csv", header = T, row.names = 1)
分析得到的数据结构是行为GO terms， 列为样品（单细胞中为单个细胞）。
group <- c(rep("control", 50), rep("test", 71)) %>% as.factor()#设置分组，对照在前
desigN <- model.matrix(~ 0 + group) #构建比较矩阵
colnames(desigN) <- levels(group)
fit = lmFit(test_control, desigN)
fit2 <- eBayes(fit)
diff=topTable(fit2,adjust='fdr',coef=2,number=Inf)#校准按照需求修改 ？topTable
write.csv(diff, file = "Diff.csv")
差异分析结果有logFC，P-value，t值等等，根据阈值筛选差异的通路即可，最后挑选与自己研究相关的或者感兴趣的通路进行可视化。可视化的类型有热图展示，个人认为不好看，因为细胞数太多，也有用平均值做的，效果一般！
这里看到了一篇文章中的可视化，感觉不错，大多数文章也是这种类型的柱状图！
1629453025(1).png
文章引用：[1] Lambrechts, D. , et al. "Phenotype molding of stromal cells in the lung tumor microenvironment." Nature Medicine 24.8(2018).
up <- c("GOBP_EGG_ACTIVATION",
"GOBP_TENDON_DEVELOPMENT",
"GOBP_SOMITE_SPECIFICATION",
"GOBP_THREONINE_CATABOLIC_PROCESS",
"GOBP_REGULATION_OF_GLUTAMATE_RECEPTOR_CLUSTERING",
"GOBP_NEGATIVE_CHEMOTAXIS",
"GOBP_NEGATIVE_REGULATION_OF_FAT_CELL_PROLIFERATION",
"GOBP_REGULATION_OF_T_HELPER_17_CELL_LINEAGE_COMMITMENT",
"GOBP_REGULATION_OF_ANTIMICROBIAL_HUMORAL_RESPONSE")
down <- c("GOBP_DETERMINATION_OF_PANCREATIC_LEFT_RIGHT_ASYMMETRY",
"GOBP_MITOTIC_DNA_REPLICATION",
"GOBP_EOSINOPHIL_CHEMOTAXIS",
"GOBP_NEUTROPHIL_MEDIATED_CYTOTOXICITY",
"GOBP_POTASSIUM_ION_EXPORT_ACROSS_PLASMA_MEMBRANE",
"GOBP_REGULATION_OF_LEUKOCYTE_MEDIATED_CYTOTOXICITY",
"GOBP_REGULATION_OF_SEQUESTERING_OF_ZINC_ION",
"GOBP_ENDOTHELIN_RECEPTOR_SIGNALING_PATHWAY",
"GOBP_PRE_REPLICATIVE_COMPLEX_ASSEMBLY_INVOLVED_IN_CELL_CYCLE_DNA_REPLICATION",
"GOBP_ESTABLISHMENT_OF_PLANAR_POLARITY_OF_EMBRYONIC_EPITHELIUM")
TEST <- c(up,down)
diff$ID <- rownames(diff)
Q <- diff[TEST,]
group1 <- c(rep("treat", 9), rep("control", 10))
df <- data.frame(ID = Q$ID, score = Q$t,group=group1 )
# 按照t score排序
sortdf <- df[order(df$score),]
sortdf$ID <- factor(sortdf$ID, levels = sortdf$ID)#增加通路ID那一列
用ggplot画图（ggplot-YYDS）
ggplot(sortdf, aes(ID, score,fill=group)) + geom_bar(stat = 'identity',alpha = 0.7) +
coord_flip() +
theme_bw() + #去除背景色
theme(panel.grid =element_blank())+
theme(panel.border = element_rect(size = 0.6))+
labs(x = "",
y="t value of GSVA score")+
scale_fill_manual(values = c("#008020","#08519C"))#设置颜色
1629453814(1).png
有闲工夫了可以自己修饰修饰即可！
GSVA分析到此结束！
欢迎分享交流！
网友评论