染色质可及性分析¶

ATAC-seq的工作原理：

1. 收集细胞（500-50000个细胞）
2. 温和裂解细胞膜（保留细胞核完整）
3. Tn5转座酶处理：
   Tn5携带测序接头 → 插入开放染色质区域
   紧密包装的染色质（异染色质）→ Tn5无法插入
   松散的染色质（真染色质）→ Tn5插入并标记
4. 提取DNA → PCR扩增 → 测序

reads长度分布特征（质控关键）：
  <150bp: 无核小体区域（nucleosome-free, NFR）→ 最有价值
  ~200bp: 单核小体保护片段（mono-nucleosome）
  ~400bp: 双核小体保护片段
  ~600bp: 三核小体保护片段

信号特征：
  - TSS（转录起始位点）附近信号高 → 活跃转录基因的启动子开放
  - 增强子区域信号高 → 远端调控元件开放
  - 异染色质区域信号低 → 基因沉默

Bulk ATAC-seq:
  - 输入：数千到数万个细胞的混合物
  - 输出：群体平均的染色质开放谱
  - 分析工具：MACS2 peak calling
  - 适合：比较不同条件下的染色质变化

scATAC-seq（单细胞ATAC-seq）:
  - 输入：单个细胞（10x Genomics平台）
  - 输出：每个细胞的染色质开放谱
  - 分析工具：ArchR / Signac
  - 适合：细胞异质性分析、细胞类型特异的调控元件

挑战：scATAC-seq数据极度稀疏
  每个细胞只捕获数千个插入位点（vs 基因组3×10^9 bp）
  → 需要特殊的降维方法（LSI而非PCA）
  → 需要特殊的peak calling策略

# 1. 质控
fastqc sample_R1.fq.gz sample_R2.fq.gz  # 质控报告

# 2. 去接头（Nextera接头）
trim_galore --paired \
  --nextera \                  # ATAC-seq用Nextera接头
  sample_R1.fq.gz \
  sample_R2.fq.gz

# 3. 比对到基因组（Bowtie2推荐）
bowtie2 \
  -x hg38_index \              # Bowtie2索引
  -1 sample_R1_trimmed.fq.gz \ # 正向reads
  -2 sample_R2_trimmed.fq.gz \ # 反向reads
  --very-sensitive \            # 高灵敏度模式
  --no-mixed \                  # 不允许单端比对
  --no-discordant \             # 不允许不一致比对
  -X 1000 \                    # 最大插入片段1000bp
  -p 8 \                       # 8个线程
  -S sample.sam                # 输出SAM

# 4. 后处理
samtools view -bS -q 30 -F 1804 sample.sam | \  # 过滤低质量和异常reads
  samtools sort -o sample_sorted.bam              # 排序
samtools index sample_sorted.bam                  # 索引

# 5. 去除PCR重复
picard MarkDuplicates \
  INPUT=sample_sorted.bam \
  OUTPUT=sample_dedup.bam \
  METRICS_FILE=dup_metrics.txt \
  REMOVE_DUPLICATES=true        # 去除重复
samtools index sample_dedup.bam

# 6. 去除线粒体reads（ATAC-seq中线粒体DNA比例常很高）
samtools view -h sample_dedup.bam | \
  grep -v "chrM" | \           # 去除线粒体
  samtools view -bS - > sample_final.bam
samtools index sample_final.bam

# 7. 偏移校正（Tn5插入偏移+4/-5bp）
# Tn5插入位点：正链+4bp，负链-5bp
alignmentSieve \
  --bam sample_final.bam \
  --ATACshift \                # ATAC-seq专用偏移校正
  --outFile sample_shifted.bam
samtools index sample_shifted.bam

# 用MACS2进行peak calling
macs2 callpeak \
  -t sample_shifted.bam \      # 处理后的BAM
  -f BAMPE \                   # 双端测序
  --nomodel \                  # 不建立模型
  --shift -75 \                # 偏移（配合extsize使用）
  --extsize 150 \              # 片段延伸
  -g hs \                      # 人类基因组大小
  -q 0.05 \                    # FDR阈值
  --keep-dup all \             # 保留所有reads（已去重）
  -n sample_atac \             # 输出前缀
  --outdir macs2_output/       # 输出目录

