70. R语言 Bioconductor 入门¶

# ===== 第一步：安装BiocManager（只需一次）=====
install.packages("BiocManager")

# ===== 检查Bioconductor版本 =====
BiocManager::version()
# [1] '3.23'

# ===== 安装单个包 =====
BiocManager::install("DESeq2")              # 安装差异表达分析包

# ===== 安装多个包 =====
BiocManager::install(c(
  "DESeq2",                                 # 差异表达
  "edgeR",                                  # 差异表达（另一种方法）
  "clusterProfiler",                        # 富集分析
  "GenomicRanges",                          # 基因组坐标操作
  "Biostrings",                             # 序列操作
  "AnnotationDbi",                          # 注释数据库接口
  "org.Hs.eg.db"                            # 人类基因注释数据库
))

# ===== 更新所有包 =====
BiocManager::install()                       # 不加参数 = 更新所有已安装的包

# ===== 检查安装是否正确 =====
BiocManager::valid()                         # 检查版本兼容性

# ===== 查看已安装的Bioconductor包 =====
BiocManager::installed()

# 问题1：版本不匹配
# 报错："Bioconductor version '3.23' requires R version '4.6'"
# 解决：升级R到对应版本，或安装旧版Bioconductor
BiocManager::install(version = "3.21")       # 安装旧版（配合旧R）

# 问题2：依赖编译失败（Linux）
# 解决：安装系统依赖
# sudo apt install libcurl4-openssl-dev libxml2-dev libssl-dev

# 问题3：国内下载慢
# 解决：设置镜像
options(BioC_mirror = "https://mirrors.tuna.tsinghua.edu.cn/bioconductor")

library(SummarizedExperiment)

# ===== 创建一个SummarizedExperiment =====
# 模拟表达数据：3个基因 × 4个样本
counts_matrix <- matrix(
  c(100, 200, 50, 150,      # 基因1在4个样本中的表达量
    300, 100, 400, 200,      # 基因2
    50,  80,  30,  60),      # 基因3
  nrow = 3, byrow = TRUE,
  dimnames = list(
    c("GeneA", "GeneB", "GeneC"),             # 行名 = 基因名
    c("Sample1", "Sample2", "Sample3", "Sample4")  # 列名 = 样本名
  )
)

# 样本信息（colData）
sample_info <- DataFrame(
  condition = c("control", "control", "disease", "disease"),  # 分组
  age = c(25, 30, 28, 35)                                     # 年龄
)

# 基因信息（rowData）
gene_info <- DataFrame(
  gene_id = c("ENSG001", "ENSG002", "ENSG003"),
  gene_name = c("GeneA", "GeneB", "GeneC"),
  chromosome = c("chr1", "chr2", "chr3")
)

# 组装成SummarizedExperiment对象
se <- SummarizedExperiment(
  assays = list(counts = counts_matrix),       # 表达矩阵
  colData = sample_info,                       # 样本信息
  rowData = gene_info                          # 基因信息
)

# ===== 访问数据 =====
assay(se, "counts")        # 获取counts矩阵
colData(se)                # 获取样本信息
rowData(se)                # 获取基因信息
se$condition               # 快捷访问样本分组

library(GenomicRanges)

# ===== 创建GRanges对象 =====
gr <- GRanges(
  seqnames = c("chr1", "chr1", "chr2", "chr3"),  # 染色体
  ranges = IRanges(                                # 坐标范围
    start = c(100, 200, 300, 400),                 # 起始位置
    end   = c(150, 250, 350, 500)                  # 终止位置
  ),
  strand = c("+", "-", "+", "*"),                  # 链方向（*=不确定）
  gene_name = c("GeneA", "GeneB", "GeneC", "GeneD"),  # 自定义列
  score = c(10, 20, 15, 30)                            # 自定义列
)

print(gr)
# GRanges object with 4 ranges and 2 metadata columns:
#       seqnames    ranges strand | gene_name     score
#   [1]     chr1   100-150      + |     GeneA        10
#   [2]     chr1   200-250      - |     GeneB        20
#   [3]     chr2   300-350      + |     GeneC        15
#   [4]     chr3   400-500      * |     GeneD        30

# ===== 区间操作（核心功能！）=====
# 找重叠区间
query <- GRanges("chr1", IRanges(120, 220))        # 查询区间
overlaps <- findOverlaps(query, gr)                 # 找与gr重叠的区间
gr[subjectHits(overlaps)]                           # 返回重叠的区间

# 取交集
intersect_gr <- GenomicRanges::intersect(gr[1], gr[2])

# 区间扩展（上下游各扩展100bp）
gr_extended <- resize(gr, width = width(gr) + 200, fix = "center")

# 取子集
gr_chr1 <- gr[seqnames(gr) == "chr1"]              # 只取chr1的区间
gr_positive <- gr[strand(gr) == "+"]                # 只取正链

library(IRanges)

# ===== IRanges是GRanges的基础（不含染色体和链信息）=====
ir <- IRanges(start = c(1, 5, 10), end = c(4, 8, 15))
print(ir)
# IRanges object with 3 ranges:
#       start end width
#   [1]     1   4     4
#   [2]     5   8     4
#   [3]    10  15     6

