主成分分析（PCA）是一种非常强大的技术，广泛应用于数据科学、生物信息学和其他领域。它最初是为了分析大量数据而开发的，以便揭示被分析的逻辑实体之间的差异/关系。它提取数据的基本结构，而不需要构建任何模型来表示它。通过一种简化过程，可以将大量变量转换成较少的不相关变量（即“主成分”），同时又能够便于对原始数据进行解释，从而得到数据的主要信息。
安装PCAtools
if (!requireNamespace('BiocManager', quietly = TRUE))
install.packages('BiocManager')
BiocManager::install('PCAtools')
#加载包
library(PCAtools)
数据预处理
library(airway)
library(magrittr)
data('airway')
#将dex列中的'untrt'（未处理）水平设为基准水平。
airway$dex %<>% relevel('untrt')
加载“airway”数据集，这个数据集包含了不同的气道平滑肌细胞在被地塞米松处理后的信息。
对数据进行归一化，并将归一化后的计数转换为方差稳定的表达水平：
library('DESeq2')
#构建了一个考虑了细胞类型和处理条件的差异表达分析的数据集对象
dds <- DESeqDataSet(airway, design = ~ cell + dex)
#执行差异表达分析
dds <- DESeq(dds)
#用vst函数对dds对象进行方差稳定转换
vst <- assay(vst(dds))
进行主成分分析
p <- pca(vst, metadata = colData(airway), removeVar = 0.1)
removeVar = 0.1: 这个参数看起来是用来设置一个阈值，移除方差小于该阈值的变量。这可以作为一种数据预处理步骤，以去除那些在各个样本中变化不大（不太重要）的基因或特征，从而可能改善PCA的结果和解释性。
散点图
用于显示主成分分析（PCA）中每个主成分的重要性或贡献度。它有助于确定应考虑多少个主成分以捕获数据中的大部分变异性。
在scree plot中，x轴代表主成分（通常按它们的特征值从大到小排序），y轴代表相应的特征值或每个主成分解释的方差。这个图表通常以陡峭的斜率开始，表明最初的几个成分解释了总方差的重要部分。随着沿x轴向右移动，斜率趋于变平，表明额外的主成分只解释了较小一部分的方差。
screeplot(p, axisLabSize = 18, titleLabSize = 22)
bi-plot
bi-plot它结合了主成分分析（PCA）得分的散点图和表示变量在主成分上载荷的向量。本质上，双标图允许你在单一图表中同时可视化得分（样本之间的关系）和载荷（变量对成分的贡献）。通常只关注前两个主成分。
biplot(p, showLoadings = TRUE,
labSize = 5, pointSize = 5, sizeLoadingsNames = 5)
从GEO数据库下载数据进行深入理解
library(Biobase)
library(GEOquery)
# 从GEO加载系列和平台数据
gset <- getGEO('GSE2990', GSEMatrix = TRUE, getGPL = FALSE)
mat <- exprs(gset[[1]])
# 移除Affymetrix控制探针
mat <- mat[-grep('^AFFX', rownames(mat)),]
# 从表型数据（pdata）中提取感兴趣的信息
idx <- which(colnames(pData(gset[[1]])) %in%
c('relation', 'age:ch1', 'distant rfs:ch1', 'er:ch1',
'ggi:ch1', 'grade:ch1', 'size:ch1',
'time rfs:ch1'))
metadata <- data.frame(pData(gset[[1]])[,idx],
row.names = rownames(pData(gset[[1]])))
# 整理列名
colnames(metadata) <- c('Study', 'Age', 'Distant.RFS', 'ER', 'GGI', 'Grade',
'Size', 'Time.RFS')
# 准备特定感兴趣的表型数据
metadata$Study <- gsub('Reanalyzed by: ', '', as.character(metadata$Study))
metadata$Age <- as.numeric(gsub('^KJ', NA, as.character(metadata$Age)))
metadata$Distant.RFS <- factor(metadata$Distant.RFS,
levels = c(0,1))
metadata$ER <- factor(gsub('\\?', NA, as.character(metadata$ER)),
levels = c(0,1))
metadata$ER <- factor(ifelse(metadata$ER == 1, 'ER+', 'ER-'),
levels = c('ER-', 'ER+'))
metadata$GGI <- as.numeric(as.character(metadata$GGI))
metadata$Grade <- factor(gsub('\\?', NA, as.character(metadata$Grade)),
levels = c(1,2,3))
metadata$Grade <- gsub(1, 'Grade 1', gsub(2, 'Grade 2', gsub(3, 'Grade 3', metadata$Grade)))
metadata$Grade <- factor(metadata$Grade, levels = c('Grade 1', 'Grade 2', 'Grade 3'))
metadata$Size <- as.numeric(as.character(metadata$Size))
metadata$Time.RFS <- as.numeric(gsub('^KJX|^KJ', NA, metadata$Time.RFS))
# 从表型数据中移除任何含有NA值的样本
discard <- apply(metadata, 1, function(x) any(is.na(x)))
metadata <- metadata[!discard,]
# 过滤表达数据以匹配我们表型数据中的样本
mat <- mat[,which(colnames(mat) %in% rownames(metadata))]
# 检查表型数据和表达数据之间的样本名称是否完全匹配
all(colnames(mat) == rownames(metadata))
# PCA分析
p <- pca(mat, metadata = metadata, removeVar = 0.1)
这段代码使用了R语言及其生物信息学包，来加载和处理来自Gene Expression Omnibus (GEO)的乳腺癌基因表达数据集。以下是逐步解释：
-
加载必要的库:
-
Biobase: 一个提供基础生物信息学数据结构和方法的R包。 -
GEOquery: 用于从NCBI的GEO数据库检索数据的R包。
-
-
从GEO加载数据:
- 使用
getGEO函数，通过指定的GEO系列编号GSE2990加载数据。GSEMatrix = TRUE表示希望以GEO表达数据矩阵的形式获取数据，getGPL = FALSE表示不需要加载平台数据（即芯片或测序平台的描述）。 -
exprs(gset[[1]])提取第一个数据集的表达矩阵。
- 使用
-
移除Affymetrix控制探针:
- 使用
grep函数找到行名（探针ID）以'AFFX'开头的行，并将这些行从表达矩阵mat中移除。这些探针通常用于技术控制和质量控制，不代表生物基因。
- 使用
-
从表型数据（phenotype data，简称pdata）中提取信息:
-
pData(gset[[1]])获取与样本相关的表型数据。 - 选择感兴趣的列，如年龄、远处复发自由生存（Distant RFS）等，并将这些信息保存在
metadata数据框中。
-
-
整理列名和数据:
- 修改
metadata的列名，使之更加清晰易读。 - 清洗和转换某些表型数据，如将年龄转换为数值型，将远处复发自由生存（Distant RFS）和雌激素受体状态（ER）转换为因子型等。
- 修改
-
移除有缺失值的样本:
- 使用
apply函数检查metadata中的每一行（即每个样本），如果任何行含有NA值，则标记为要丢弃。 - 从
metadata中移除这些含有缺失值的样本。
- 使用
-
过滤表达数据以匹配
metadata中的样本:- 将
mat矩阵中的列（代表样本）过滤，只保留那些在metadata中出现的样本。
- 将
-
检查样本名是否完全匹配:
- 确认经过过滤后的表达矩阵
mat的列名（样本名）与metadata的行名完全匹配，以确保表达数据和表型数据是对应的。
- 确认经过过滤后的表达矩阵
bi-plot
biplot(p, showLoadings = TRUE, lab = NULL)
探针205225_at指向下方，它靶向ESR1基因。这已经是一个有用的验证，因为部分由ESR1编码的雌激素受体在PC2（y轴）上有很强的表达，其受体状态从上到下由负变正。
pairs plot
pairsplot(p)
eigencor plot
特征相关图）用于展示主成分（特征值）与其他变量（例如，基因表达水平、表型特征等）之间的相关性。
eigencorplot(p,
metavars = c('Study','Age','Distant.RFS','ER',
'GGI','Grade','Size','Time.RFS'))
高级功能
加载数据
suppressMessages(require(hgu133a.db))
newnames <- mapIds(hgu133a.db,
keys = rownames(p$loadings),
column = c('SYMBOL'),
keytype = 'PROBEID')
# tidy up for NULL mappings and duplicated gene symbols
newnames <- ifelse(is.na(newnames) | duplicated(newnames),
names(newnames), newnames)
rownames(p$loadings) <- newnames
将原始数据中的探针ID转换为更易于理解和识别的基因名。
suppressMessages(require(hgu133a.db)): 这行代码加载hgu133a.db包，该包包含Affymetrix Human Genome U133A Array的注释数据。suppressMessages函数用于抑制加载包时可能产生的任何消息，使输出更加清洁。-
mapIds(hgu133a.db, ...):mapIds函数用于从hgu133a.db注释包中获取映射信息。该函数属于AnnotationDbi包，是处理生物信息学注释数据的通用工具。keys = rownames(p$loadings):keys参数接收要映射的ID，这里使用了p$loadings的行名，假设p$loadings包含了探针ID作为行名，这表示你可能在进行PCA分析后，想要将PCA加载（loadings）中的探针ID映射到基因符号。column = c('SYMBOL'): 指定返回的注释类型，这里是'SYMBOL'，即基因符号。keytype = 'PROBEID': 指明keys中提供的ID类型，这里是'PROBEID'，表示提供的是探针ID。
newnames: 映射结果被赋值给newnames变量，这个变量现在包含一个从探针ID到基因符号的映射，可以用于更新数据集，使其包含更直观的基因符号而不是原始的探针ID。-
处理NULL映射和重复的基因符号:
-
ifelse(is.na(newnames) | duplicated(newnames), names(newnames), newnames): 这行代码使用ifelse函数对每个newnames中的元素进行判断。如果元素是NA（表示没有找到对应的基因符号）或者是重复的基因符号，则使用names(newnames)（即原始探针ID）作为替代。否则，保持原始的基因符号映射。
-
-
更新PCA加载（loadings）的行名:
-
rownames(p$loadings) <- newnames: 将处理后的newnames（包含了基因符号或在必要时回退到探针ID）设置为p$loadings的行名。这一步骤是在将映射结果应用到PCA结果上，将PCA加载矩阵中的探针ID替换为对应的基因符号或在无法映射时保留探针ID。
-
通过这样的处理，即使在遇到无法映射的探针ID或基因符号重复的情况下，也能保持数据的一致性和清晰度。这对于后续的数据分析和解释非常重要，特别是在需要将PCA结果与生物学意义联系起来时。
确定保留的最佳主成分数量
horn <- parallelPCA(mat)
horn$n
elbow <- findElbowPoint(p$variance)
elbow
library(ggplot2)
screeplot(p,
components = getComponents(p, 1:20),
vline = c(horn$n, elbow)) +
geom_label(aes(x = horn$n + 1, y = 50,
label = 'Horn\'s', vjust = -1, size = 8)) +
geom_label(aes(x = elbow + 1, y = 50,
label = 'Elbow method', vjust = -1, size = 8))
比较pc1与pc2
biplot(p,
lab = paste0(p$metadata$Age, ' años'),
colby = 'ER',
hline = 0, vline = 0,
legendPosition = 'right')
自定义颜色并突出显示
biplot(p,
colby = 'ER', colkey = c('ER+' = 'forestgreen', 'ER-' = 'purple'),
colLegendTitle = 'ER-\nstatus',
# encircle config
encircle = TRUE,
encircleFill = TRUE,
hline = 0, vline = c(-25, 0, 25),
legendPosition = 'top', legendLabSize = 16, legendIconSize = 8.0)
biplot(p,
colby = 'ER', colkey = c('ER+' = 'forestgreen', 'ER-' = 'purple'),
colLegendTitle = 'ER-\nstatus',
# encircle config
encircle = TRUE, encircleFill = FALSE,
encircleAlpha = 1, encircleLineSize = 5,
hline = 0, vline = c(-25, 0, 25),
legendPosition = 'top', legendLabSize = 16, legendIconSize = 8.0)
加95%置信水平的椭圆
biplot(p,
colby = 'ER', colkey = c('ER+' = 'forestgreen', 'ER-' = 'purple'),
# ellipse config
ellipse = TRUE,
ellipseType = 't',
ellipseLevel = 0.95,
ellipseFill = TRUE,
ellipseAlpha = 1/4,
ellipseLineSize = 1.0,
xlim = c(-125,125), ylim = c(-50, 80),
hline = 0, vline = c(-25, 0, 25),
legendPosition = 'top', legendLabSize = 16, legendIconSize = 8.0)
修改线型、移除网格线并增大点的大小
biplot(p,
lab = NULL,
colby = 'ER', colkey = c('ER+'='royalblue', 'ER-'='red3'),
hline = c(-25, 0, 25), vline = c(-25, 0, 25),
vlineType = c('dotdash', 'solid', 'dashed'),
gridlines.major = FALSE, gridlines.minor = FALSE,
pointSize = 5,
legendPosition = 'left', legendLabSize = 14, legendIconSize = 8.0,
shape = 'Grade', shapekey = c('Grade 1'=15, 'Grade 2'=17, 'Grade 3'=8),
drawConnectors = FALSE,
title = 'PCA bi-plot',
subtitle = 'PC1 versus PC2',
caption = '27 PCs ≈ 80%')
根据连续变量进行着色并绘制其他主成分
# add ESR1 gene expression to the metadata
p$metadata$ESR1 <- mat['205225_at',]
biplot(p,
x = 'PC2', y = 'PC3',
lab = NULL,
colby = 'ESR1',
shape = 'ER',
hline = 0, vline = 0,
legendPosition = 'right') +
scale_colour_gradient(low = 'gold', high = 'red2')
通过成对图快速探索可能具有信息量的主成分
pairsplot(p,
components = getComponents(p, c(1:10)),
triangle = TRUE, trianglelabSize = 12,
hline = 0, vline = 0,
pointSize = 0.4,
gridlines.major = FALSE, gridlines.minor = FALSE,
colby = 'Grade',
title = 'Pairs plot', plotaxes = FALSE,
margingaps = unit(c(-0.01, -0.01, -0.01, -0.01), 'cm'))
pairsplot(p,
components = getComponents(p, c(4,33,11,1)),
triangle = FALSE,
hline = 0, vline = 0,
pointSize = 0.8,
gridlines.major = FALSE, gridlines.minor = FALSE,
colby = 'ER',
title = 'Pairs plot', titleLabSize = 22,
axisLabSize = 14, plotaxes = TRUE,
margingaps = unit(c(0.1, 0.1, 0.1, 0.1), 'cm'))
组合绘图
pscree <- screeplot(p, components = getComponents(p, 1:30),
hline = 80, vline = 27, axisLabSize = 14, titleLabSize = 20,
returnPlot = FALSE) +
geom_label(aes(20, 80, label = '80% explained variation', vjust = -1, size = 8))
ppairs <- pairsplot(p, components = getComponents(p, c(1:3)),
triangle = TRUE, trianglelabSize = 12,
hline = 0, vline = 0,
pointSize = 0.8, gridlines.major = FALSE, gridlines.minor = FALSE,
colby = 'Grade',
title = '', plotaxes = FALSE,
margingaps = unit(c(0.01, 0.01, 0.01, 0.01), 'cm'),
returnPlot = FALSE)
pbiplot <- biplot(p,
# loadings parameters
showLoadings = TRUE,
lengthLoadingsArrowsFactor = 1.5,
sizeLoadingsNames = 4,
colLoadingsNames = 'red4',
# other parameters
lab = NULL,
colby = 'ER', colkey = c('ER+'='royalblue', 'ER-'='red3'),
hline = 0, vline = c(-25, 0, 25),
vlineType = c('dotdash', 'solid', 'dashed'),
gridlines.major = FALSE, gridlines.minor = FALSE,
pointSize = 5,
legendPosition = 'none', legendLabSize = 16, legendIconSize = 8.0,
shape = 'Grade', shapekey = c('Grade 1'=15, 'Grade 2'=17, 'Grade 3'=8),
drawConnectors = FALSE,
title = 'PCA bi-plot',
subtitle = 'PC1 versus PC2',
caption = '27 PCs ≈ 80%',
returnPlot = FALSE)
ploadings <- plotloadings(p, rangeRetain = 0.01, labSize = 4,
title = 'Loadings plot', axisLabSize = 12,
subtitle = 'PC1, PC2, PC3, PC4, PC5',
caption = 'Top 1% variables',
shape = 24, shapeSizeRange = c(4, 8),
col = c('limegreen', 'black', 'red3'),
legendPosition = 'none',
drawConnectors = FALSE,
returnPlot = FALSE)
peigencor <- eigencorplot(p,
components = getComponents(p, 1:10),
metavars = c('Study','Age','Distant.RFS','ER',
'GGI','Grade','Size','Time.RFS'),
cexCorval = 1.0,
fontCorval = 2,
posLab = 'all',
rotLabX = 45,
scale = TRUE,
main = "PC clinical correlates",
cexMain = 1.5,
plotRsquared = FALSE,
corFUN = 'pearson',
corUSE = 'pairwise.complete.obs',
signifSymbols = c('****', '***', '**', '*', ''),
signifCutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1),
returnPlot = FALSE)
library(cowplot)
library(ggplotify)
top_row <- plot_grid(pscree, ppairs, pbiplot,
ncol = 3,
labels = c('A', 'B Pairs plot', 'C'),
label_fontfamily = 'serif',
label_fontface = 'bold',
label_size = 22,
align = 'h',
rel_widths = c(1.10, 0.80, 1.10))
bottom_row <- plot_grid(ploadings,
as.grob(peigencor),
ncol = 2,
labels = c('D', 'E'),
label_fontfamily = 'serif',
label_fontface = 'bold',
label_size = 22,
align = 'h',
rel_widths = c(0.8, 1.2))
plot_grid(top_row, bottom_row, ncol = 1,
rel_heights = c(1.1, 0.9))
以上内容来自PCAtools帮助文档
