본문 바로가기
AI 및 Data Analysis/Code

GSE62452 (microarray) Analysis Results

by doraemin_dev 2025. 5. 7.

GSE62452 (microarray) Dataset

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62452

 

GEO Accession viewer

 

www.ncbi.nlm.nih.gov

 

Analysis Code

https://github.com/doraemon49/lab-data-ai/tree/main/GSE_Analysis

 

lab-data-ai/GSE_Analysis at main · doraemon49/lab-data-ai

연구실 기록. Contribute to doraemon49/lab-data-ai development by creating an account on GitHub.

github.com

 

Analysis Results

 

Analysis Code

test <- c(-1, 0, 1)
test

# install.packages("installr")
# library(installr)
# check.for.updates.R()
# install.R()
# .libPaths()

# 파일 경로 설정
GSE62452 <- "/home/kim89/GES62452.csv"

# CSV 파일 읽기
data <- read.csv(GSE62452, header = TRUE)

# Data 정상화 (Optional)
# RNA-seq일 경우, 데이터가 로그 변환되어 있다면 이 과정은 생략 가능
# For array data or non-log-transformed data:
# data <- log2(data + 1)

# 데이터 분포 확인
hist(data$X.E005.Tc8., main="Distribution of X.E005.Tc8.", xlab="Expression levels")
summary(data$X.E005.Tc8.)

# 데이터 확인
print(head(data))

#     GeneSymbol  ID_REF X.E005.Tc8. X.E005.Tn. X.E012.Tp2. X.E12.Tn. X.E021.Tc8.
# 1       OR4F17 7896740     1.98519    2.22654     2.20386   1.89503     2.22880
# 2 LOC100134822 7896742     4.27115    5.08438     4.40239   5.53027     5.13266
# 3       OR4F29 7896744     4.10218    4.26139     4.20362   4.47237     3.84170
# 4 LOC100287934 7896754     5.37197    5.68217     6.17415   4.83325     4.60658
# 5       FAM87B 7896756     2.66610    3.39605     2.86317   2.79406     2.82701
# 6    LINC01128 7896759     4.59326    4.53481     4.31149   4.60532     4.06294

####  limma   ######################################################
# 필요한 패키지 설치 및 로드
if (!requireNamespace("limma", quietly = TRUE)) {
    install.packages("limma")
}
library(limma)

# 데이터 불러오기
data <- read.csv(GSE62452, header = TRUE)
# GeneSymbol 및 ID_REF 열 제외한 발현량 데이터만 선택
# 세 번째 열부터가 발현 데이터이므로 [, -c(1, 2)]로 첫 번째, 두 번째 열을 제외
data_expression <- data[, -c(1, 2)]

# GeneSymbol에 중복이 있는지 확인 
sum(duplicated(data$GeneSymbol)) # 2193개 중복...
gene_symbols <- data$GeneSymbol  # GeneSymbol 열만 따로 저장
# GeneSymbol을 기준으로 발현량 평균 계산
data_aggregated <- aggregate(data_expression, by = list(GeneSymbol = gene_symbols), FUN = mean)
sum(duplicated(data_aggregated$GeneSymbol)) # 0 # 이제 중복된 GeneSymbol 없음

# 결과에서 GeneSymbol 열을 추출해 rownames로 설정
rownames(data_aggregated) <- data_aggregated$GeneSymbol

# 더 이상 필요 없는 GeneSymbol 열을 제거
data_numeric <- data_aggregated[, -1]
nrow(data_numeric) # 20002 # 22195개에서 줄어듦

# head(rownames(data_numeric))
# nrow(data_numeric) # 22195 rows
# head(colnames(data_numeric))
# ncol(data_numeric) # 130 columns

# 샘플 조건 정의
# group 정보는 각 샘플이 속한 조건을 지정하는 벡터로 정의
# 샘플의 순서대로, Tumor인지 Nomal인지 기입해준다.
group <- factor(c(
  "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N",
  "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N",
  "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N",
  "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N",
  "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "T", "N",
  "T", "N", "T", "N", "T", "N", "T", "N", "T", "N", "N", "T", "N", "T", "T", "T",
  "N", "T", "N", "T", "T", "N", "T", "N", "T", "T", "N", "T", "T", "N", "T", "T",
  "N", "T", "T", "T", "N", "N", "T", "N", "T", "N", "T", "T", "N", "T", "N", "T",
  "N", "T"
))
# Design Matrix 생성
design <- model.matrix(~ 0 + group)
# head(design)
#   groupN groupT
# 1      0      1
# 2      1      0