# 区间合并（合并重叠的区间）
ir2 <- IRanges(start = c(1, 3, 10), end = c(5, 7, 15))
reduce(ir2)                                         # 合并重叠 → (1,7) (10,15)

# 计算覆盖度
coverage(ir)                                        # 每个位置被几个区间覆盖

library(Biostrings)

# ===== 创建DNA序列 =====
dna <- DNAString("ATGATCTCGTAA")                    # 单条DNA序列
print(dna)
# 12-letter DNAString object
# seq: ATGATCTCGTAA

# ===== 基本操作 =====
complement(dna)                                     # 互补链：TACTAGAGCATT
reverseComplement(dna)                              # 反向互补：TTACGAGATCAT
translate(dna)                                      # 翻译：MIS*
alphabetFrequency(dna)                              # 碱基频率统计
letterFrequency(dna, "GC") / length(dna)            # GC含量

# ===== 序列集合 =====
dna_set <- DNAStringSet(c(
  gene1 = "ATGATCTCGTAA",
  gene2 = "GCGCATATGGCC",
  gene3 = "TTAAGCGCGATA"
))
width(dna_set)                                      # 各序列长度
letterFrequency(dna_set, "GC") / width(dna_set)     # 批量计算GC含量

# ===== 模式匹配 =====
matchPattern("ATG", dna)                            # 查找ATG在序列中的位置
vmatchPattern("ATG", dna_set)                       # 在序列集合中查找

library(DESeq2)

# ===== 从counts矩阵创建DESeq2对象 =====
dds <- DESeqDataSetFromMatrix(
  countData = counts_matrix,                        # 基因×样本的counts矩阵
  colData = sample_info,                            # 样本信息（必须包含分组列）
  design = ~ condition                              # 差异分析的公式：按condition分组
)

# ===== 运行差异分析 =====
dds <- DESeq(dds)                                   # 一行搞定（标准化+统计检验）

# ===== 提取结果 =====
res <- results(dds, contrast = c("condition", "disease", "control"))
res_df <- as.data.frame(res)                        # 转为数据框方便查看

# ===== 查看显著差异基因 =====
sig_genes <- subset(res_df, padj < 0.05 & abs(log2FoldChange) > 1)
print(sig_genes)
# 列说明：
#   baseMean      —— 所有样本的平均表达量
#   log2FoldChange —— log2(disease/control)，正值=上调，负值=下调
#   padj           —— 校正后p值（BH方法）

library(edgeR)

# ===== 创建DGEList对象 =====
group <- factor(c("control", "control", "disease", "disease"))
dge <- DGEList(counts = counts_matrix, group = group)

# ===== 标准化 + 差异分析 =====
dge <- calcNormFactors(dge)                         # TMM标准化
design <- model.matrix(~ group)                     # 设计矩阵
dge <- estimateDisp(dge, design)                    # 估计离散度
fit <- glmQLFit(dge, design)                        # 拟合模型
res <- glmQLFTest(fit, coef = 2)                    # 差异检验
topTags(res, n = 10)                                # 查看top 10差异基因

library(clusterProfiler)
library(org.Hs.eg.db)                               # 人类基因注释数据库

# ===== 准备基因列表（Entrez ID）=====
gene_list <- c("7157", "672", "675", "5290", "1956") # 示例：TP53, BRCA1, BRCA2等

# ===== GO富集分析 =====
go_result <- enrichGO(
  gene = gene_list,
  OrgDb = org.Hs.eg.db,                             # 物种注释库
  ont = "BP",                                        # BP/MF/CC
  pAdjustMethod = "BH",                              # p值校正方法
  pvalueCutoff = 0.05,
  readable = TRUE                                    # 将ID转为基因名
)
head(go_result)

# ===== KEGG富集分析 =====
kegg_result <- enrichKEGG(
  gene = gene_list,
  organism = "hsa",                                  # hsa=人类
  pvalueCutoff = 0.05
)
head(kegg_result)

# ===== 可视化 =====
dotplot(go_result, showCategory = 10)                # 气泡图
barplot(go_result, showCategory = 10)                # 柱状图

library(GenomicRanges)

# 详见上面3.2节的GRanges代码示例
# 核心功能：findOverlaps, intersect, reduce, coverage
# 常用于ChIP-seq peak分析、变异注释、区间比较

library(AnnotationDbi)
library(org.Hs.eg.db)

# ===== 基因ID转换（最常用！）=====
# Entrez ID → Gene Symbol
symbols <- mapIds(
  org.Hs.eg.db,
  keys = c("7157", "672", "675"),                    # 输入ID
  column = "SYMBOL",                                 # 目标格式
  keytype = "ENTREZID",                              # 输入格式
  multiVals = "first"                                # 多个结果取第一个
)
print(symbols)
# 7157  672  675
# "TP53" "BRCA1" "BRCA2"

# Gene Symbol → Ensembl ID
ensembl_ids <- mapIds(
  org.Hs.eg.db,
  keys = c("TP53", "BRCA1"),
  column = "ENSEMBL",
  keytype = "SYMBOL"
)

