今晚回到办公室,看到关于绘制FPI的计算方法的文章,里面展示的图也是比较喜欢,那么就分享一下吧。其实这样的图是非常实用的,可以使用在各个分析结果中。
一、数据的输入`
#'@加载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()
三、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值
往期文章(总汇,点击链接即可进入)
小杜的生信筆記 ,主要发表或收录生物信息学的教程,以及基于R的分析和可视化(包括数据分析,图形绘制等);分享感兴趣的文献和学习资料!