1# 清理工作環境中所有現存的物件，確保一個乾淨的開始
 2rm(list = ls(all = T))
 3
 4# 載入分析所需的核心 R 套件
 5library(dplyr)         # 提供強大的資料處理與操作功能 (如 filter, mutate, arrange)
 6library(Seurat)        # 單細胞 RNA 定序數據分析的核心工具
 7library(ggplot2)       # 繪製高品質圖形的基礎套件
 8library(patchwork)     # 組合多個 ggplot 圖形的工具
 9library(sctransform)   # Seurat 中用於正規化數據的另一種方法
10library(clusterProfiler) # 進行基因功能富集分析 (GO, KEGG) 的主要工具
11library(cellcall)      # (此腳本中載入但未使用，推測為其他分析流程所需)
12
13# 設定工作目錄，所有檔案的讀取與儲存都將以此路徑為基準
14# 請務必將 "C:/Users/TPOW31714/Desktop/..." 修改為您自己的檔案路徑
15setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS/")
16
17# 從 RDS 檔案中讀取 Seurat 物件
18# RDS 是 R 語言特有的檔案格式，可以儲存任何 R 物件
19da0 <- readRDS("2. GPL17586_12_sample_250625.rds", refhook = NULL)

 1# 1. 正規化 (Normalization)
 2#    目的：消除因測序深度不同造成的技術性差異
 3da0 <- NormalizeData(da0)
 4
 5# 2. 尋找高變異基因 (Find Variable Features)
 6#    目的：找出在不同細胞間表現量差異最大的基因，這些基因最可能用來區分細胞類型
 7da0 <- FindVariableFeatures(da0)
 8
 9# 3. 規模化 (Scaling)
10#    目的：將每個基因的表現量進行標準化 (平均值為0，變異數為1)，避免高表現量基因主導後續分析
11da0 <- ScaleData(da0, features = rownames(da0))

 1# 設定預設要分析的 Assay (通常是 'RNA')
 2DefaultAssay(da0) <- 'RNA'
 3
 4# 設定細胞的身分 (Identity)，這裡我們使用元數據 (metadata) 中的 "Group_1" 欄位
 5# "Group_1" 欄位應包含 "2D" 和 "3D" 的標籤
 6Idents(da0) <- factor(da0$Group_1)
 7
 8# 檢查各組的細胞數量
 9table(da0$Group_1) %>% as.data.frame()
10
11# 定義要比較的目標群組
12cell_groups <- c("3D")
13
14# 使用 for 迴圈來找出差異基因
15# 雖然這裡只比較一組，但迴圈的寫法提供了未來擴充的彈性
16for (cell_group in cell_groups) {
17  # 核心步驟：使用 FindMarkers 找差異基因
18  da0_markers_1 <- FindMarkers(da0, 
19                               ident.1 = cell_group, # 實驗組
20                               ident.2 = "2D",       # 對照組 (此處腳本原為NULL，改為"2D"更明確)
21                               only.pos = FALSE)     # only.pos=F 表示同時找出上調與下調基因
22
23  # 設定輸出目錄
24  setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/3. Data/20250625 data_ORA/1. 3Dvs2D DEGs_12samples")
25  
26  # 將結果寫入文字檔
27  write.table(da0_markers_1, 
28              paste0("1. Findmarkers_", cell_group , " vs 2D_20250625.txt"), 
29              sep = "\t",           # 使用 tab 作為分隔符
30              row.names = TRUE,    # 保留基因名稱作為列名
31              col.names = TRUE,    # 保留欄位名稱
32              quote = FALSE)       # 不對字串加上引號
33}
34
35# 清理不再需要的變數
36rm(da0_markers_1)

 1# 清理環境，重新載入必要的套件
 2rm(list = ls(all = T))
 3library(dplyr)
 4library(Seurat)
 5library(ggplot2)
 6library(patchwork)
 7library(clusterProfiler)
 8library(enrichplot)
 9library(org.Hs.eg.db) # 人類基因註解資料庫
