前言
本教學文件旨在詳細解說一份用於分析人類間質幹細胞 (Mesenchymal Stem Cells, MSCs) 在 2D 與 3D 培養環境下基因表現差異的 R 腳本。我們將使用 Affymetrix HTA 2.0 微陣列平台的數據，並透過一系列生物資訊學工具，從原始數據讀取、標準化、註解，到最終的數據可視化，一步步完成整個分析流程。
本教學適合對象：
- 對微陣列數據分析有興趣的生物學家或研究人員。
- 正在學習 R 語言進行生物資訊分析的學生。
- 希望了解
oligo和Seurat套件在微陣列分析中應用的使用者。
1. 原始數據讀取與標準化 (Oligo 套件)
在分析的第一步，我們需要處理最原始的 .CEL 檔案。.CEL 檔案是 Affymetrix 微陣列掃描儀產生的原始數據檔案，包含了每個探針的螢光強度值。我們將使用 oligo 這個專為處理 Affymetrix 晶片數據而設計的 R 套件來讀取並進行標準化。
1-1. 讀取所有樣本的 .CEL 檔案
這個區塊的目標是：
- 載入所有需要的 R 套件。
- 定義一個函數，確保
.CEL檔案能按照檔案名稱中的數字順序被正確讀取。 - 找到指定路徑下的所有
.CEL檔案。 - 使用
oligo套件的read.celfiles函數將這些檔案讀取成一個ExpressionFeatureSet物件。 - 將讀取後的原始數據物件儲存為
.rds檔案，方便未來快速取用，無需重複讀取。
1r 1-1-load-raw-data, eval=FALSE}
2# 清空當前工作環境中的所有變數，確保一個乾淨的開始
3rm(list=ls(all=TRUE))
4
5# 載入本次分析所需的核心套件
6library(affy) # 傳統 Affymetrix 分析工具
7library(oligo) # 新一代 Affymetrix 分析工具，支援較新的晶片
8library(oligoData) # oligo 套件所需的範例數據
9library(AnnotationDbi) # 基因註解的核心工具
10library(affycoretools) # Affy 分析的輔助工具
11library(Biobase) # 生物資訊分析的基礎套件，定義了 eSet 等核心物件
12library(pd.hta.2.0) # HTA 2.0 晶片專用的平台設計文件 (Platform Design)
13
14# --- 自定義函數：按檔名中的數字排序 ---
15# 這個函數的目的是為了解決 list.celfiles() 可能不會按數字順序回傳檔名的問題
16# 例如，它可能會回傳 "1.cel", "10.cel", "2.cel" 而不是 "1.cel", "2.cel", "10.cel"
17numeric_order <- function(filepaths) {
18 # 從完整路徑中提取檔案名稱
19 filenames <- basename(filepaths)
20 # 使用正規表示式從檔名開頭提取數字
21 numbers <- sapply(filenames, function(filename) {
22 matches <- regexpr("^[0-9]+", filename)
23 if (matches[1] == -1) {
24 # 如果檔名不是以數字開頭，則回傳 NA
25 NA_integer_
26 } else {
27 # 將提取出的字串轉為整數
28 as.integer(substr(filename, matches[1], matches[1] + attr(matches, "match.length") - 1))
29 }
30 })
31 # 回傳根據提取出的數字排序後的索引
32 order(numbers, na.last = TRUE)
33}
34
35# 設定 .CEL 檔案所在的資料夾路徑，並列出所有 .CEL 檔案
36cel_files <- list.celfiles("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/1. Original data/CEL", full.names=T)
37
38# 使用我們自定義的函數來對檔案列表進行排序
39cel_files <- cel_files[numeric_order(cel_files)]
40
41# 使用 oligo 套件的 read.celfiles 函數讀取所有 .CEL 檔
42GPL17586_all_sample_raw <- read.celfiles(cel_files)
43
44# 設定儲存 .rds 檔案的路徑
45setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS")
46
47# 將讀取的原始數據物件儲存起來，方便下次直接載入
48saveRDS(GPL17586_all_sample_raw, file = "0. GPL17586_15_sample_raw_250625.rds")
1-2. 數據標準化與探針註解
讀取原始數據後，下一步是「標準化」(Normalization)。標準化的目的是消除不同晶片之間由於實驗技術差異（如螢光標記效率、掃描儀敏感度等）造成的系統性誤差，使得樣本之間的基因表現量具有可比性。這裡我們使用最廣泛的 RMA (Robust Multi-array Average) 演算法。
完成標準化後，我們會得到每個探針組 (probeset) 的表現量，但我們通常更關心的是基因層面的表現。因此，我們需要對這些探針組進行「註解」(Annotation)，也就是找出它們對應的基因是什麼。
1r 1-2-rma-normalization, eval=FALSE}
2# 清空環境
3rm(list=ls(all=TRUE))
4
5# 載入所需套件
6library(oligo)
7library(oligoData)
8library(pd.hta.2.0)
9library(AnnotationDbi)
10library(affycoretools)
11library(Biobase)
12
13# 設定工作目錄到 .rds 檔案所在位置
14setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS")
15
16# 讀取上一部儲存的原始數據物件
17raw_data <- readRDS("0. GPL17586_15_sample_raw_250625.rds", refhook = NULL)
18
19# 使用 RMA 演算法進行標準化
20# RMA 包含三個步驟：背景校正、Quantile Normalization、探針匯總
21rma_data <- rma(raw_data)
22
23# 使用 annotateEset 函數為標準化後的數據添加註解
24# pd.hta.2.0 包含了 HTA 2.0 晶片上每個探針的詳細資訊
25rma_data_anno1 <- annotateEset(rma_data, pd.hta.2.0)
26
27# 從註解後的物件中，分別提取特徵資料 (featureData) 和表現量矩陣 (assayData)
28x1 <- as.matrix(rma_data_anno1@featureData@data) # 包含了探針ID、對應的基因符號等
29x2 <- as.matrix(rma_data_anno1@assayData[["exprs"]]) # 標準化後的表現量 (log2尺度)
30
31# 將特徵資料和表現量矩陣合併成一個大的表格
32RMA_GPL17586_All_sample <- cbind(x1, x2)
33
34# 移除暫存變數
35rm(x1, x2)
36
37# 設定輸出路徑
38setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS/1. Extend file")
39
40# 將合併後的數據寫入一個文字檔，方便查看或給其他軟體使用
41write.table(RMA_GPL17586_All_sample, "1. GPL17586_15_sample_RMA_250625.txt", sep="\t", row.names=T, col.names=T, quote=F)
1-3. 處理多探針映射問題
微陣列晶片上，常常會有多個不同的探針組 (probeset) 對應到同一個基因。在進行下游分析前，我們通常需要將這些探針的表現量整合成一個單一的數值來代表該基因的表現。
常見的策略有：
- 取平均值。
- 取中位數。
- 選擇表現量最高的那個探針。
- 選擇與同基因其他探針相關性最高的那個。
在這個腳本中，我們採用的策略是：對於有多個探針對應的基因，我們選擇在樣本間表現量差異絕對值最大的那個探針。這背後的假設是，差異越大的探針可能越能反映真實的生物學變化。
1r 1-3-handle-multi-probes, eval=FALSE}
2# 清空環境
3rm(list=ls(all=TRUE))
4
5# 載入所需套件
6library(dplyr)
7library(Seurat) # 雖然主要用於單細胞，但其強大的數據結構和工具也可用於微陣列
8# ... (其他套件載入)
9library(readr)
10
11# 設定工作目錄
12setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS/1. Extend file")
13
14# 讀取上一步產生的 RMA 標準化後的數據
15da0 <- read.table('1. GPL17586_15_sample_RMA_250625.txt', row.names = 1, sep = '\t', check.names = FALSE)
16
17# --- 準備數據以篩選探針 ---
18# 複製 da0，並增加兩欄用於計算
19da0_1 <- da0[,c(1:length(colnames(da0)), 1, 1)]
20names(da0_1)[length(colnames(da0))+1] <- "Value"
21names(da0_1)[length(colnames(da0))+2] <- "Abs Value"
22
23# 計算第6個樣本和第5個樣本表現量的差異 (這裡的欄位選擇可能需要根據實際實驗設計來決定，此處為範例)
24da0_1[,length(colnames(da0))+1] <- da0_1[,6] - da0_1[,5]
25# 計算該差異的絕對值
26da0_1[,length(colnames(da0))+2] <- abs(da0_1[,length(colnames(da0))+1])
27
28# 移除那些沒有對應到任何基因符號 (SYMBOL) 的探針
29da0_1 <- da0_1[!is.na(da0_1$SYMBOL), ]
30
31# --- 區分單一映射和多重映射的基因 ---
32da <- da0_1
33# 計算每個基因符號 (SYMBOL) 出現的次數
34pps1 <- table(da[,3]) # 第3欄是 SYMBOL
35# 找出只出現一次的基因 (單一探針)
36pp1 <- names(pps1[pps1==1])
37# 找出出現超過一次的基因 (多個探針)
38pp2 <- names(pps1[pps1>1])
39
40# 將單一探針的數據直接存入 box1
41box1 <- da[da[,3] %in% pp1,]
42
43# --- 迴圈處理多重映射的基因 ---
44# 預先建立一個空的 dataframe 來存放篩選結果
45box2 <- da[1:length(pp2),]
46# 紀錄開始時間
47t0 <- Sys.time()
48# 迴圈遍歷每一個有多個探針的基因
49for (i in 1:length(pp2)) {
50 # 找到這個基因對應的所有探針
51 xx <- (da[,3] == pp2[i])
52 qaz <- da[xx,]
53 # 根據我們之前計算的 "Abs Value" 欄位，找到值最大的那一個探針的索引
54 qa1 <- order(abs(as.numeric(qaz[,ncol(qaz)])), decreasing = T)[1]
55 # 將這個探針的數據存入 box2
56 box2[i,] <- qaz[qa1,]
57}
58# 紀錄結束時間並計算耗時
59tn <- Sys.time()
60difftime(tn, t0, tz=" ", unit='mins')
61
62# --- 合併結果 ---
63# 將單一探針的結果和多探針篩選後的結果合併
64box3 <- rbind(box1, box2)
65# 移除我們為了計算而增加的 "Value" 和 "Abs Value" 兩欄
66box3 <- box3[,c(-(length(colnames(da0))+1), -(length(colnames(da0))+2))]
67
68# 設定輸出路徑
69setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS/1. Extend file")
70# 將最終整理好的基因表現矩陣寫入檔案
71write.table(box3, "2. Muti-symbols_largest_15_sample_250625.txt", sep="\t", row.names=F, col.names=T, quote=F)
2. 創建 Seurat 物件並進行下游分析
Seurat 是當前單細胞 RNA 定序 (scRNA-seq) 分析領域最流行的 R 套件。儘管我們的數據來自微陣列，但我們仍然可以借用 Seurat 強大的數據結構和分析功能。將數據轉換為 Seurat 物件可以讓我們方便地進行正規化、降維 (PCA, UMAP)、分群以及可視化。
2-1. 創建 Seurat 物件
這個步驟的目標是：
- 讀取上一步產生的、已經處理好多探針問題的基因表現矩陣。
- 讀取包含樣本資訊的元數據 (metadata) 檔案。
- 對數據進行整理，使其符合
CreateSeuratObject的輸入格式。 - 特別注意：微陣列的 RMA 數據是 log2 尺度的，而
CreateSeuratObject通常預期輸入原始計數 (raw counts)。因此，我們需要先用2^da0將數據轉換回線性尺度。 - 執行
Seurat的標準分析流程：正規化、尋找高變異基因、數據縮放、PCA。 - 最後，將樣本的元數據添加到
Seurat物件中，並儲存物件。
1r 2-1-create-seurat-object, eval=FALSE}
2# 清空環境
3rm(list=ls(all=T))
4
5# 載入所需套件
6library(dplyr)
7library(Seurat)
8library(ggplot2)
9# ... (其他套件)
10library(readxl) # 用於讀取 Excel 檔案
11
12# 設定工作目錄
13setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS/1. Extend file")
14
15# 讀取基因表現矩陣 (注意這裡讀取的是 Excel 檔)
16da0 <- read_excel("3. Muti-symbols_largest_15_sample_250625.escaped.xlsx", sheet = "1. Original RMA")
17# 讀取樣本元數據
18meta_info <- read_excel("0. 20250625 MSC info.xlsx", sheet = "2. Edited_1")
19
20# --- 數據整理 ---
21# 將數據轉為矩陣格式
22da0 <- as.matrix(da0)
23# 使用基因符號 (SYMBOL) 作為行名 (rownames)
24rownames(da0) <- da0[,3]
25# 移除前面幾欄的註解資訊，只保留樣本的表現量數據
26da0 <- da0[,c(-1:-4)]
27# 轉回 data frame 格式
28da0 <- as.data.frame(da0)
29# 將欄名 (colnames) 設置為元數據中對應的樣本名
30colnames(da0) <- meta_info$`Sample Name`
31
32# 將所有表現量數據從字串轉為數值格式
33da0 <- da0 %>% mutate(across(everything(), as.numeric))
34
35# !!! 重要步驟 !!!
36# RMA 數據是 log2 尺度的，Seurat 的標準流程預期是 counts，所以我們先把它轉回線性尺度
37da0 <- 2^da0
38
39# 使用整理好的數據創建 Seurat 物件
40da0 <- CreateSeuratObject(counts = da0)
41
42# --- 執行 Seurat 標準分析流程 ---
43# 官方教學: https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
44
45# 正規化數據 (這裡使用的是 LogNormalize)
46da0 <- NormalizeData(da0)
47# 尋找樣本間高變異的基因
48da0 <- FindVariableFeatures(da0)
49# 縮放數據，使其平均值為0，變異數為1，這是 PCA 的前提
50da0 <- ScaleData(da0)
51
52# 執行主成分分析 (PCA)
53da0 <- RunPCA(da0, features = VariableFeatures(da0), npcs = 14) # npcs 設為 14
54# 繪製 PCA 圖，觀察樣本分佈
55DimPlot(da0, reduction = "pca")
56# 繪製 "Elbow Plot"，幫助判斷要選擇多少個主成分 (PC) 進行下游分析
57ElbowPlot(da0)
58
59# 根據選擇的 PC 數 (1到14) 來建構 KNN 圖
60da0 <- FindNeighbors(da0, dims = 1:14)
61# 根據 KNN 圖進行分群
62da0 <- FindClusters(da0, resolution = 0.3, verbose = FALSE)
63# 執行 UMAP 降維分析
64da0 <- RunUMAP(da0, dims = 1:14, n.neighbors = 5)
65# 繪製 UMAP 圖，並標示分群結果
66DimPlot(da0, label = TRUE)
67
68# --- 添加元數據 ---
69# 重新讀取元數據 (確保萬無一失)
70setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS/1. Extend file")
71meta_info <- read_excel("0. 20250625 MSC info.xlsx", sheet = "2. Edited_1")
72# 將元數據的行名設置為 Seurat 物件中細胞 (樣本) 的名稱，以確保能正確匹配
73rownames(meta_info) <- rownames(da0@meta.data)
74# 將我們的元數據和 Seurat 物件原有的元數據合併
75meta_info_cbind <- cbind(meta_info, da0@meta.data)
76# 更新 Seurat 物件的元數據
77da0@meta.data <- meta_info_cbind
78
79# 設定儲存路徑
80setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS/")
81# 將處理好的 Seurat 物件儲存起來
82saveRDS(da0, file = "1. GPL17586_15_sample_250625.rds")
3. 數據可視化
數據可視化是探索性數據分析的關鍵一步。透過圖形，我們可以直觀地了解樣本之間的關係、實驗分組的效果以及特定基因的表現模式。
3-1. UMAP 可視化
UMAP (Uniform Manifold Approximation and Projection) 是一種強大的非線性降維技術，能將高維度的數據投影到二維或三維空間中，同時盡可能地保持數據的全局結構。在圖中，相似的樣本會被聚集在一起。
這個區塊展示了如何：
- 根據不同的實驗分組 (
Group_3,Group_1,Sample Name) 來對 UMAP 圖上的點進行著色或塑形。 - 客製化
ggplot2圖形，例如調整圖例、標題等。 - 創建互動式的 UMAP 圖 (
plotly)，讓你可以用滑鼠懸停在點上查看樣本資訊。 - 使用
ggrepel套件為每個點加上文字標籤，並防止標籤重疊。
1r 3-1-umap-visualization, eval=FALSE}
2# --- 繪製基礎 UMAP 圖 ---
3# ... (前面省略了清空環境和載入套件的程式碼) ...
4setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS/")
5da0 <- readRDS("1. GPL17586_15_sample_250625.rds", refhook = NULL)
6
7# 將樣本的 "身份" (Identity) 設置為 Group_3
8Idents(da0) <- factor(da0$Group_3)
9# 根據 Group_3 著色，並根據 Group_1 分割成多個子圖
10p1 <- DimPlot(da0, reduction = 'umap', pt.size = 4, label.size = 4,
11 group.by = "Group_3", label = F, split.by = "Group_1", ncol = 2) + NoLegend()
12p1
13
14# --- 繪製更複雜、可發表的 UMAP 圖 ---
15# 這次我們將 Group_1 作為顏色，Group_3 (ident) 作為形狀
16Idents(da0) <- factor(da0$Group_3)
17p2 <- DimPlot(da0, reduction = 'umap', pt.size = 4, label.size = 4, label = F) + theme(legend.position = "right")
18p2 <- p2 + guides(color = guide_legend(override.aes = list(size=4), ncol=1))
19p2[[1]][["labels"]][["title"]] <- "" # 移除圖例標題
20p2
21
22# --- 使用 ggplot2 從頭繪製 UMAP ---
23# 有時候 Seurat 的 DimPlot 無法滿足所有的客製化需求，
24# 這時我們可以提取出 UMAP 的座標和元數據，自己用 ggplot2 來畫。
25
26# 提取 UMAP 座標數據
27p2_data <- p2$data
28# 提取元數據
29metadata <- da0@meta.data
30# 將 UMAP 座標和我們感興趣的元數據 (Group_1, Sample Name) 合併
31p2_data <- cbind(p2_data, metadata$Group_1, metadata$`Sample Name`)
32
33# 使用 ggplot2 繪圖
34p2_redraw <- ggplot(data = p2_data, aes(x = umap_1, y = umap_2)) +
35 # 點的顏色由 Group_1 決定，形狀由 ident (即 Group_3) 決定
36 geom_point(aes(color = p2_data$`metadata$Group_1`, shape = p2_data$ident), size = 3) +
37 # 設定主題，美化圖形，例如移除背景格線，設定背景為透明
38 theme(
39 panel.border = element_blank(),
40 panel.grid.major = element_blank(),
41 panel.grid.minor = element_blank(),
42 axis.line.x = element_line(color = "black"),
43 axis.line.y = element_line(color = "black"),
44 panel.background = element_rect(fill = "transparent", colour = NA),
45 # ... (其他美化設定) ...
46 ) +
47 # 設定座標軸和圖例的標題
48 labs(x = "UMAP_1", y = paste0("UMAP_2"), color = '', shape = '')
49# 調整圖例中點的大小和排列順序
50p2_redraw <- p2_redraw + guides(color = guide_legend(override.aes = list(size=4), ncol=1, order = 2),
51 shape = guide_legend(override.aes = list(size=4), ncol=1, order = 1))
52p2_redraw
53
54# --- 添加互動式標籤 ---
55# 使用 ggrepel 為每個點加上樣本名稱標籤
56library(ggrepel)
57p2_redraw <- p2_redraw +
58 geom_text_repel(aes(label = p2_data$`metadata$`Sample Name``), size = 3, show.legend = FALSE,
59 box.padding = unit(0.5, 'lines'))
60p2_redraw
61
62# --- 轉換為互動式圖形 ---
63library(plotly)
64p2_interactive <- ggplotly(p2_redraw)
65p2_interactive
3-2. 篩選樣本並重新分析
有時我們需要專注於一部分樣本進行更深入的分析。Seurat 的 subset 函數可以很方便地實現這一點。
1r 3-2-subset-samples, eval=FALSE}
2# ... (省略載入套件) ...
3setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS/")
4da0 <- readRDS("1. GPL17586_15_sample_250625.rds", refhook = NULL)
5
6# 讀取一個只包含我們想保留的樣本資訊的 Excel 檔案
7Meta_info_1 <- read_excel("1. Extend file/0. 20250625 MSC info.xlsx", sheet = "3. Edited_2")
8
9# 使用 subset 函數，根據樣本名稱 (cells) 篩選 Seurat 物件
10da1 <- subset(da0, cells = Meta_info_1$`Sample Name`)
11
12# 將篩選後的 Seurat 物件另存為一個新的 .rds 檔案
13setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS/")
14saveRDS(da1, file = "2. GPL17586_12_sample_250625.rds")
15
16# 接下來，可以對這個新的 da1 物件重複 3-1 中的 UMAP 繪圖步驟，
17# 觀察只包含這12個樣本時的數據分佈。
18# (此處省略了重複的繪圖程式碼)
3-3. 特定基因表現量的盒狀圖 (Barplot)
當我們發現一些有趣的基因時，通常會想看看它們在不同樣本或不同組別中的表現量分佈。盒狀圖 (Boxplot) 或小提琴圖 (Violin Plot) 是展示這類資訊的絕佳方式。
這個區塊的程式碼會：
- 遍歷一個預設的基因列表 (例如
c("IL1A", "IL1B"))。 - 為列表中的每一個基因，從
Seurat物件中提取其在所有樣本中的表現量。 - 使用
ggpubr套件中的ggboxplot函數繪製盒狀圖。 - 將每個基因的圖儲存為一個獨立的 .png 檔案。
1r 3-3-boxplots, eval=FALSE}
2# ... (省略載入套件) ...
3setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS/")
4# 這次我們讀取只包含12個樣本的 Seurat 物件
5da0 <- readRDS("2. GPL17586_12_sample_250625.rds", refhook = NULL)
6
7# 確保我們使用的是 RNA assay 的數據
8DefaultAssay(da0) <- "RNA"
9# 提取元數據
10metainfo <- da0@meta.data
11
12# 從 Seurat 物件中提取出基因表現矩陣
13# 注意：這裡提取的是 `counts` 層，即 2^RMA 的值
14da1 <- t(as.data.frame(da0@assays[["RNA"]]@layers[["counts"]]))
15da1 <- as.data.frame(da1)
16rownames(da1) <- colnames(da0) # 行名是樣本
17colnames(da1) <- rownames(da0) # 欄名是基因
18
19# 將樣本名稱作為新的一欄 "CellType" 加入矩陣
20da1 <- cbind(da1, metainfo$`Sample Name`)
21colnames(da1)[ncol(da1)] <- "CellType"
22# 將 CellType 轉為因子 (factor)，並指定其順序，確保 X 軸按我們想要的順序排列
23ordered_names <- unique(metainfo$`Sample Name`)
24da1$CellType <- factor(da1$CellType, levels = ordered_names)
25
26# 定義我們感興趣的基因列表
27my_gene_list <- intersect(colnames(da1), c("IL1A", "IL1B"))
28
29# 建立一個空的列表來存放畫好的圖
30my_plot_list <- list()
31
32# 迴圈遍歷基因列表
33for(i in 1:length(my_gene_list)) {
34 # 提取當前基因和分組資訊
35 new <- da1[, c(my_gene_list[i], "CellType")]
36 colnames(new) <- c("interest", "Group")
37
38 # 使用 ggboxplot 繪製盒狀圖
39 bxp <- ggboxplot(new, x = "Group", y = "interest",
40 color = "Group",
41 add = "jitter", legend = "none") # add = "jitter" 會加上散點
42 # 美化圖形，加上標題、座標軸標籤等
43 bxp <- bxp + ggtitle(my_gene_list[i]) +
44 xlab("") + ylab("Expression Level (Linear Scale)") +
45 theme(plot.title = element_text(size = 16, face = "bold"),
46 axis.text.x = element_text(angle = 45, hjust = 1)) # X軸標籤旋轉45度
47
48 # 將畫好的圖存入列表
49 my_plot_list[[i]] <- bxp
50
51 # 設定儲存路徑並儲存圖片
52 setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/3. Data/20250625 data_seurat/1. Bar plot")
53 ggsave(file = paste0(my_gene_list[i], "_boxplot.png"),
54 plot = bxp, width = 10, height = 5, dpi = 600)
55}
3-4. 主成分熱圖 (PC Heatmap)
主成分分析 (PCA) 不僅可以用於降維，還能告訴我們哪些基因對樣本間的主要差異貢獻最大。Seurat 的 DimHeatmap 函數可以將貢獻最大的基因（無論是正貢獻還是負貢獻）以熱圖的形式展示出來。
1r 3-4-pc-heatmap, eval=FALSE}
2# ... (省略載入套件) ...
3setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS/")
4da0 <- readRDS("2. GPL17586_12_sample_250625.rds", refhook = NULL)
5
6# --- 使用 Seurat 內建函數繪製 PC 熱圖 ---
7# https://satijalab.org/seurat/reference/dimheatmap
8Idents(da0) <- da0$`Sample Name`
9# 繪製對第一個主成分 (PC1) 貢獻最大的基因的熱圖
10# cells = 12 表示顯示所有12個樣本
11# balanced = TRUE 表示同時顯示正貢獻和負貢獻最大的基因
12PC_1 <- DimHeatmap(da0, dims = 1, cells = 12, balanced = TRUE, fast = F)
13PC_1
14
15# --- 提取數據並使用 pheatmap 重新繪製 ---
16# DimHeatmap 很方便，但客製化選項有限。我們可以提取其數據，用功能更強的 pheatmap 套件重畫。
17PC_1_data <- PC_1$data # 提取熱圖數據
18
19# ... (中間省略了使用 dplyr 和 tidyr 將數據從長格式轉換為寬格式的程式碼) ...
20# 最終我們會得到一個名為 df2 的矩陣，行是基因，欄是樣本
21
22# 使用 pheatmap 繪製熱圖
23library(pheatmap)
24pheatmap(df2)
25
26# --- 使用 tinyarray 繪製熱圖 ---
27# tinyarray 是另一個方便繪製熱圖的套件
28library(tinyarray)
29group_list <- colnames(df2)
30draw_heatmap(df2, group_list, cluster_cols = F, # 不對欄 (樣本) 進行聚類
31 legend = T, show_rownames = T)
結論
本教學文件完整地走過了一個典型的微陣列數據分析流程。從最原始的 .CEL 檔案開始，我們學習了如何使用 oligo 套件進行數據讀取和 RMA 標準化，如何處理多探針映射到單一基因的問題，以及如何將傳統的微陣列數據導入 Seurat 物件以利用其強大的分析和可視化工具。最後，我們展示了多種可視化方法，包括 UMAP、盒狀圖和熱圖，來幫助我們探索和解釋數據。
希望這份教學能對您的研究和學習有所幫助！
Comments