##############
# 이제, Linear Model fitting
fit <- lmFit(data_numeric, design)

# Contrasts 설정 
contrast.matrix <- makeContrasts(Tumor_vs_Nomal = groupT - groupN, levels = design)

# Contrasts 적용
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)

# limma DEGs 결과 확인
deg_limma <- topTable(fit2, adjust.method = "BH", number = Inf)
head(deg_limma, 10) # 상위 10개의 DEGs 출력
#            logFC  AveExpr        t      P.Value    adj.P.Val        B
# LAMC2   2.646110 6.061879 11.32408 3.550988e-21 7.102686e-17 37.45794
# LAMB3   2.084578 5.486133 10.92686 3.525125e-20 3.525478e-16 35.22765
# TSPAN1  2.384041 5.796048 10.72418 1.136131e-19 7.574961e-16 34.09015

# adf.P.Val < 0.05 및  abs(log2 FC) > 1.5인 DEG 필터링
deg_limma_filtered <- deg_limma[deg_limma$adj.P.Val < 0.05 & abs(deg_limma$logFC) > 1.5, ]
nrow(deg_limma_filtered)  # 필터링된 DEG 개수 확인 # 88 # 20002개 중에서 선별됨.
rownames(deg_limma_filtered)
head(deg_limma_filtered)   # 필터링된 상위 DEGs 확인
#           logFC  AveExpr        t      P.Value    adj.P.Val        B
# LAMC2  2.646110 6.061879 11.32408 3.550988e-21 7.102686e-17 37.45794
# LAMB3  2.084578 5.486133 10.92686 3.525125e-20 3.525478e-16 35.22765
# TSPAN1 2.384041 5.796048 10.72418 1.136131e-19 7.574961e-16 34.09015


# 결과 파일로 저장
# write.csv(deg_limma_filtered, file = "deg_limma_results.csv")


####  limma - Paired   ######################################################
# 파일 경로 설정
GSE62452 <- "/home/kim89/GES62452.csv"

# CSV 파일 읽기
data <- read.csv(GSE62452, header = TRUE)

# 데이터 분포 확인
hist(data$X.E005.Tc8., main="Distribution of X.E005.Tc8.", xlab="Expression levels")
summary(data$X.E005.Tc8.)

# 데이터 확인
# print(head(data))

#     GeneSymbol  ID_REF X.E005.Tc8. X.E005.Tn. X.E012.Tp2. X.E12.Tn. X.E021.Tc8.
# 1       OR4F17 7896740     1.98519    2.22654     2.20386   1.89503     2.22880
# 2 LOC100134822 7896742     4.27115    5.08438     4.40239   5.53027     5.13266
# 3       OR4F29 7896744     4.10218    4.26139     4.20362   4.47237     3.84170

# 필요한 패키지 설치 및 로드
if (!requireNamespace("limma", quietly = TRUE)) {
    install.packages("limma")
}
library(limma)

# GeneSymbol 및 ID_REF 열 제외한 발현량 데이터만 선택
data_expression <- data[, -c(1, 2)]

# GeneSymbol에 중복이 있는지 확인 # 2193개 중복...
sum(duplicated(data$GeneSymbol))
gene_symbols <- data$GeneSymbol      # GeneSymbol 열만 따로 저장
# GeneSymbol을 기준으로 발현량 평균 계산
data_aggregated <- aggregate(data_expression, by = list(GeneSymbol = gene_symbols), FUN = mean)
sum(duplicated(data_aggregated$GeneSymbol)) # 0 # 이제 중복된 GeneSymbol 없음

# 결과에서 GeneSymbol 열을 추출해 rownames로 설정
rownames(data_aggregated) <- data_aggregated$GeneSymbol

# 더 이상 필요 없는 GeneSymbol 열을 제거
data_numeric <- data_aggregated[, -1]
nrow(data_numeric) # 20002 # 22195에서 줄어듦.

# Pair 정보를 읽기
# 필요한 패키지 설치 및 로드
if (!require("readxl")) install.packages("readxl")
library(readxl)
pair_data <- read_excel("/home/kim89/GSE62452_paired_sample_information.xlsx")
pair_info <- pair_data[,c("Sample", "Paired Sample ID")]
# nrow(pair_info) # 130 rows
# 샘플 이름을 맞추기 위해 열 이름을 수정
colnames(pair_info) <- c("Sample", "Pair") 