10
11# 重新讀取與前處理資料 (此為重複步驟，在實際流程中可省略)
12setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/4. RDS/")
13da0 <- readRDS("2. GPL17586_12_sample_250625.rds", refhook = NULL)
14da0 <- NormalizeData(da0)
15da0 <- FindVariableFeatures(da0)
16da0 <- ScaleData(da0, features = rownames(da0))
17DefaultAssay(da0) <- 'RNA'
18Idents(da0) <- factor(da0$Group_1)
19
20# 核心步驟：使用 FindAllMarkers 找出各群組的標記基因
21# only.pos = TRUE 表示我們只關心在該群組中表現量顯著較高的基因
22da0_markers <- FindAllMarkers(da0, only.pos = TRUE)

 1# 使用 dplyr 篩選每個 cluster 中 avg_log2FC 最高的 500 個基因
 2TOP_DEGs <- da0_markers %>%
 3  group_by(cluster) %>% # 根據 cluster (即 2D, 3D) 進行分組
 4  slice_max(order_by = avg_log2FC, n = 500, with_ties = FALSE) %>% # 選出 avg_log2FC 最大的前 500 個
 5  ungroup() # 取消分組
 6
 7# 使用 DoHeatmap 繪製熱圖
 8p1 <- DoHeatmap(da0, features = TOP_DEGs$gene) + NoLegend() # NoLegend() 移除圖例
 9p1 # 顯示圖形
10
11# 設定存圖目錄並儲存
12setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/3. Data/20250625 data_ORA/1. 3Dvs2D DEGs_12samples/RDS/")
13ggsave(file = "HEATMAP_12samples_DEG_TOP500_250625v1.png",
14       plot = p1, width = 9, height = 15, dpi = 300, limitsize = FALSE)
15
16# 清理變數
17rm(da0_markers, TOP_DEGs, p1)

 1# 建立輸出目錄
 2dir.create("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/3. Data/20250625 data_ORA/1. 3Dvs2D DEGs_12samples/RDS/", recursive = TRUE)
 3
 4# 重新執行 FindAllMarkers，這次包含下調基因 (only.pos = F)
 5da0_markers <- FindAllMarkers(da0, only.pos = FALSE)
 6
 7# 定義要處理的細胞群組
 8cell_groups <- c("2D", "3D")
 9# 初始化一個空的 list 來儲存結果
10results_list <- list()
11
12# 透過迴圈處理每個細胞群組
13for (cell_group in cell_groups) {
14  
15  # 1. 從總表中篩選出當前群組的 markers
16  markers <- da0_markers %>% filter(cluster == cell_group)
17  
18  # 2. 對 markers 進行處理
19  filtered_markers <- markers %>%
20    mutate(change = ifelse(avg_log2FC > 0, "UP", "DOWN")) %>% # 增加 UP/DOWN 標籤
21    arrange(desc(avg_log2FC)) # 根據 logFC 降序排列
22  
23  # 3. 將基因 SYMBOL 轉換為 ENTREZID
24  de <- filtered_markers$gene
25  ids <- bitr(de, 
26              fromType = "SYMBOL",   # 來源ID類型
27              toType = c("ENTREZID"),# 目標ID類型
28              OrgDb = "org.Hs.eg.db") # 指定物種為人類
29  
30  # 4. 將轉換後的 ID 加回資料框，並過濾掉無法轉換的基因
31  filtered_markers <- filtered_markers %>%
32    left_join(ids, by = c("gene" = "SYMBOL")) %>% # 將 ids 合併進來
33    filter(!is.na(ENTREZID)) # 移除沒有對應 ENTREZID 的基因
34  
35  # 5. 建立最終的基因列表 (只取上調基因)
36  gene_list <- filtered_markers %>%
37    filter(avg_log2FC > 0) %>%
38    arrange(desc(avg_log2FC)) %>%
39    head(500) %>% # 使用前500上調基因
40    dplyr::select(ENTREZID, avg_log2FC)
41  
42  # 6. 建立 ORA 所需的 "named vector"
43  gene_list_values <- gene_list$avg_log2FC
44  names(gene_list_values) <- gene_list$ENTREZID
45  
46  # 7. 將結果存入 results_list
47  results_list[[cell_group]] <- gene_list_values
48}
49
50# 將準備好的基因列表儲存為 RDS 檔案
51setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/3. Data/20250625 data_ORA/1. 3Dvs2D DEGs_12samples/RDS/")
52saveRDS(results_list, file = "1. Combined_list_12samples_DEGs_POS_TOP500_250625.rds")

 1# 清理環境並載入套件
 2rm(list = ls(all = T))
 3library(dplyr)
 4library(ggplot2)
 5library(clusterProfiler)
 6library(enrichplot)
 7library(org.Hs.eg.db)
 8
 9# 設定工作目錄並讀取基因列表
10setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/3. Data/20250625 data_ORA/1. 3Dvs2D DEGs_12samples/RDS/")
11markers_list <- readRDS("1. Combined_list_12samples_DEGs_POS_TOP500_250625.rds", refhook = NULL)
12
13# 將 list 中的 named vector 轉換成只包含基因 ID 的 character vector
14markers_list_converted <- lapply(markers_list, function(x) names(x))
15str(markers_list_converted) # 檢查轉換後的格式
16
17# 核心步驟：執行 GO 富集分析
18ora.out <- compareCluster(geneClusters = markers_list_converted,
19                          fun = "enrichGO",         # 指定使用 GO 分析
20                          OrgDb = org.Hs.eg.db,     # 指定物種為人類
21                          ont = "BP",               # 指定 GO 的本體 (BP, MF, CC, ALL)
22                          pvalueCutoff = 0.05,      # p-value 閾值
23                          pAdjustMethod = "none")   # p-value 校正方法
24
25# 儲存分析結果
26setwd("C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/3. Data/20250625 data_ORA/1. 3Dvs2D DEGs_12samples/RDS/")
27saveRDS(ora.out, file = "1-2. GOBP_compare_12samples_DEGs_POS_TOP500_250625.rds")

 1# ... (此處省略重複的 rm 和 library 載入) ...
 2# 讀取 ORA 結果
 3xx <- readRDS("1-2. GOBP_compare_12samples_DEGs_POS_TOP500_250625.rds", refhook = NULL)
 4
 5# 1. 繪製點圖 (Dot Plot)
 6p1 <- dotplot(xx, showCategory = 10, label_format = 80)
 7p1 <- p1 + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
 8p1
 9# 存檔
10ggsave("Dotplot_GOBP_TOP10_ORA_DEGs_POS_TOP500_250625.png", p1, width = 12, height = 16, dpi = 600)
11
12# 2. 將基因 ID 轉回可讀的 Symbol 並匯出文字檔
13xx2 <- setReadable(xx, 'org.Hs.eg.db', 'ENTREZID')
14write.table(xx2@compareClusterResult, "1-2. GOBP_compare_12samples_DEGs_POS_TOP500_250625.txt", sep = "\t", row.names = T, col.names = T, quote = F)
15
16# 3. 繪製 Emapplot (Enrichment Map)
17#    Emapplot 可以將功能相似的 GO terms 聚合在一起，形成網絡圖
18#    首先需要計算 GO terms 之間的相似性
19xx2_sim <- enrichplot::pairwise_termsim(xx2)
20
21# 篩選出特定群組 (例如 "2D") 的結果來繪圖
22Three_Dim_result <- xx2_sim
23Three_Dim_result@compareClusterResult <- Three_Dim_result@compareClusterResult %>%
24                                         filter(Cluster == "2D")
25
26# 當結果只剩一個群組時，可以像處理 enrichResult 物件一樣繪圖
27p_emap <- emapplot(Three_Dim_result, layout = "kk", cex_category = 1.5)
28p_emap

 1# ... (省略 rm 和 library) ...
 2library(readxl)
 3library(openxlsx)
 4library(scales)
 5
 6cell_groups <- c("2D", "3D")
 7plots_list <- list()
 8data_list <- list()
 9
