今晚回到办公室，看到关于绘制FPI的计算方法的文章，里面展示的图也是比较喜欢，那么就分享一下吧。其实这样的图是非常实用的，可以使用在各个分析结果中。

image.png

一、数据的输入`

#'@加载R包
library(data.table)
library(GSVA)
library(ggplot2)
Sys.setenv(LANGUAGE ='en')  # #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor

 df <- read.csv("input.data.csv",header = T, row.names = 1)

##---
>head(df)
                       fpi tissue            p   p.lab TCGA
TCGA-BP-4989-01 -2.0539160  Tumor 4.724241e-22 4.7e-22 KIRC
TCGA-CZ-5459-01  0.2892758  Tumor 4.724241e-22 4.7e-22 KIRC
TCGA-A3-3335-01 -0.2936776  Tumor 4.724241e-22 4.7e-22 KIRC
TCGA-CJ-6028-01  0.4542248  Tumor 4.724241e-22 4.7e-22 KIRC
TCGA-CZ-5456-01  1.0573740  Tumor 4.724241e-22 4.7e-22 KIRC
TCGA-A3-3383-01  0.8313845  Tumor 4.724241e-22 4.7e-22 KIRC

二、绘图

ggplot(df, aes(TCGA, fpi, fill=tissue)) + 
  geom_boxplot(outlier.shape = NA) + 
  geom_text(aes(TCGA, y = 4,
                label = paste0("P = ", p.lab)),
            data = df,
            inherit.aes = F) +
  scale_fill_manual(values = c(green, yellow)) + 
  scale_y_continuous(breaks = c(-2.5,0,2.5,5), labels = c(-2.5,0,2.5,5), limits = c(-3,5)) + 
  xlab(NULL) + ylab("FPI") + coord_flip() +
  theme_classic()

image.png

三、FPI的计算

#'@输入你的数据及表达量和验证集数据
> head(gse.expr)
                       GSM3443357 GSM3443358 GSM3443359 GSM3443360 GSM3443361 GSM3443362 GSM3443363 GSM3443364 GSM3443365
APITD1 /// APITD1-CORT   8.441097   8.475046   8.655390   8.516031   8.418268   8.460168   8.282255   8.287272   8.186783
SAA2-SAA4 /// SAA4       5.252111   5.430918   5.420259   5.496871   5.673105   5.437256   5.603959   5.337072   5.346030
ORM1 /// ORM2            3.078793   3.091741   3.110379   2.997581   3.046607   2.879500   3.129926   3.055256   2.965808
T                        5.100341   5.147388   5.031716   5.212217   5.170536   5.037246   5.128167   5.176053   5.243752
DHFR                    10.410754  10.326273  10.434939  10.458324  10.341043  10.319331  10.233846   9.912388   9.979284
LEFTY1 /// LEFTY2        4.722003   4.799446   4.825212   4.819980   4.912926   4.549273   4.723768   4.376253   4.768935
                       GSM3443366 GSM3443367 GSM3443368
APITD1 /// APITD1-CORT   8.287317   8.325607   8.459985
SAA2-SAA4 /// SAA4       5.307920   5.499212   5.281267
ORM1 /// ORM2            2.997028   2.966939   2.993494
T                        5.285633   5.204160   5.137498
DHFR                    10.008167  10.054794  10.005322
LEFTY1 /// LEFTY2        4.855704   4.642899   4.699271

四、开始计算

# ssgsea计算enrichment score
es.gse <- gsva(expr = as.matrix(gse.expr),
               gset.idx.list = fpi.sig,
               method = "ssgsea",
               ssgsea.norm = TRUE)
fpi.gse <- as.numeric(scale(es.gse[1,] - es.gse[2,])); names(fpi.gse) <- colnames(gse.expr)
wt <- wilcox.test(fpi.gse[dmso.sam],fpi.gse[erastin.sam]) # 秩和检验

五、绘图

boxplot(fpi.gse[dmso.sam],
        fpi.gse[erastin.sam],
        col = c("#2E8049","#DD8606"),
        ylab = "FPI (you data file name)",
        xlab = "",
        ylim = c(-2,2), # 控制y轴区间
        names = c("DMSO","Erastin"))
lines(c(1,2),c(2.2,2.2)) # 在两个箱子之间加横线
lines(c(1,1),c(2.1,2.2)) # 垂直于第一个箱子加短竖线
lines(c(2,2),c(2.1,2.2)) # 垂直于第二个箱子加短竖线
text(1.5, 2.3, paste0("p = ",round(wt$p.value, 3))) # 添加p值

image.png

往期文章（总汇，点击链接即可进入）

小杜的生信筆記 ，主要发表或收录生物信息学的教程，以及基于R的分析和可视化（包括数据分析，图形绘制等）；分享感兴趣的文献和学习资料！