主成分分析（PCA）是一種非常強大的技術暖混，廣泛應用于數據科學、生物信息學和其他領域晾咪。它最初是為了分析大量數據而開發(fā)的，以便揭示被分析的邏輯實體之間的差異/關系塞赂。它提取數據的基本結構昼蛀，而不需要構建任何模型來表示它。通過一種簡化過程仇哆，可以將大量變量轉換成較少的不相關變量（即“主成分”）夫植，同時又能夠便于對原始數據進行解釋，從而得到數據的主要信息辟拷。

安裝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”數據集，這個數據集包含了不同的氣道平滑肌細胞在被地塞米松處理后的信息隅俘。

對數據進行歸一化笤喳，并將歸一化后的計數轉換為方差穩(wěn)定的表達水平：

  library('DESeq2')

#構建了一個考慮了細胞類型和處理條件的差異表達分析的數據集對象
  dds <- DESeqDataSet(airway, design = ~ cell + dex)
#執(zhí)行差異表達分析
  dds <- DESeq(dds)
#用vst函數對dds對象進行方差穩(wěn)定轉換
  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)的乳腺癌基因表達數據集魄缚。以下是逐步解釋：

1. 加載必要的庫:

   • Biobase: 一個提供基礎生物信息學數據結構和方法的R包焚廊。
   • GEOquery: 用于從NCBI的GEO數據庫檢索數據的R包。

2. 從GEO加載數據:

   • 使用getGEO函數徙硅，通過指定的GEO系列編號GSE2990加載數據搞疗。GSEMatrix = TRUE表示希望以GEO表達數據矩陣的形式獲取數據须肆，getGPL = FALSE表示不需要加載平臺數據（即芯片或測序平臺的描述）桩皿。
   • exprs(gset[[1]])提取第一個數據集的表達矩陣。

3. 移除Affymetrix控制探針:

   • 使用grep函數找到行名（探針I(yè)D）以'AFFX'開頭的行拒贱，并將這些行從表達矩陣mat中移除佛嬉。這些探針通常用于技術控制和質量控制，不代表生物基因斜做。

4. 從表型數據（phenotype data湾揽，簡稱pdata）中提取信息:

   • pData(gset[[1]])獲取與樣本相關的表型數據。
   • 選擇感興趣的列霸旗，如年齡戚揭、遠處復發(fā)自由生存（Distant RFS）等，并將這些信息保存在metadata數據框中蔬啡。

5. 整理列名和數據:

   • 修改metadata的列名镀虐，使之更加清晰易讀。
   • 清洗和轉換某些表型數據空猜，如將年齡轉換為數值型恨旱，將遠處復發(fā)自由生存（Distant RFS）和雌激素受體狀態(tài)（ER）轉換為因子型等。

6. 移除有缺失值的樣本:

   • 使用apply函數檢查metadata中的每一行（即每個樣本）谆沃，如果任何行含有NA值仪芒，則標記為要丟棄耕陷。
   • 從metadata中移除這些含有缺失值的樣本哟沫。

7. 過濾表達數據以匹配metadata中的樣本:

   • 將mat矩陣中的列（代表樣本）過濾锌介，只保留那些在metadata中出現的樣本。

8. 檢查樣本名是否完全匹配:

   • 確認經過過濾后的表達矩陣mat的列名（樣本名）與metadata的行名完全匹配隆敢，以確保表達數據和表型數據是對應的崔慧。

bi-plot

  biplot(p, showLoadings = TRUE, lab = NULL)

探針205225_at指向下方，它靶向ESR1基因。這已經是一個有用的驗證拇涤，因為部分由ESR1編碼的雌激素受體在PC2（y軸）上有很強的表達誉结，其受體狀態(tài)從上到下由負變正。

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

將原始數據中的探針I(yè)D轉換為更易于理解和識別的基因名。

1. suppressMessages(require(hgu133a.db)): 這行代碼加載hgu133a.db包滥沫，該包包含Affymetrix Human Genome U133A Array的注釋數據。suppressMessages函數用于抑制加載包時可能產生的任何消息世分，使輸出更加清潔缀辩。

2. mapIds(hgu133a.db, ...): mapIds函數用于從hgu133a.db注釋包中獲取映射信息踪央。該函數屬于AnnotationDbi包杯瞻，是處理生物信息學注釋數據的通用工具炫掐。

   • keys = rownames(p$loadings): keys參數接收要映射的ID，這里使用了p$loadings的行名募胃，假設p$loadings包含了探針I(yè)D作為行名痹束，這表示你可能在進行PCA分析后，想要將PCA加載（loadings）中的探針I(yè)D映射到基因符號祷嘶。

   • column = c('SYMBOL'): 指定返回的注釋類型论巍，這里是'SYMBOL'，即基因符號嘉汰。

   • keytype = 'PROBEID': 指明keys中提供的ID類型鞋怀，這里是'PROBEID'，表示提供的是探針I(yè)D密似。

3. newnames: 映射結果被賦值給newnames變量，這個變量現在包含一個從探針I(yè)D到基因符號的映射薄扁，可以用于更新數據集废累，使其包含更直觀的基因符號而不是原始的探針I(yè)D。

4. 處理NULL映射和重復的基因符號:

   • ifelse(is.na(newnames) | duplicated(newnames), names(newnames), newnames): 這行代碼使用ifelse函數對每個newnames中的元素進行判斷邑滨。如果元素是NA（表示沒有找到對應的基因符號）或者是重復的基因符號，則使用names(newnames)（即原始探針I(yè)D）作為替代匣距。否則，保持原始的基因符號映射尚卫。

5. 更新PCA加載（loadings）的行名:
   • rownames(p$loadings) <- newnames: 將處理后的newnames（包含了基因符號或在必要時回退到探針I(yè)D）設置為p$loadings的行名尸红。這一步驟是在將映射結果應用到PCA結果上，將PCA加載矩陣中的探針I(yè)D替換為對應的基因符號或在無法映射時保留探針I(yè)D怎爵。

通過這樣的處理盅蝗，即使在遇到無法映射的探針I(yè)D或基因符號重復的情況下，也能保持數據的一致性和清晰度芙委。這對于后續(xù)的數據分析和解釋非常重要狂秦，特別是在需要將PCA結果與生物學意義聯(lián)系起來時。

確定保留的最佳主成分數量

  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%')

根據連續(xù)變量進行著色并繪制其他主成分

  # 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幫助文檔