########--------GSE98177二代测序数据的分析--------########
# GSE98177数据分析：识别调控IL2RA的增强子
# 1. 加载必要包 --------------------------------------------------------------
library(GEOquery)
## 载入需要的程序包：Biobase
## 载入需要的程序包：BiocGenerics
##
## 载入程序包：'BiocGenerics'
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, table,
## tapply, union, unique, unsplit, which.max, which.min
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
## Setting options('download.file.method.GEOquery'='auto')
## Setting options('GEOquery.inmemory.gpl'=FALSE)
library(DESeq2)
## 载入需要的程序包：S4Vectors
## Warning: 程序包'S4Vectors'是用R版本4.4.1 来建造的
## 载入需要的程序包：stats4
##
## 载入程序包：'S4Vectors'
## The following object is masked from 'package:utils':
##
## findMatches
## The following objects are masked from 'package:base':
##
## expand.grid, I, unname
## 载入需要的程序包：IRanges
## Warning: 程序包'IRanges'是用R版本4.4.1 来建造的
##
## 载入程序包：'IRanges'
## The following object is masked from 'package:grDevices':
##
## windows
## 载入需要的程序包：GenomicRanges
## Warning: 程序包'GenomicRanges'是用R版本4.4.1 来建造的
## 载入需要的程序包：GenomeInfoDb
## 载入需要的程序包：SummarizedExperiment
## 载入需要的程序包：MatrixGenerics
## 载入需要的程序包：matrixStats
##
## 载入程序包：'matrixStats'
## The following objects are masked from 'package:Biobase':
##
## anyMissing, rowMedians
##
## 载入程序包：'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars
## The following object is masked from 'package:Biobase':
##
## rowMedians
library(tidyverse)
## Warning: 程序包'tidyverse'是用R版本4.4.2 来建造的
## Warning: 程序包'lubridate'是用R版本4.4.2 来建造的
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ lubridate::%within%() masks IRanges::%within%()
## ✖ dplyr::collapse() masks IRanges::collapse()
## ✖ dplyr::combine() masks Biobase::combine(), BiocGenerics::combine()
## ✖ dplyr::count() masks matrixStats::count()
## ✖ dplyr::desc() masks IRanges::desc()
## ✖ tidyr::expand() masks S4Vectors::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first() masks S4Vectors::first()
## ✖ dplyr::lag() masks stats::lag()
## ✖ ggplot2::Position() masks BiocGenerics::Position(), base::Position()
## ✖ purrr::reduce() masks GenomicRanges::reduce(), IRanges::reduce()
## ✖ dplyr::rename() masks S4Vectors::rename()
## ✖ lubridate::second() masks S4Vectors::second()
## ✖ lubridate::second<-() masks S4Vectors::second<-()
## ✖ dplyr::slice() masks IRanges::slice()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2)
library(pheatmap)
library(clusterProfiler)
## Warning: 程序包'clusterProfiler'是用R版本4.4.1 来建造的
##
## clusterProfiler v4.12.6 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
##
## Please cite:
##
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan,
## X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal
## enrichment tool for interpreting omics data. The Innovation. 2021,
## 2(3):100141
##
## 载入程序包：'clusterProfiler'
##
## The following object is masked from 'package:purrr':
##
## simplify
##
## The following object is masked from 'package:IRanges':
##
## slice
##
## The following object is masked from 'package:S4Vectors':
##
## rename
##
## The following object is masked from 'package:stats':
##
## filter
library(org.Hs.eg.db)
## 载入需要的程序包：AnnotationDbi
##
## 载入程序包：'AnnotationDbi'
##
## The following object is masked from 'package:clusterProfiler':
##
## select
##
## The following object is masked from 'package:dplyr':
##
## select
# 2. 数据下载与加载 ----------------------------------------------------------
cat("Step 1: 下载并加载GSE98177数据...\n")
## Step 1: 下载并加载GSE98177数据...
gse <- getGEO("GSE98177", GSEMatrix = TRUE)
## Found 1 file(s)
## GSE98177_series_matrix.txt.gz
# 获取表达矩阵和样本信息
expr_data <- exprs(gse[[1]]) # 标准化表达矩阵(如log2FPKM)
pdata <- pData(gse[[1]]) # 样本元数据
count_data <- read.table("GSE98177_raw_counts_tsv.gz", header=TRUE)
count_data <- count_data %>% rename(ENTREZID = GeneID)
count_data <- count_data %>% mutate(ENTREZID = as.character(ENTREZID))
# 获取所有可能的映射
all_mappings <- select(org.Hs.eg.db,
keys = as.character(count_data$ENTREZID),
columns = "SYMBOL",
keytype = "ENTREZID")
## 'select()' returned 1:1 mapping between keys and columns
# 处理一对多关系 - 选择表达量最高的那个样本的符号
library(dplyr)
final_mapping <- count_data %>%
left_join(all_mappings, by = "ENTREZID") %>%
group_by(ENTREZID) %>%
mutate(SYMBOL = ifelse(all(is.na(SYMBOL)),
ENTREZID,
na.omit(SYMBOL)[1])) %>%
# 添加计数使SYMBOL唯一
group_by(SYMBOL) %>%
mutate(dup_count = seq_along(SYMBOL),
UNIQUE_SYMBOL = ifelse(dup_count > 1,
paste0(SYMBOL, "_", dup_count),
as.character(SYMBOL))) %>%
ungroup() %>%
distinct(ENTREZID, .keep_all = TRUE)
data = as.data.frame(final_mapping)
row.names(data) = data$UNIQUE_SYMBOL # 使用唯一的符号作为行名
# 移除临时列
dat = data[,c(2:19)]
# 3. 样本分组设计 -----------------------------------------------------------
cat("Step 2: 设置样本分组...\n")
## Step 2: 设置样本分组...
pdata$condition <- factor(c(
rep("NonTargeting", 2), rep("TSS_sg1", 2), rep("TSS_sg2", 2),
rep("CaRE3_sg1", 2), rep("CaRE3_sg2", 2), rep("CaRE4_sg1", 2),
rep("CaRE4_sg2", 2), rep("PBS", 2), rep("Stimulated", 2)
))
pdata$condition <- factor(pdata$condition,
levels = c("NonTargeting", "TSS_sg1", "TSS_sg2",
"CaRE3_sg1", "CaRE3_sg2", "CaRE4_sg1",
"CaRE4_sg2", "PBS", "Stimulated"))
# 4. 差异表达分析 ----------------------------------------------------------
cat("Step 3: 进行差异表达分析...\n")
## Step 3: 进行差异表达分析...
# 差异分析
dds <- DESeqDataSetFromMatrix(
countData = dat,
colData = pdata,
design = ~ condition
)
# 运行DESeq2
dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
# 比较组分析
res_care3_sg1 <- results(dds, contrast = c("condition", "CaRE3_sg1", "NonTargeting"))
res_care3_sg2 <- results(dds, contrast = c("condition", "CaRE3_sg2", "NonTargeting"))
res_care4_sg1 <- results(dds, contrast = c("condition", "CaRE4_sg1", "NonTargeting"))
# 5. IL2RA表达分析 ---------------------------------------------------------
cat("Step 4: 分析IL2RA表达模式...\n")
## Step 4: 分析IL2RA表达模式...
il2ra_id <- "IL2RA" # 假设IL2RA存在于数据中
# 提取归一化表达数据
il2ra_expr <- log2(counts(dds, normalized=TRUE)[il2ra_id, ] + 1)
il2ra_df <- data.frame(
expr = as.numeric(il2ra_expr),
condition = pdata$condition
)
# 绘制表达箱线图
ggplot(il2ra_df, aes(x = condition, y = expr, fill = condition)) +
geom_boxplot() +
geom_jitter(width = 0.2) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(
title = "IL2RA Expression Across Conditions",
y = "Normalized Counts (log2)",
x = ""
)
# 6. 差异基因热图 ----------------------------------------------------------
cat("Step 5: 生成差异基因热图...\n")
## Step 5: 生成差异基因热图...
# 获取显著差异基因(padj < 0.05 & |log2FC| > 1)
sig_genes <- rownames(res_care3_sg1)[
which(res_care3_sg1$padj < 0.05 & abs(res_care3_sg1$log2FoldChange) > 1)
]
# 确保sig_genes存在于vsd数据中
vsd <- vst(dds, blind=FALSE)
sig_genes <- sig_genes[sig_genes %in% rownames(assay(vsd))]
# 如果找到显著基因则绘制热图
if(length(sig_genes) > 0){
vsd <- vst(dds, blind=FALSE)
mat <- assay(vsd)[sig_genes, , drop=FALSE]
# 检查mat是否有有效数据
if(all(is.finite(mat))) {
# 对数据进行缩放
mat_scaled <- t(scale(t(mat)))
# 设置热图颜色
heatmap_colors <- colorRampPalette(c("blue", "white", "red"))(100)
# 只有当基因和样本数都≥2时才绘制聚类热图
if(nrow(mat_scaled) >= 2 && ncol(mat_scaled) >= 2){
pheatmap(mat_scaled,
annotation_col = pdata["condition"],
show_rownames = ifelse(length(sig_genes) < 50, TRUE, FALSE),
scale = "none",
color = heatmap_colors,
main = "Differentially Expressed Genes in CaRE3 vs Control",
cluster_rows = TRUE,
cluster_cols = TRUE,
fontsize_row = 8,
fontsize_col = 8)
} else if (nrow(mat_scaled) >= 1) {
# 如果只有1个基因但有足够样本，绘制非聚类热图
cat("注意: 只有", nrow(mat_scaled), "个差异基因 - 绘制简化热图\n")
pheatmap(mat_scaled,
annotation_col = pdata["condition"],
show_rownames = TRUE,
scale = "none",
color = heatmap_colors,
main = paste("Single DEG:", rownames(mat_scaled)),
cluster_rows = FALSE,
cluster_cols = FALSE,
fontsize_row = 8,
fontsize_col = 8)
}
} else {
cat("警告: 热图数据包含非有限值(NA/NaN/Inf)\n")
problematic_genes <- apply(mat, 1, function(x) any(!is.finite(x)))
print(names(which(problematic_genes)))
}
} else {
cat("未检测到满足条件的显著差异基因\n")
}
## 注意: 只有 1 个差异基因 - 绘制简化热图
# 7. 功能富集分析 ----------------------------------------------------------
cat("Step 6: 进行功能富集分析...\n")
## Step 6: 进行功能富集分析...
# 示例基因列表（实际应使用差异分析结果）
example_genes <- c("IL2RA", "CD28", "CD3D", "CD4", "FOXP3", "STAT5A")
# 转换基因ID
entrez_ids <- mapIds(org.Hs.eg.db,
keys = example_genes,
column = "ENTREZID",
keytype = "SYMBOL")
## 'select()' returned 1:1 mapping between keys and columns
# GO富集分析
go_enrich <- enrichGO(
gene = na.omit(entrez_ids),
OrgDb = org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "BH",
readable = TRUE
)
# 绘制富集结果
if(nrow(go_enrich) > 0){
dotplot(go_enrich, title = "GO Enrichment Analysis")
} else {
cat("未检测到显著富集的GO terms\n")
}
# 8. 结果汇总 --------------------------------------------------------------
cat("\n===== 分析完成 =====\n")
##
## ===== 分析完成 =====
cat("关键结论检查:\n")
## 关键结论检查:
cat("1. 检查IL2RA在各处理组中的表达变化\n")
## 1. 检查IL2RA在各处理组中的表达变化
cat("2. 验证差异基因是否富集于T细胞激活相关通路\n")
## 2. 验证差异基因是否富集于T细胞激活相关通路
cat("3. 比较CaRE3和CaRE4处理组的差异基因重叠情况\n")
## 3. 比较CaRE3和CaRE4处理组的差异基因重叠情况
cat("Step 4: 分析IL2RA表达模式...\n")
## Step 4: 分析IL2RA表达模式...
cat("Step 5: 生成差异基因热图...\n")
## Step 5: 生成差异基因热图...
cat("Step 6: 进行功能富集分析...\n")
## Step 6: 进行功能富集分析...