# 查看支持的ID类型
columns(org.Hs.eg.db)
# "ENTREZID" "ENSEMBL" "SYMBOL" "GENENAME" "GO" "PATH" "UNIPROT" ...

library(ggplot2)

# 假设res_df是DESeq2的结果数据框，含log2FoldChange和padj列
# res_df <- as.data.frame(results(dds))

# ===== 添加颜色分类列 =====
res_df$significance <- "Not Sig"                                     # 默认不显著
res_df$significance[res_df$padj < 0.05 & res_df$log2FoldChange > 1] <- "Up"    # 上调
res_df$significance[res_df$padj < 0.05 & res_df$log2FoldChange < -1] <- "Down" # 下调

# ===== 画火山图 =====
ggplot(res_df, aes(x = log2FoldChange, y = -log10(padj), color = significance)) +
  geom_point(alpha = 0.6, size = 1.5) +                             # 散点
  scale_color_manual(values = c("Up" = "red", "Down" = "blue", "Not Sig" = "grey")) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +      # p值阈值线
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +          # FC阈值线
  labs(
    title = "Volcano Plot",
    x = "log2(Fold Change)",
    y = "-log10(adjusted p-value)"
  ) +
  theme_bw() +                                                      # 白色背景
  theme(legend.position = "top")

library(pheatmap)    # 或 ComplexHeatmap

# ===== 用pheatmap画热图 =====
# 假设 expr_matrix 是标准化后的表达矩阵（行=基因，列=样本）

# 准备注释信息
annotation_col <- data.frame(
  Group = c("Control", "Control", "Disease", "Disease"),
  row.names = colnames(counts_matrix)
)

pheatmap(
  counts_matrix,                                     # 表达矩阵
  scale = "row",                                     # 按行标准化（Z-score）
  clustering_method = "ward.D2",                     # 聚类方法
  annotation_col = annotation_col,                   # 列注释（样本分组）
  color = colorRampPalette(c("blue", "white", "red"))(100),  # 颜色
  show_rownames = TRUE,                              # 显示基因名
  fontsize_row = 8,                                  # 基因名字体大小
  main = "Gene Expression Heatmap"                   # 标题
)

library(ggplot2)

# MA图：x轴=平均表达量（A），y轴=差异倍数（M）
# res_df需要有baseMean和log2FoldChange列

ggplot(res_df, aes(x = log10(baseMean + 1), y = log2FoldChange, color = significance)) +
  geom_point(alpha = 0.5, size = 1) +
  scale_color_manual(values = c("Up" = "red", "Down" = "blue", "Not Sig" = "grey")) +
  geom_hline(yintercept = 0, color = "black") +                     # 零线
  geom_hline(yintercept = c(-1, 1), linetype = "dashed", color = "grey50") +
  labs(
    title = "MA Plot",
    x = "log10(Mean Expression)",
    y = "log2(Fold Change)"
  ) +
  theme_bw()

答：两者都是基于负二项分布的RNA-seq差异表达分析工具。
- DESeq2：用中值比率法（median of ratios）标准化，对小样本量（<5）表现更稳健，
  自动做离群值检测和独立过滤。代码更简洁（一个DESeq()搞定）。
- edgeR：用TMM标准化，计算速度稍快，灵活性更高（支持更复杂的实验设计）。
- 实际差异不大，大多数情况下结果高度一致。推荐：常规分析用DESeq2，
  复杂设计（多因素、配对样本）两个都试试看结果是否一致。

答：直接用matrix的问题是：矩阵只有数值，样本信息和基因信息要另外维护，
容易行列不对齐。SummarizedExperiment把三者打包在一起：
- assay：表达矩阵（counts/TPM/FPKM）
- colData：样本元数据（分组、批次、临床信息）
- rowData：基因元数据（基因名、坐标、功能注释）
取子集时三者自动同步，不会出错。DESeq2/edgeR等包也都以它为输入格式。

答：GRanges是GenomicRanges包定义的基因组区间对象，存储"染色体+起止坐标+链方向"。
常用场景：
1. ChIP-seq peak文件（BED格式）的导入和操作
2. 基因组注释区间的重叠比较（findOverlaps）
3. 变异位点（SNP/InDel）落在哪些基因区域
4. 合并重叠区间（reduce）、计算覆盖度（coverage）

答：核心是超几何检验（Fisher精确检验的推广）。
给定一个差异基因列表和一个功能集合（如GO term包含哪些基因），
计算"差异基因中属于这个GO term的比例"是否显著高于"背景基因中的比例"。
p值用BH方法校正为padj，padj < 0.05的term被认为显著富集。
GSEA则不需要设阈值，用排序后的全基因列表做富集打分。

答：用AnnotationDbi包的mapIds()函数，配合物种注释库（如org.Hs.eg.db）。
支持的ID类型包括：ENTREZID、ENSEMBL、SYMBOL、UNIPROT、GO等。
示例：mapIds(org.Hs.eg.db, keys="TP53", column="ENTREZID", keytype="SYMBOL")
注意：有些基因可能没有对应ID，需要处理NA值。
多对一映射时，用multiVals参数控制取第一个还是全部返回。