DeepSeek，让代码更丝滑！GSE98177，多样本、多分组，二代测序数据！如何提问是个技术活～

########--------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 ( ) 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)     } elseif (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: 进行功能富集分析...