# 输出文件：
# sample_atac_peaks.narrowPeak — peak位置
# sample_atac_summits.bed — peak顶点

# ATAC-seq质控（R语言）
library(ATACseqQC)

# 1. 片段长度分布
fragSizeDist <- fragSizeDist(bamFiles = "sample_shifted.bam")
# 应该看到NFR峰(<150bp)和核小体阶梯（~200bp, ~400bp, ~600bp）

# 2. TSS富集
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)

peaks <- readPeakFile("sample_atac_peaks.narrowPeak")
promoter <- getPromoters(TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene,
                          upstream = 3000, downstream = 3000)
tagMatrix <- getTagMatrix(peaks, windows = promoter)  # 在TSS周围计数
tagHeatmap(tagMatrix, xlim = c(-3000, 3000),           # 热图
           title = "TSS Enrichment")

# 3. 基因组区域注释
peakAnno <- annotatePeak(peaks,
                          TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene)
plotAnnoPie(peakAnno)  # 饼图：启动子、内含子、基因间区等占比

# 安装Signac
install.packages("Signac")     # 安装
library(Signac)                 # 加载
library(Seurat)                 # 加载Seurat
library(GenomicRanges)

# 读取10x CellRanger ATAC输出
counts <- Read10X_h5("filtered_peak_bc_matrix.h5")  # peak×cell矩阵
metadata <- read.csv("singlecell.csv", row.names = 1)  # 元数据

# 创建ChromatinAssay
chrom_assay <- CreateChromatinAssay(
  counts = counts,
  sep = c(":", "-"),           # peak命名格式（chr:start-end）
  genome = "hg38",             # 基因组版本
  fragments = "fragments.tsv.gz",  # 片段文件
  min.cells = 10,              # peak至少在10个细胞中出现
  min.features = 200           # 细胞至少200个peak
)

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

# 质控过滤
seurat_obj <- subset(seurat_obj,
  peak_region_fragments > 3000 &    # 最少3000个片段
  peak_region_fragments < 100000 &  # 最多100000（排除双细胞）
  nucleosome_signal < 4 &           # 核小体信号<4
  TSS.enrichment > 2                # TSS富集>2
)

# 1. TF-IDF归一化（不用LogNormalize）
seurat_obj <- RunTFIDF(seurat_obj)  # Term Frequency-Inverse Document Frequency

# 2. 选择高变特征
seurat_obj <- FindTopFeatures(seurat_obj, min.cutoff = "q0")  # 所有特征

# 3. SVD降维（不用PCA）
seurat_obj <- RunSVD(seurat_obj)  # 奇异值分解（=LSI第二步）

# 注意：LSI第一个成分通常与测序深度相关，需要排除
DepthCor(seurat_obj)  # 检查各成分与深度的相关性
# 如果component_1和深度高度相关（r>0.7），后续分析排除它

# 4. 聚类
seurat_obj <- FindNeighbors(seurat_obj,
                             reduction = "lsi",
                             dims = 2:30)  # 排除第1个成分
seurat_obj <- FindClusters(seurat_obj, resolution = 0.5)

# 5. UMAP
seurat_obj <- RunUMAP(seurat_obj,
                       reduction = "lsi",
                       dims = 2:30)
DimPlot(seurat_obj, label = TRUE)

# Gene Activity——从染色质开放性推断基因表达
gene_activities <- GeneActivity(seurat_obj)  # 计算基因活性
seurat_obj[["RNA"]] <- CreateAssayObject(counts = gene_activities)
seurat_obj <- NormalizeData(seurat_obj, assay = "RNA")

# 可视化已知标记基因的Gene Activity
DefaultAssay(seurat_obj) <- "RNA"
FeaturePlot(seurat_obj,
            features = c("MS4A1", "CD3D", "CD14"),  # 标记基因
            ncol = 3)

# Motif分析——在开放区域中搜索TF结合基序
library(JASPAR2020)
library(TFBSTools)
library(motifmatchr)

# 获取motif数据库
pfm <- getMatrixSet(
  x = JASPAR2020,
  opts = list(collection = "CORE",
              tax_group = "vertebrates",
              all_versions = FALSE)
)

