单细胞RNA-ATAC多组学测序数据分析教程¶

# 加载必要的R包
library(Seurat)      # 单细胞分析核心包
library(Signac)      # 单细胞ATAC分析包
library(EnsDb.Hsapiens.v86)  # 人类基因组注释
library(BSgenome.Hsapiens.UCSC.hg38)  # 人类参考基因组
library(Pando)       # 基因调控网络重建包

# 设置工作目录
setwd("path/to/your/data")

# 读取10x Genomics输出的多组学数据
# 注意：数据包含三个文件：过滤后的特征条形码矩阵、片段文件和峰值文件
input_dir <- "path/to/10x_output"
counts <- Read10X_h5(file.path(input_dir, "filtered_feature_bc_matrix.h5"))
fragments <- file.path(input_dir, "atac_fragments.tsv.gz")
peaks <- file.path(input_dir, "atac_peaks.bed")

# 创建Seurat对象
# 注释：RNA和ATAC数据分别存储在Assay中
seurat_obj <- CreateSeuratObject(
  counts = counts$`Gene Expression`,
  assay = "RNA",
  project = "scMultiome"
)

# 添加ATAC数据
gr_peaks <- StringToGRanges(peaks, sep = c(":", "-"))
seurat_obj[["ATAC"]] <- CreateChromatinAssay(
  counts = counts$Peaks,
  sep = c(":", "-"),
  fragments = fragments,
  genome = "hg38"
)

# 添加基因注释信息
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) <- "UCSC"
Annotation(seurat_obj[["ATAC"]]) <- annotations

# RNA数据预处理
# 注释：使用SCTransform进行标准化，比LogNormalize效果更好
seurat_obj <- SCTransform(seurat_obj, assay = "RNA", verbose = FALSE)

# 降维与聚类
seurat_obj <- RunPCA(seurat_obj, assay = "SCT", verbose = FALSE)
seurat_obj <- FindNeighbors(seurat_obj, reduction = "pca", dims = 1:30)
seurat_obj <- FindClusters(seurat_obj, resolution = 0.5)
seurat_obj <- RunUMAP(seurat_obj, reduction = "pca", dims = 1:30)

# 可视化RNA聚类结果
DimPlot(seurat_obj, group.by = "seurat_clusters", label = TRUE) + 
  ggtitle("RNA-based Clustering")

# 鉴定差异表达基因
# 注释：比较每个cluster与其他所有cluster
rna_markers <- FindAllMarkers(seurat_obj, assay = "SCT", only.pos = TRUE)
head(rna_markers)

# ATAC数据预处理
# 注释：使用TF-IDF标准化，适合稀疏的ATAC数据
seurat_obj <- RunTFIDF(seurat_obj, assay = "ATAC")

# 使用LSI降维
seurat_obj <- FindTopFeatures(seurat_obj, assay = "ATAC", min.cutoff = "q0")
seurat_obj <- RunSVD(seurat_obj, assay = "ATAC")

# 去除与测序深度相关的成分（通常为第1个成分）
seurat_obj <- RunUMAP(seurat_obj, reduction = "lsi", dims = 2:30)

# 可视化ATAC聚类结果
DimPlot(seurat_obj, group.by = "seurat_clusters", label = TRUE) + 
  ggtitle("ATAC-based Clustering")

# 鉴定差异可及性区域
# 注释：使用Logistic回归或Wilcoxon检验
da_peaks <- FindAllMarkers(
  seurat_obj,
  assay = "ATAC",
  only.pos = TRUE,
  test.use = "LR",
  latent.vars = "nCount_ATAC"
)

# 计算加权最近邻（WNN）
# 注释：WNN自动为每个细胞分配RNA和ATAC的权重
seurat_obj <- FindMultiModalNeighbors(
  seurat_obj,
  reduction.list = list("pca", "lsi"),
  dims.list = list(1:30, 2:30),
  modality.weight.name = "RNA.weight"
)

# 基于WNN进行聚类和UMAP
seurat_obj <- FindClusters(seurat_obj, graph.name = "wsnn", resolution = 0.5)
seurat_obj <- RunUMAP(seurat_obj, nn.name = "weighted.nn", reduction.name = "wnn.umap")

# 可视化整合结果
DimPlot(seurat_obj, reduction = "wnn.umap", group.by = "seurat_clusters", label = TRUE) +
  ggtitle("WNN-based Clustering")

# 比较不同模态的聚类一致性
# 注释：计算调整兰德指数（ARI）
library(mclust)
ari_rna <- adjustedRandIndex(
  seurat_obj$seurat_clusters,
  seurat_obj$wsnn_clusters
)
cat("RNA vs WNN ARI:", ari_rna)

# 准备Pando分析所需数据
# 注释：Pando需要基因表达矩阵和染色质开放性矩阵
seurat_obj <- initiate_grn(
  seurat_obj,
  rna_assay = "SCT",
  atac_assay = "ATAC",
  genome = "hg38"
)

# 扫描转录因子结合位点
# 注释：使用JASPAR数据库中的motif
seurat_obj <- find_motifs(
  seurat_obj,
  pfm = motifs,
  genome = BSgenome.Hsapiens.UCSC.hg38
)

# 构建调控模型
# 注释：使用负二项回归拟合基因表达
seurat_obj <- infer_grn(
  seurat_obj,
  peak_to_gene_method = "signac",
  method = "glm",
  family = "negative.binomial"
)

# 提取显著调控关系
# 注释：筛选p值小于0.05的调控对
grn_network <- Network(seurat_obj)
significant_edges <- grn_network[grn_network$pval < 0.05, ]
head(significant_edges)

# 可视化特定转录因子的调控网络
# 注释：例如查看转录因子SOX9的靶基因
plot_network(seurat_obj, tf = "SOX9")