# 데이터를 pair 정보와 매칭
paired_samples <- pair_info[complete.cases(pair_info), ]  # 쌍이 있는 데이터만 선택
# nrow(paired_samples) # 120 rows
# colnames(data_numeric)에서 "X." 제거 및 마침표를 하이픈으로 변환
adjusted_colnames <- gsub("^X\\.", "", colnames(data_numeric))  # "X." 제거
adjusted_colnames <- gsub("\\.(?=\\w)", "-", adjusted_colnames, perl = TRUE)  # 중간 마침표만 하이픈으로 변환 (lookahead 사용)
adjusted_colnames <- gsub("\\.$", "", adjusted_colnames)  # 마지막 마침표 제거
# 샘플 이름을 맞추기 위해 수정된 열 이름으로 비교
matched_samples <- adjusted_colnames %in% paired_samples$Sample
# print(sum(matched_samples)) # 120개. 

# matched_samples의 TRUE인 열만 선택되어 paired_data에 저장
paired_data <- data_numeric[, matched_samples]
# length(paired_data) # pair이 있는 120개
# length(data_numeric) # 전체 130개

# 그룹과 쌍 정보 생성 # 샘플 이름에 "Tc","Tp" 또는 마지막에 "T"가 가 포함된 경우 그 샘플은 "Tumor"로, 그렇지 않은 경우는 "Normal"로 분류
group <- factor(ifelse(grepl("Tc|Tp|T$", paired_samples$Sample), "Tumor", "Normal"))
# group 데이터뿐만아니라, pair 데이터까지 고려함.
pair <- factor(paired_samples$Pair)
# head(group)
# [1] Tumor  Normal Tumor  Normal Tumor  Normal
# Levels: Normal Tumor
# head(pair)
# [1] 30 30 35 35 46 46
# 60 Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 ... 60

# 디자인 매트릭스 생성 (paired analysis) # 절편은 제외 ~0
design <- model.matrix(~ 0 + group + pair)
# head(design) # Tumor이면 'groupTumor' 값이 1. # pair n 의 값이 1
  # groupNormal groupTumor pair2 pair3 pair4 pair5 pair6 pair7 pair8 pair9 pair10 ... pair30 pair31 pair32 pair33 pair34 pair35
# 1          0          1     0     0     0     0     0     0     0     0      0 ...      1      0      0      0      0      0
# 2          1          0     0     0     0     0     0     0     0     0      0 ...      1      0      0      0      0      0
# 3          0          1     0     0     0     0     0     0     0     0      0 ...      0      0      0      0      0      1
# 4          1          0     0     0     0     0     0     0     0     0      0 ...      0      0      0      0      0      1

# Linear Model fitting
fit <- lmFit(paired_data, design)

# Contrasts 설정 (Tumor vs Normal 비교)
contrast.matrix <- makeContrasts(Tumor_vs_Normal = groupTumor - groupNormal, levels = design)

# Contrasts 적용
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)

# DEGs 결과 확인
deg_limma_pair <- topTable(fit2, adjust.method = "BH", number = Inf)
head(deg_limma_pair, 10)  # 상위 10개의 DEGs 출력
nrow(deg_limma_pair)  # 20002

# adj.P.Val < 0.05 및 abs(logFC) > 1.5인 DEG 필터링
deg_limma_pair_filtered <- deg_limma_pair[deg_limma_pair$adj.P.Val < 0.05 & abs(deg_limma_pair$logFC) > 1.5, ]

nrow(deg_limma_pair_filtered)  # 필터링된 DEG 개수 확인 # 88 # 20002개 중에서 선별됨.
rownames(deg_limma_pair_filtered)
head(deg_limma_pair_filtered)  # 필터링된 상위 DEGs 확인
#            logFC  AveExpr        t      P.Value    adj.P.Val        B
# TSPAN1  2.511219 5.789988 13.48984 3.157358e-20 4.406567e-16 35.51121
# TMPRSS4 2.165359 4.779879 13.42411 3.970774e-20 4.406567e-16 35.29028
# LAMC2   2.757130 6.049321 13.09952 1.241606e-19 9.185815e-16 34.19055

# 결과 파일로 저장
# write.csv(deg_limma_pair, file = "deg_limma_paired_results_all.csv")
# write.csv(deg_limma_pair_filtered, file = "deg_limma_paired_results.csv")


############ 방법 1과 2의 공통 결과 ##########
# 공통된 유전자 이름 찾기
common_genes <- intersect(rownames(deg_limma_filtered), rownames(deg_limma_pair_filtered))
# 공통된 유전자의 데이터프레임 생성
common_data_limma <- deg_limma_filtered[common_genes, ]
# 공통된 유전자 이름 출력
print(common_genes)
length(common_genes) # 83개

