PCA
在RNA-seq中，主成分分析（PCA）是最常见的多元数据分析类型之一。基因表达定量后获得了各样本中所有基因的表达值信息，随后我们通常会期望比较样本之间在基因表达值的整体相似性或者差异程度。基因数量成千上万，肯定不能对每个基因的表达都作个比较，这时候就要用到“降维”算法，PCA分析因此派上用场。PCA设法将N维（N=基因数量）的表达矩阵降维成只有少数几个主成分的形式，这几个主成分就可以代表所有基因的整体表达格局，进而据此描述样本差异。
在PCA图中，如果两个样本间整体基因表达值更为相似，那么它们的距离就接近；反之若整体基因表达值差异很大，则它们的距离就会更远。据此，我们就可以从中评估，不同组间样本的基因表达是否相差更为明显，组内样本的基因表达是否均匀或一致性较高，哪些处理组之间引起了相似的基因表达变化趋势等。
本篇教程就让我们来学习如何使用R语言画带基因名字标签PCA图。
本文涉及到的测试数据及代码已经给大家打包好，如下：
image.png
1. 示例数据
1. “PCA.data.txt” 为基因表达值矩阵。其中第一列为基因名称，这里以ensembl id作为指代；其余各列记录了RNA-seq获得的各基因在各样本中的表达量信息。
All_gene_fpkm.list
2. “group.txt” 则为样本分组文件，记录了样本所属的不同分组。
group.txt
2. R代码部分
代码重要部分已被删除，完整代码+数据需要私信要！！！
代码重要部分已被删除，完整代码+数据需要私信要！！！
代码重要部分已被删除，完整代码+数据需要私信要！！！
#读取基因表达值矩阵
#推荐使用 log 转化后的基因表达值，降低不同基因表达水平数量级相差过大的问题
gene <- read.delim('All_gene_fpkm.list', row.names = 1, sep = '\t',check.names = FALSE)
#将基因表达值矩阵作个转置，使行为样本，列为基因
gene <- t(gene)
#我们使用 FactoMineR 包中的方法，实现 PCA 分析和聚类添加
library(FactoMineR)
#样本中基因表达值的 PCA 分析
gene.pca <- PCA(gene, ncp = 2, scale.unit = TRUE, graph = FALSE)
#plot(gene.pca) #PCA 简图
#提取样本在 PCA 前两轴中的坐标
pca_sample <- data.frame(gene.pca$ind$coord[ ,1:2])
head(pca_sample)
#提取 PCA 前两轴的贡献度
pca_eig1 <- round(gene.pca$eig[1,2], 2)
pca_eig2 <- round(gene.pca$eig[2,2],2 )
#读取并合并样本分组信息
group <- read.delim('group.txt', row.names = 1, sep = '\t', check.names = FALSE)
group <- group[rownames(pca_sample), ]
#ggplot2 绘制二维散点图
library(ggplot2)
p <- ggplot(data = pca_sample, aes(x = Dim.1, y = Dim.2)) +
geom_point(aes(color = group), size = 5) + #根据样本坐标绘制二维散点图
scale_color_manual(values = c('orange', 'purple','blue','black','green','yellow','red')) + #自定义颜色
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'),
legend.key = element_rect(fill = 'transparent')) + #去除背景和网格线
labs(x = paste('PCA1:', pca_eig1, '%'), y = paste('PCA2:', pca_eig2, '%'), color = '') #将 PCA 轴贡献度添加到坐标轴标题中
#输出图1
#标识样本名称，使用 ggplot2 的拓展包 ggrepel 来完成
#带样本名称标签
library(ggrepel)
p
#输出图2
#添加 95% 置信椭圆，可用于表示对象分类，但只能适用于各组样本数大于 5 个的情况
#置信区间线条
p + stat_ellipse(aes(color = group), level = 0.95, show.legend = FALSE)
#输出图3
#置信区间阴影
p + stat_ellipse(aes(fill = group), geom = 'polygon', level = 0.95, alpha = 0.1, show.legend = FALSE) +
scale_fill_manual(values = c('orange', 'purple','blue','black','green','yellow','red'))
#输出图4
#多边形连接同类别对象边界的样式，适用于各组样本数大于 3 个的情况
#多边形线条
library(plyr)
cluster_border <- ddply(pca_sample, 'group', function(df) df[chull(df[[1]], df[[2]]), ])
p + geom_polygon(data = cluster_border, aes(color = group), fill = NA, show.legend = FALSE)
#输出图5
#多边形阴影
p + geom_polygon(data = cluster_border, aes(fill = group), alpha = 0.2, show.legend = FALSE) +
scale_fill_manual(values = c('orange', 'purple','blue','black','green','yellow','red'))
3. 图片输出结果
pca1.png
pca2.png
pca3.png
pca4.png
pca5.png
网友评论