# 添加motif信息
seurat_obj <- AddMotifs(seurat_obj,
                         genome = BSgenome.Hsapiens.UCSC.hg38,
                         pfm = pfm)

# 差异可及性分析 → Motif富集
da_peaks <- FindMarkers(seurat_obj, ident.1 = "cluster_1",
                         only.pos = TRUE, min.pct = 0.1)

# 在差异peak中富集的motif
enriched_motifs <- FindMotifs(seurat_obj,
                               features = rownames(da_peaks))
head(enriched_motifs, 10)  # 查看Top10富集motif

# 安装ArchR
devtools::install_github("GreenleafLab/ArchR")  # 安装
library(ArchR)                                    # 加载

# ArchR使用Arrow文件（HDF5格式）处理大数据
# 直接从BAM/fragments创建Arrow文件
ArrowFiles <- createArrowFiles(
  inputFiles = "fragments.tsv.gz",    # 10x片段文件
  sampleNames = "PBMC",              # 样本名
  minTSS = 4,                        # 最小TSS富集
  minFrags = 1000,                   # 最小片段数
  addTileMat = TRUE,                 # 添加tile矩阵（500bp窗口）
  addGeneScoreMat = TRUE             # 添加基因评分矩阵
)

# 创建ArchR项目
proj <- ArchRProject(
  ArrowFiles = ArrowFiles,
  outputDirectory = "ArchR_output",
  copyArrows = TRUE
)

# 质控
proj <- filterDoublets(proj)  # 去除双细胞

# LSI降维
proj <- addIterativeLSI(
  ArchRProj = proj,
  useMatrix = "TileMatrix",    # 使用tile矩阵
  name = "IterativeLSI",
  iterations = 2,              # 迭代次数
  clusterParams = list(resolution = 0.2)
)

# 聚类
proj <- addClusters(
  input = proj,
  reducedDims = "IterativeLSI",
  resolution = 0.8
)

# UMAP
proj <- addUMAP(proj, reducedDims = "IterativeLSI")

# Peak calling（ArchR使用MACS2）
proj <- addGroupCoverages(proj, groupBy = "Clusters")
proj <- addReproduciblePeakSet(proj, groupBy = "Clusters")
proj <- addPeakMatrix(proj)

# TF Motif分析
proj <- addMotifAnnotations(proj, motifSet = "cisbp")
proj <- addBgdPeaks(proj)  # 添加背景peak
enrichMotifs <- peakAnnoEnrichment(proj, peakAnnotation = "Motif")

# TF足迹分析（footprinting）
proj <- addFootprints(proj, positions = motifPositions)
plotFootprints(proj, features = "CTCF", group = "Clusters")

# Bulk ATAC-seq流程
去接头(trim_galore --nextera) → 比对(Bowtie2)
→ 去重(Picard) → 去线粒体 → Tn5偏移校正
→ MACS2 peak calling → 质控(TSS富集/片段分布)
→ 差异peak → Motif富集 → 基因注释

# scATAC-seq Signac流程
Read10X_h5() → CreateChromatinAssay() → QC过滤
→ RunTFIDF() → FindTopFeatures() → RunSVD()
→ FindNeighbors(dims=2:30) → FindClusters()
→ RunUMAP() → GeneActivity() → AddMotifs()

# scATAC-seq ArchR流程
createArrowFiles() → ArchRProject()
→ filterDoublets() → addIterativeLSI()
→ addClusters() → addUMAP()
→ addReproduciblePeakSet() → addMotifAnnotations()

# 工具对比
ArchR: 速度快，内存低，百万级细胞，功能全面
Signac: Seurat生态，与scRNA-seq整合方便
选择: 大数据集→ArchR; Seurat用户→Signac

# 质控标准
片段数: >3000/细胞（scATAC）
TSS富集: >2（越高越好）
核小体信号: <4
线粒体比例: <5%
NFR比例: >50%

# LSI vs PCA
scATAC用LSI（TF-IDF + SVD），不用PCA
原因：ATAC-seq数据是二值/稀疏的，不适合PCA
注意：LSI第1个成分通常与深度相关，要排除