单细胞表观基因组分析方法¶

# ===== R: scATAC-seq 分析（使用 Signac + Seurat） =====
library(Signac)   # 专门处理 scATAC-seq 的 R 包
library(Seurat)   # 单细胞分析主力框架
library(GenomicRanges)

# 1. 读取 10x Genomics scATAC-seq 数据
# 需要三个文件：fragments.tsv.gz、peaks.bed、metadata
counts <- Read10X_h5("filtered_peak_bc_matrix.h5")
metadata <- read.csv("singlecell.csv", row.names = 1)

# 创建 ChromatinAssay 对象
chrom_assay <- CreateChromatinAssay(
    counts = counts,
    sep = c(":", "-"),          # peak 名称分隔符
    fragments = "fragments.tsv.gz",  # 片段文件
    min.cells = 10,            # 至少 10 个细胞有信号
    min.features = 200)        # 每个细胞至少 200 个 peaks

# 创建 Seurat 对象
pbmc <- CreateSeuratObject(
    counts = chrom_assay,
    assay = "peaks",
    meta.data = metadata)

# 2. 质控
# 计算 QC 指标
pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments /
    pbmc$passed_filters * 100  # peaks 内 reads 比例
pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments /
    pbmc$peak_region_fragments  # 黑名单区域比例

# 过滤低质量细胞
pbmc <- subset(pbmc,
    peak_region_fragments > 3000 &    # 至少 3000 fragments
    peak_region_fragments < 20000 &   # 不超过 20000（排除多细胞）
    pct_reads_in_peaks > 15 &         # peaks 内比例 > 15%
    blacklist_ratio < 0.05 &          # 黑名单比例 < 5%
    nucleosome_signal < 4 &           # 核小体信号合理
    TSS.enrichment > 3)               # TSS 富集 > 3

# 3. 标准化 + 降维
# TF-IDF 标准化（Term Frequency - Inverse Document Frequency）
pbmc <- RunTFIDF(pbmc)
# 选择高变特征
pbmc <- FindTopFeatures(pbmc, min.cutoff = "q0")
# LSI 降维（Latent Semantic Indexing）
pbmc <- RunSVD(pbmc)
# 注意：第一个 LSI 分量通常与测序深度相关，需检查并排除

# 4. UMAP 可视化 + 聚类
pbmc <- RunUMAP(pbmc, reduction = "lsi", dims = 2:30)  # 跳过第1维
pbmc <- FindNeighbors(pbmc, reduction = "lsi", dims = 2:30)
pbmc <- FindClusters(pbmc, resolution = 0.5,
                     algorithm = 3)  # SLM 算法

# 5. 计算基因活性评分（Gene Activity Score）
# 用基因体 + 上游 2kb 启动子区域的可及性估计表达
gene_activities <- GeneActivity(pbmc)
pbmc[["RNA"]] <- CreateAssayObject(counts = gene_activities)
pbmc <- NormalizeData(pbmc, assay = "RNA",
    normalization.method = "LogNormalize",
    scale.factor = median(pbmc$nCount_RNA))

# 6. 可视化标记基因
DefaultAssay(pbmc) <- "RNA"
# 用已知标记基因辅助标注细胞类型
FeaturePlot(pbmc,
    features = c("MS4A1",   # B 细胞
                 "CD3D",    # T 细胞
                 "CD14",    # 单核细胞
                 "FCGR3A"), # NK 细胞
    max.cutoff = "q95")