说在前面
Immugent在之前一篇推文照葫芦画图之PCA中介绍了PCA的分析，无独有偶，最近又看了一篇文章中用的是PcoA分析，今天就来给大家介绍一下如何做这种图。
PCoA（Principal Co-ordinates Analysis）是一种研究数据相似性或差异性的可视化方法，通过一系列的特征值和特征向量进行排序后，选择主要排在前几位的特征值；通过PCoA分析 可以观察个体或群体间的差异。
图片
原文发表在Cell Metab(IF:27.29)，主要是用PcoA对肠道菌的信息进行降维，此外还使用PERMANOVA test来检测了各组间的差异，这个流程非常值得大家借鉴。
数据准备
因为大家研究肿瘤的较多，今天小编就来利用肿瘤的数据来进行展示，所用数据是乳腺癌的TCGA数据和Christina Yau整合的乳腺癌芯片数据。
library(ape) #做PCoA
library(RColorBrewer) #配色
library(usedist) #PERMANOVA test时处理距离矩阵
Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor
# 读取表达矩阵
expr <- read.table("easy_input_expr.txt",sep = "\t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)
# 读取样本信息
sinfo <- read.table("easy_input_sinfo.txt",sep = "\t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)
# 查看亚型
table(sinfo$PAM50) # 这里有5个亚型
# 查看队列
table(sinfo$Cohort) # 2个队列
## TCGA Yau
## 1025 682
代码实现
首先画一个PCA图
# 创建绘图信息
plotinfo <- cbind.data.frame(x = pc[,1],
y = pc[,2],
PAM50 = sinfo[rownames(pc),"PAM50"],
Cohort = sinfo[rownames(pc),"Cohort"],
Class = paste(sinfo[rownames(pc),"PAM50"], # 根据class信息确定颜色
sinfo[rownames(pc),"Cohort"],
sep = "_"))
# 为了使来自两个队列相同亚型的颜色相近，像这样修改因子，使得配对颜色可以对应上
plotinfo$Class <- factor(plotinfo$Class,
levels = c("Basal_TCGA","Basal_Yau",
"Her2_TCGA","Her2_Yau",
"LumA_TCGA","LumA_Yau",
"LumB_TCGA","LumB_Yau",
"Normal_TCGA","Normal_Yau"))
# 10个组（5个亚型*2个队列），定义10种颜色，采用配对颜色分别表示TCGA和Yau对应亚型
mycol <- brewer.pal(n = 10, name = "Paired") # Paired最多可设置12组颜色
plotinfo$color <- mycol[plotinfo$Class] # 匹配颜色
# 保存到文件
write.csv(plotinfo, "output_PCA1_PCA2.csv", quote = F)
# 产生基本图看一下效果，可以看到不同的亚型还是可以分开的
pdf("basic scatter plot.pdf", width = 5,height = 5)
par(bty="o", mgp = c(1.9,.33,0), mar=c(3.1,3.1,2.1,2.1)+.1, las=1, tcl=-.25)
plot(plotinfo$x,
plotinfo$y,
pch = 19,
col = plotinfo$color,
xlab = xlab.text,
ylab = ylab.text)
dev.off()
图片
可以看出分群的效果并不好，于是尝试做PcoA分析，并且对计算相同组的取均值，以及算出标准误后再画图。
plotinfo2 <- NULL
for (i in unique(plotinfo$Class)) {
tmp <- plotinfo[plotinfo$Class == i,] # 取出当前亚型当前来源下的数据
avgx <- mean(tmp$x) # 计算横坐标均值
avgy <- mean(tmp$y) # 计算纵坐标均值
sdx <- sd(tmp$x) # 计算横坐标标准差
sdy <- sd(tmp$y) # 计算纵坐标标准差
plotinfo2 <- rbind.data.frame(plotinfo2,
data.frame(Class = i,
color = unique(tmp$color), # 添加颜色
shape = ifelse(unique(tmp$Cohort) == "TCGA", # 如果是TCGA就实心圆
"closed","opened"), # 如果是Yau就为空心圆
label = switch(i, # 创建简易标签以显示在圆中
"Her2_TCGA" = "HT",
"Normal_TCGA" = "NT",
"Basal_TCGA" = "BT",
"LumA_TCGA" = "LaT",
"LumB_TCGA" = "LbT",
"Her2_Yau" = "HY",
"Normal_Yau" = "NY",
"Basal_Yau" = "BY",
"LumA_Yau" = "LaY",
"LumB_Yau" = "LbY"),
avgx = avgx, # 添加圆的x位置
avgy = avgy, # 添加圆的y位置
sdx = sdx, # 添加圆的水平标准差
sdy = sdy, # 添加圆的垂直标准差
stringsAsFactors = F),
stringsAsFactors = F)
}
# 保存到文件
write.csv(plotinfo2, "output_av_sd.csv", quote = F)
# 绘图
pdf("PoCA.pdf", width = 5,height = 5)
par(bty="o", mgp = c(1.9,.33,0), mar=c(3.1,3.1,2.1,2.1)+.1, las=1, tcl=-.25)
# 根据上面的基础图像，调整横轴纵轴的宽度，生成一张空白图像
plot(NULL,NULL,
xlab = xlab.text,
ylab = ylab.text,
xlim = c(-25,25),
ylim = c(-20,20))
# 先产生误差线，这样可以被后面的圆挡住
for (i in 1:nrow(plotinfo2)) {
tmp <- plotinfo2[i,]
# 产生横向误差线
lines(x = c(tmp$avgx - tmp$sdx,
tmp$avgx + tmp$sdx),
y = c(tmp$avgy, tmp$avgy),
lty = ifelse(tmp$shape == "closed",1,2), # 如果是closed就是实线，否则为虚线
col = tmp$color,
lwd = 2)
# 产生横向误差线的封口线
lines(x = c(tmp$avgx - tmp$sdx,
tmp$avgx - tmp$sdx),
y = c(tmp$avgy - 0.5, # 宽度根据情况调整
tmp$avgy + 0.5),
col = tmp$color,
lwd = 2)
lines(x = c(tmp$avgx + tmp$sdx,
tmp$avgx + tmp$sdx),
y = c(tmp$avgy - 0.5,
tmp$avgy + 0.5),
col = tmp$color,
lwd = 2)
# 产生纵向误差线
lines(x = c(tmp$avgx, tmp$avgx),
y = c(tmp$avgy - tmp$sdy,
tmp$avgy + tmp$sdy),
lty = ifelse(tmp$shape == "closed",1,2),
col = tmp$color,
lwd = 2)
# 产生纵向误差线的封口线
lines(x = c(tmp$avgx - 0.5,
tmp$avgx + 0.5),
y = c(tmp$avgy + tmp$sdy,
tmp$avgy + tmp$sdy),
col = tmp$color,
lwd = 2)
lines(x = c(tmp$avgx - 0.5,
tmp$avgx + 0.5),
y = c(tmp$avgy - tmp$sdy,
tmp$avgy - tmp$sdy),
col = tmp$color,
lwd = 2)
}
# 后添加圆，以挡住误差线
points(plotinfo2$avgx,
plotinfo2$avgy,
pch = ifelse(plotinfo2$shape == "closed",19,21), # 如果为closed就是实心圆，否则为空心圆
bg = ifelse(plotinfo2$shape == "closed",plotinfo2$color,"white"), # 填充背景色，如果为closed就是为该颜色，否则为白色
col = plotinfo2$color, # 填充边框颜色
lwd = 2, # 边框粗细
cex = 5) # 圆的大小
# 添加文本
text(plotinfo2$avgx,
plotinfo2$avgy,
plotinfo2$label,
col = ifelse(plotinfo2$shape == "closed","white",plotinfo2$color), # 如果是实心圆文字为白色，否则为对应颜色
cex = 1.1)
dev.off()
图片
这样看各亚型之间就能很好的区分了，下面再按照作者推荐的PERMANOVA test检验一下它们之间的差异。
## TCGA队列数据
tcga.sam <- rownames(sinfo[sinfo$Cohort == "TCGA",])
# 抽取tcga样本的距离矩阵
tcga.dist <- dist_subset(brca.dist,tcga.sam)
# PERMANOVA test
brca.TCGA <- adonis2(tcga.dist~PAM50, data=sinfo[tcga.sam,],
permutations = 1000, method="bray")
print(brca.TCGA)
## adonis2(formula = tcga.dist ~ PAM50, data = sinfo[tcga.sam, ], permutations = 1000, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## PAM50 4 152749 0.17185 52.916 0.000999 ***
## Yau队列数据，也得到类似的结果
yau.sam <- rownames(sinfo[sinfo$Cohort == "Yau",])
yau.dist <- dist_subset(brca.dist,yau.sam)
brca.Yau <- adonis2(yau.dist~PAM50, data=sinfo[yau.sam,],
permutations = 1000, method="bray")
print(brca.Yau)
## adonis2(formula = yau.dist ~ PAM50, data = sinfo[yau.sam, ], permutations = 1000, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## PAM50 4 213707 0.16525 33.505 0.000999 ***
总结
PERMANOVA test是一种生态学中常用的统计方法，主要是对样品组间差异显著性进行检验。这是一种非参数检验，不需要满足正态分布，大家可以放心使用。本推文分别以两个队列（TCGA和Yau）的各亚组之间比较为例，简单介绍了它的用法。
本期到这就结束啦，后面生信宝库还会推出系列对高分SCI图片的复现，敬请期待！