# 결과 파일로 저장
# write.csv(common_genes, file = "deg_limma_common_results.csv")
# write.csv(common_data_limma, file = "deg_limma_common.csv")

# ggplot2 패키지 로드
library(ggplot2)

# common gene의 volcano plot
ggplot(common_data_limma, aes(x = logFC, y = -log10(adj.P.Val))) +
  geom_point(aes(color = logFC > 0), size = 4) +
  scale_color_manual(values = c("blue", "red")) + 
  theme_minimal() +
  labs(title = "Volcano Plot of Common DEGs",
       x = "Log2 (Fold Change)",
       y = "-log10 (Adjusted P-Value)") +
  geom_hline(yintercept = -log10(0.05), color = "red", linetype = "dashed") + # p-value 기준선 추가
  theme(plot.title = element_text(hjust = 0.5))

  # 상위 10개 유전자를 adj.P.Val 값으로 정렬 후 선택
top_genes <- head(common_data_limma[order(common_data_limma$adj.P.Val), ], 30)

# Volcano plot에 상위 10개 유전자 이름을 라벨링하는 코드
ggplot(common_data_limma, aes(x = logFC, y = -log10(adj.P.Val))) +
  geom_point(aes(color = logFC > 0), size = 4) +
  scale_color_manual(values = c("blue", "red")) + 
  theme_minimal() +
  labs(title = "Volcano Plot of Common DEGs",
       x = "Log2 (Fold Change)",
       y = "-log10 (Adjusted P-Value)") +
  geom_hline(yintercept = -log10(0.05), color = "red", linetype = "dashed") + # p-value 기준선 추가
  theme(plot.title = element_text(hjust = 0.5)) +
    theme(
    plot.title = element_text(hjust = 0.5, size = 35),  # 제목 글씨 크기
    axis.title.x = element_text(size = 25),             # x축 제목 글씨 크기
    axis.title.y = element_text(size = 25),             # y축 제목 글씨 크기
    axis.text.x = element_text(size = 25),              # x축 숫자 크기
    axis.text.y = element_text(size = 25)               # y축 숫자 크기
  ) +
  # 상위 10개 유전자에 라벨 추가
  geom_text(data = top_genes, aes(label = rownames(top_genes)), 
            vjust = 1, hjust = 1.2, size = 8, color = "black")

# # 히트맵 필요한 패키지 설치
# if (!requireNamespace("pheatmap", quietly = TRUE)) {
#   install.packages("pheatmap")
# }
# library(pheatmap)

# # 히트맵에 사용할 데이터를 준비합니다.
# # common gene들의 발현 데이터를 준비합니다.
# # 데이터는 common_data_limma에서 logFC 값들로 사용한다고 가정합니다.
# heatmap_data <- common_data_limma[, "logFC", drop = FALSE]  # 열로 구성된 logFC 데이터 선택
# rownames(heatmap_data) <- rownames(common_data_limma)       # rownames 설정
# # 스케일링을 제거하고 히트맵을 다시 그립니다
# pheatmap(heatmap_data, 
#          cluster_rows = TRUE,    # 행 기준 클러스터링 여부
#          cluster_cols = FALSE,   # 열 기준 클러스터링 여부 (열이 하나이므로 비활성화)
#          scale = "none",         # 스케일링 제거
#          show_rownames = TRUE,   # 유전자 이름 표시
#          color = colorRampPalette(c("blue", "white", "red"))(100),  # 색상 팔레트 설정
#          main = "Heatmap of Common Genes (logFC)")  # 히트맵 제목
# # ggplot2를 사용한 막대 그래프
# library(ggplot2)

# # logFC 막대 그래프
# ggplot(heatmap_data, aes(x = rownames(heatmap_data), y = logFC)) +
#   geom_bar(stat = "identity", fill = "steelblue") +
#   theme_minimal() +
#   labs(title = "Bar Plot of Common Genes (logFC)", x = "Genes", y = "Log2 Fold Change") +
#   theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))  # 유전자 이름 회전

# 예시 데이터: DAVID에서 추출한 GO term 데이터 프레임 (예: df)
# df는 GO Term, PValue, Count 등의 열을 포함한다고 가정

