【后续来了】有了这个包，猪的GSEA和GSVA分析也不在话下（第二集）

书接上文，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值等等，根据阈值筛选差异的通路即可，最后挑选与自己研究相关的或者感兴趣的通路进行可视化。可视化的类型有热图展示，个人认为不好看，因为细胞数太多，也有用平均值做的，效果一般！
这里看到了一篇文章中的可视化，感觉不错，大多数文章也是这种类型的柱状图！

文章引用：[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"))#设置颜色

有闲工夫了可以自己修饰修饰即可！
GSVA分析到此结束！
欢迎分享交流！