10# 讀取 enrichGO 分析結果 (這裡假設已將 txt 轉存為 Excel)
11# 注意：原腳本此處讀取 Excel，若無此檔案會報錯。
12# 我們直接使用上一步的 R 物件 `xx`
13go_analysis_1 <- xx@compareClusterResult
14
15for (cell_group in cell_groups) {
16  
17  # 篩選出當前群組的結果
18  go_analysis_2 <- go_analysis_1 %>%
19    filter(Cluster == cell_group) %>%
20    arrange(p.adjust) %>%
21    mutate(Description = factor(Description, levels = rev(unique(Description))))
22  
23  # 繪製點圖
24  p1 <- ggplot(go_analysis_2, aes(x = -log10(p.adjust), y = Description)) +
25    geom_point(aes(size = Count, color = p.adjust)) +
26    theme_bw() +
27    scale_color_gradient(low = "red", high = "blue") +
28    labs(title = paste("GO BP Enrichment for", cell_group),
29         x = "-log10(p.adjust)", y = "")
30  
31  plots_list[[cell_group]] <- p1
32  data_list[[cell_group]] <- go_analysis_2
33  
34  # 顯示圖形
35  print(p1)
36}
37
38# 將各群組結果寫入同一個 Excel 檔案的不同工作表
39output_file <- "C:/Users/TPOW31714/Desktop/20250224 Human microarray for CY_add cell/3. Data/20250625 data_ORA/1. 3Dvs2D DEGs_12samples/RDS/1. Origianl GO_DEGs_POS/POS_results.xlsx"
40wb <- createWorkbook()
41for (cell_group in names(data_list)) {
42  addWorksheet(wb, sheetName = cell_group)
43  writeData(wb, sheet = cell_group, x = data_list[[cell_group]])
44}
45saveWorkbook(wb, file = output_file, overwrite = TRUE)

 1# ... (省略 rm 和 library) ...
 2# 讀取 ORA 結果
 3xx <- readRDS("1-2. GOBP_compare_12samples_DEGs_POS_TOP500_250625.rds", refhook = NULL)
 4
 5# 1. 取得資料框
 6df <- xx@compareClusterResult
 7
 8# 2. 找出共同的 pathway ID
 9all_clusters <- unique(df$Cluster)
10num_clusters <- length(all_clusters)
11common_ids <- df %>%
12  group_by(ID) %>%
13  summarise(n = n_distinct(Cluster)) %>%
14  filter(n == num_clusters) %>%
15  pull(ID)
16
17# 3. 分割資料
18df_common <- df %>% filter(ID %in% common_ids)
19df_noncommon <- df %>% filter(!ID %in% common_ids)
20
21# 4. 對非共同部分進行 simplify
22xx_noncommon <- xx
23xx_noncommon@compareClusterResult <- df_noncommon
24xx_noncommon_simplified <- simplify(xx_noncommon,
25                                    cutoff = 0.7, # 相似度閾值
26                                    by = "p.adjust", # 根據 p.adjust 選擇代表性 term
27                                    select_fun = min,
28                                    measure = "Wang") # 計算語義相似度的方法
29
30# 5. 合併結果
31df_combined <- rbind(xx_noncommon_simplified@compareClusterResult, df_common)
32xx_combined <- xx
33xx_combined@compareClusterResult <- df_combined
34
35# 儲存簡化後的結果
36saveRDS(xx_combined, file = "1-3. GOBP_SIMPLIFY_12samples_DEGs_POS_TOP500_250625.rds")

1# 讀取簡化後的結果
2xx_simplified <- readRDS("1-3. GOBP_SIMPLIFY_12samples_DEGs_POS_TOP500_250625.rds", refhook = NULL)
3
4# 繪製點圖
5p_simplified <- dotplot(xx_simplified, showCategory = 10, label_format = 80)
6p_simplified
7
8# 同樣可以分組繪製 Emapplot 或匯出成 Excel
9# ... (程式碼與步驟 4.2, 4.3 類似，此處省略) ...