go_bp <- read.table("/home/kim89/DAVID_BP.txt", header = TRUE, sep = "\t")
go_cc <- read.table("/home/kim89/DAVID_CC.txt", header = TRUE, sep = "\t")
go_mf <- read.table("/home/kim89/DAVID_MF.txt", header = TRUE, sep = "\t")
go_all <- rbind(go_bp, go_cc, go_mf)
# Bubble Plot (버블 플롯):
# 점 플롯과 유사하지만, 점의 크기와 색상을 동시에 활용하여 유전자 수와 P-value를 함께 시각화합니다.
ggplot(go_all, aes(x = Category, y = reorder(Term, -Count), size = Count, color = PValue)) + 
    geom_point(alpha = 0.7) + 
    theme_minimal() +
    labs(title = "Bubble Plot for GO Term Results",
         x = "Category",
         y = "GO Term") +
    scale_color_gradient(low = "blue", high = "red") +
    scale_size(range = c(3, 10)) +  # 버블 크기 조절
    theme(
    plot.title = element_text(hjust = 0.5, size = 20),  # 제목 글씨 크기
    axis.title.x = element_text(size = 15),             # x축 제목 글씨 크기
    axis.title.y = element_text(size = 15),             # y축 제목 글씨 크기
    axis.text.x = element_text(size = 15),              # x축 숫자 크기
    axis.text.y = element_text(size = 15)               # y축 숫자 크기
  )

# 범주별로 구분하는 열 추가
go_bp$Category <- "BP"
go_cc$Category <- "CC"
go_mf$Category <- "MF"

# P-value와 FDR 기준으로 필터링
filtered_go_bp <- go_bp[go_bp$PValue < 0.05 & go_bp$FDR < 0.05, ]
filtered_go_cc <- go_cc[go_cc$PValue < 0.05 & go_cc$FDR < 0.05, ]
filtered_go_mf <- go_mf[go_mf$PValue < 0.05 & go_mf$FDR < 0.05, ]
filtered_go_all <- rbind(filtered_go_bp, filtered_go_cc, filtered_go_mf)

# 필요한 라이브러리 설치 및 불러오기
# install.packages("ggplot2")
library(ggplot2)

# 데이터 병합

# 시각화
ggplot(filtered_go_all, aes(x = reorder(Term, -Count), y = Count, fill = PValue)) + 
    geom_bar(stat = "identity") + 
    coord_flip() +  # 수평 바 차트
    theme_minimal() +
    labs(title = "DAVID GO Term Results (BP, CC, MF)",
         x = "GO Term",
         y = "Gene Count") +
    facet_wrap(~ Category, scales = "free") +  # 범주별로 나눠서 표시
    scale_fill_gradient(low = "blue", high = "red")

# Dot Plot (점 그래프):
ggplot(filtered_go_all, aes(x = reorder(Term, -Count), y = Count, color = PValue, size = Count)) + 
    geom_point() + 
    coord_flip() +  # 수평으로 정렬
    theme_minimal() +
    labs(title = "Dot Plot for GO Term Results",
         x = "GO Term",
         y = "Gene Count") +
    scale_color_gradient(low = "blue", high = "red")

# Bubble Plot (버블 플롯):
# 점 플롯과 유사하지만, 점의 크기와 색상을 동시에 활용하여 유전자 수와 P-value를 함께 시각화합니다.
ggplot(filtered_go_all, aes(x = Category, y = reorder(Term, -Count), size = Count, color = PValue)) + 
    geom_point(alpha = 0.7) + 
    theme_minimal() +
    labs(title = "Bubble Plot for GO Term Results",
         x = "Category",
         y = "GO Term") +
    scale_color_gradient(low = "blue", high = "red") +
    scale_size(range = c(3, 10)) +  # 버블 크기 조절
    theme(
    plot.title = element_text(hjust = 0.5, size = 20),  # 제목 글씨 크기
    axis.title.x = element_text(size = 15),             # x축 제목 글씨 크기
    axis.title.y = element_text(size = 15),             # y축 제목 글씨 크기
    axis.text.x = element_text(size = 15),              # x축 숫자 크기
    axis.text.y = element_text(size = 15)               # y축 숫자 크기
  )

library(pheatmap)
# GO term과 유전자 수를 기반으로 히트맵 생성
gene_count_matrix <- as.matrix(filtered_go_all$Count)
rownames(gene_count_matrix) <- filtered_go_all$Term
pheatmap(gene_count_matrix, cluster_rows = TRUE, cluster_cols = FALSE)

library(clusterProfiler)
cnetplot(your_go_result, categorySize = "gene_count", foldChange = NULL)

# KEGG
kegg_data <- read.table("/home/kim89/KEGG_PATHWAY.txt", header = TRUE, sep = "\t")
library(ggplot2)

# 점 그래프 (Dot plot)
ggplot(kegg_data, aes(x = reorder(Term, -Count), y = Count, color = PValue, size = Count)) + 
  geom_point() + 
  coord_flip() +  # 수평 정렬
  theme_minimal() +
  labs(title = "KEGG Pathway Enrichment Results",
       x = "KEGG Pathway",
       y = "Gene Count") +
  scale_color_gradient(low = "blue", high = "red") +  # 색상 그라데이션
  scale_size(range = c(3, 10)) +  # 점 크기 조절
   theme(
    plot.title = element_text(hjust = 0.5, size = 20),  # 제목 글씨 크기
    axis.title.x = element_text(size = 15),             # x축 제목 글씨 크기
    axis.title.y = element_text(size = 15),             # y축 제목 글씨 크기
    axis.text.x = element_text(size = 15),              # x축 숫자 크기
    axis.text.y = element_text(size = 15)               # y축 숫자 크기
  )

 

GO, KEGG 

# DESeq2 패키지 설치 (설치되어 있지 않은 경우)
if (!requireNamespace("DESeq2", quietly = TRUE)) {
  install.packages("BiocManager")
  BiocManager::install("DESeq2")
}



# Install necessary packages if not already installed
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

#BiocManager::install('clusterProfiler') # 설치완료
#BiocManager::install("org.Hs.eg.db") # 설치 완료료

# Load the required libraries
library(clusterProfiler)
library(org.Hs.eg.db)

common_data_limma <- read.csv("C:\\Users\\LG\\Documents\\MJU\\연구실\\24.09.27_GSE62452\\deg_limma_common.csv", row.names = 1, header = TRUE)

# Extract gene symbols from your dataset `common_data_limma`
gene_symbols <- rownames(common_data_limma)

# Convert gene symbols to ENTREZ IDs using org.Hs.eg.db
eg <- bitr(gene_symbols, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)

# Save the conversion result to a CSV file
write.csv(eg, file = "C:\\Users\\LG\\Documents\\MJU\\연구실\\24.09.27_GSE62452\\common_deg_limma_entrez.csv")


# Display the result
head(eg)

# DAVID로 가서 결과 받아오자.

# DAVID에서 결과 받아옴
go_bp <- read.table("C:\\Users\\LG\\Documents\\MJU\\연구실\\24.09.27_GSE62452\\GO, KEGG\\DAVID_BP.txt", header = TRUE, sep = "\t")
go_cc <- read.table("C:\\Users\\LG\\Documents\\MJU\\연구실\\24.09.27_GSE62452\\GO, KEGG\\DAVID_CC.txt", header = TRUE, sep = "\t")
go_mf <- read.table("C:\\Users\\LG\\Documents\\MJU\\연구실\\24.09.27_GSE62452\\GO, KEGG\\DAVID_MF.txt", header = TRUE, sep = "\t")


# 범주별로 구분하는 열 추가
go_bp$Category <- "BP"
go_cc$Category <- "CC"
go_mf$Category <- "MF"

# P-value와 FDR 기준으로 필터링
filtered_go_bp <- go_bp[go_bp$PValue < 0.05 & go_bp$FDR < 0.05, ]
filtered_go_cc <- go_cc[go_cc$PValue < 0.05 & go_cc$FDR < 0.05, ]
filtered_go_mf <- go_mf[go_mf$PValue < 0.05 & go_mf$FDR < 0.05, ]
filtered_go_all <- rbind(filtered_go_bp, filtered_go_cc, filtered_go_mf)

library(clusterProfiler)

# GO 분석 수행 (BP는 Biological Process)
go_enrich <- enrichGO(gene = eg$ENTREZID, 
                      OrgDb = org.Hs.eg.db, 
                      keyType = "ENTREZID", 
                      ont = "BP", 
                      pvalueCutoff = 0.05)

# 네트워크 플롯 그리기
cnetplot(go_enrich, categorySize = "gene_count", color.params = list(foldChange = NULL))

'AI 및 Data Analysis > Code' 카테고리의 다른 글

Data Size Optimization  (0) 2025.07.07
[ScRAT] customized dataset  (0) 2025.07.01
GSE86982 Analysis Summary  (0) 2025.04.30
[ProtoCell 4P] scRNA Analysis  (0) 2025.04.10
[Hierarchical MIL] Preprocessing Create '.h5ad'  (0) 2025.03.28