参考地址：
https://stuartlab.org/signac/articles/pbmc_vignette.html
http://events.jianshu.io/p/9f5e0d79ff67
数据下载：
https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5
https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_singlecell.csv
https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz
https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz.tbi

1. Signac的安装

Signac包可以直接使用install.packages("Signac")进行安装。但是我在安装时出现报错，显示缺少Rsamtools包和BiocGenerics包。这两个包是在Bioconductor上可以搜到的。可以使用BiocManager::install("Rsamtools") 和BiocManager::install("BiocGenerics")进行下载。

提示：此处若出现上述的包无法下载的情况，可以直接在对应网站上复制source的地址进行下载，或者将按照包下载到本地进行安装。我在下载的时候BiocGenerics是出现了这种情况，可能是我的bioconda的出现了某些问题。代码如下：

install.packages("https://www.bioconductor.org/packages/release/bioc/src/contrib/GenomeInfoDb_1.34.9.tar.gz",repos = NULL,type = "source")

2.加载流程中所有的安装包

library(EnsDb.Hsapiens.v75)  ##这个安装包可以用于人类的基因组注释文件
library(patchwork)
library(Signac)
library(Seurat)
library(ggplot2)
set.seed(1234)

3.分析前准备工作

Signac需要两个输入文件，这两个输入文件都可以由CellRanger获得
（1）Peak/Cell matrix：

• 横坐标为每个peak区域，纵坐标为每个细胞。类似于单细胞的gene expression count matrix 。
• 每个peak区域代表预测的每一个染色质开放区域。
• 矩阵内的数值代表该细胞在这个个peak位置中Tn5酶结合的个数
  （2）Fragment file：
• 该文件为一个列表，其中包含了所有单细胞中的所有的unique—fragement
• 这个文件很大，通常储存在磁盘（disk）上，而不需要导入内存（memory）中。
• 这个文件的作用就是展示一个细胞中所有的fragment片段，而不是仅包含peak里面的片段

4.导入数据并创建Signac和Seurat对象

#首先使用peak matrix 和meta.data创建Seurat对象，
#然后将储存fragment的路径导入seurat对象中
setwd("D:/R-data/learning/signac_data")  ##设置工作路径
getwd()

counts <- Read10X_h5(filename = "../signac_data/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")

metadata <- read.csv(
  file = "../signac_data/atac_v1_pbmc_10k_singlecell.csv",
  header = TRUE,
  row.names = 1
)

chrom_assay <- CreateChromatinAssay(
  counts = counts,
  sep = c(":", "-"),
  genome = 'hg19',
  fragments = '../signac_data/atac_v1_pbmc_10k_fragments.tsv.gz',
  min.cells = 10,
  min.features = 200
)
##ATAC-seq的数据被储存在一种特殊的assays（ChromatinAssay）
pbmc <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks",
  meta.data = metadata
)

pbmc

##通过打印这个peaks的assays，
#可以看见其中包含的其他的附加信息motif information、 gene annotations、genome information
pbmc[['peaks']]

#ChromatinAssay这种assays的结构讲解
#https://stuartlab.org/signac/articles/data_structures.html


#可以使用grangs函数应用于Seurat对象或者ChromatinAssay assays直接获取每个peaks在基因组上的ranges（范围）
granges(pbmc)
pbmc@assays$peaks@ranges


##我们可以向pbmc对象中添加人类基因组的gene的注释，添加后后面的一些函数（FUNCTION）就可以直接从seurat对象中提取gene的注释用于分析（如，画track）

#BiocManager::install("AnnotationHub")
##最上面人的基因组注释文件R包没有从biocondor上下载下来，我们从annotationhub上导入一样
library(AnnotationHub)
hub <- AnnotationHub() ##从网上下载基因组信息到hub中

query(hub, c("Homo sapiens","ensdb"))##检索ensdb数据中人类基因组
homo_ensdb <- hub[["AH53211"]]##提取某个注释文件

##设置注释
#library(EnsDb.Hsapiens.v75)#从该包中获取hg19版本人基因组注释
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)#将注释文件由EnsDb格式转换为GRanges格式
seqlevelsStyle(annotations) <- "UCSC"
Annotation(pbmc) <- annotations##给pbmc对象设置注释文件

5.计算QC标准（Computing QC Metrics）

计算QC标准（Computing QC Metrics）
不同研究的QC指标各不相同，据样本的生物系统、细胞活力和其他因素而有所不同。
（1）核小体条带分布模式 (Nucleosome banding pattern)：
DNA fragment的大小分布直方图应该能够展示包裹在一个核小体上的长度的强（核小体banding带型）。我们将每个细胞的这个值，并量化（单个核小体片段）和（无核小体片段）的比率（称之为核小体信号NS）。

（2）转录起始位点富集分数(Transcriptional start site (TSS) enrichment score) ：
ENCODE 项目根据将 TSS 区域中心的片段与 TSS 侧翼的片段的比率定义为 （ATAC-seq targating分数）（gene得分）
较差的ATAC-seq通常含有较低的TSS分数，TSSEnrichment() 函数可以计算每个细胞的TSS分数，并储存在meta.data中的TSS.enrichment一列。

（3）peaks中的fragments总数目(Total number of fragments in peaks)：
这个指标用于评估细胞测序深度（复杂性），reads数较少的细胞需要去除（测序深度低），reads数较高的细胞也需要去除（可能是由于双胞doublets, 细胞核团nuclei clumps, or other artefacts.）

（4）peaks中片段的比例 (Fraction of fragments in peaks)：
peaks中reads比例小于15~20%的细胞应当被去除，这些细胞通常代表低质量的细胞或其他技术错误。
！注意：这个指标可能对于我们所选的peaks较为敏感。

（5）比对到基因组"blacklist"区域中的reads比率 (Ratio reads in genomic blacklist regions)：
ENCODE 项目提供了基因组上的"blacklist"区域列表，该区域通常与人工信号相关的读取。reads比对到这些区域的比例较高的细胞（与reads比对到peaks区域比例相比）通常代表了技术误差，应将其删除。 Signac包中包含针对人类（hg19 和 GRCh38）、小鼠 (mm10)、果蝇 (dm3) 和秀丽隐杆线虫 (ce10)物种 在ENCODE “blacklist”区域。

（3）（4）（5）三个指标可以直接由CellRanger直接获取并存储在meta.data中。而一般使用signac为非10X的数据计算这三个指标

##计算QC指标
pbmc <- NucleosomeSignal(object = pbmc)##计算核小体条带分布模式

pbmc$nucleosome_group <- ifelse(pbmc$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
FragmentHistogram(object = pbmc, group.by = 'nucleosome_group')
#可以看到，单核小体/无核小体比率（基于上面的图）的异常值细胞具有不同的核小体条带模式。
#剩余的细胞呈现出成功ATAC seq实验的典型模式。

图片.png

pbmc <- TSSEnrichment(object = pbmc, fast = FALSE)##计算每个细胞的TSS富集分数
##当fast=T时，仅计算TSS分数，而不存储每个细胞Tn5插入频率的位置矩阵。？？
##这样后续将不能使用TSSPlot画图。
##*基于TSS分数将不同的细胞进行分组以及画所有TSS位点的可及性信号图来检查TSS富集分数
pbmc$high.tss <- ifelse(pbmc$TSS.enrichment > 2, 'High', 'Low')
TSSPlot(pbmc, group.by = 'high.tss') + NoLegend()

到TSS的距离.png

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 ##将blacklist区域中reads比例加进来
pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragments



##画上述5个指标的小提琴图
VlnPlot(
  object = pbmc,
  features = c('pct_reads_in_peaks', 'peak_region_fragments',
               'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal'),
  pt.size = 0.1,
  ncol = 5
)

图片.png

##去除不符合上述5个指标的细胞
pbmc <- subset(
  x = pbmc,
  subset = peak_region_fragments > 3000 &
    peak_region_fragments < 20000 &
    pct_reads_in_peaks > 15 &
    blacklist_ratio < 0.05 &
    nucleosome_signal < 4 &
    TSS.enrichment > 2
)
pbmc

6.标准化及线性降维：

（1）标准化/归一化（Normalization）：

• Signac是基于词频-逆文档频次 (term frequency-inverse document frequency，TF-IDF) 方法进行的归一化。它是将进行两步归一化操作，既跨细胞归一化以校正细胞测序深度的差异，又跨峰（peaks）归一化以为更罕见的峰提供更高的值。

（TF-IDF是一种用于信息检索与文本挖掘的常用加权技术。TF-IDF是一种统计方法，用以评估一字词对于一个文件集或一个语料库中的其中一份文件的重要程度。字词的重要性随着它在文件中出现的次数成正比增加，但同时会随着它在语料库中出现的频率成反比下降。）

（2）特征选择 (Feature selection)：

• 由于scATAC-seq的低变化性（因为测的是DNA，所有每个细胞在每个peak能测到的数目要么是1要么是0，最多可能测的在复制的细胞为4），使我们很难像scRNA-seq数据那样筛选高可变特征。
• 作为替代方案，我们仅选择前n%的peaks以进行降维。或者去除仅在少于n个细胞中检测到的peaks（使用FindTopFeatures()函数）
• 在这里我们一般选择降维的方法来加快运行速度（将min.cutoff 设置为‘q75’或者使用前25%的peaks ）（由于发现使用降维后的子集与全部的peaks得到的结果是相似的）
• 使用Seurat中的VariableFeatures() 可以自动用于降维来筛选的peaks。

（3）降维 (Dimension reduction)：

• 接下来，我们将对TD-IDF矩阵采用singular value decomposition (SVD)方法进行降维，并使用上一步筛选出的peaks（features），这将返回这个对象（object）的低维LSI代表（类似于scRNA-seq里面的PCA）

• 这种在TF-IDF后继续进行SVD的方法就是著名的latent semantic indexing (LSI)方法。

pbmc <- RunTFIDF(pbmc)
pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0')
pbmc <- RunSVD(pbmc)

##通常在进行完线性降维后，得到的第一个LSI所代表的通常为技术差异，而非生物学差异，我们一般将第一个LSI去除。
##可以使用DepthCor函数评估每个LSI成分与测序深度之间的相关性
DepthCor(pbmc)
# 从这个图上可以发现，第一个LSI与细胞的总计数counts 之间存在非常强的相关性，因此，我们将去除这个组分后再执行下游步骤

图片.png

7.非线性降维（ Non-linear dimension reduction and clustering）：

• 接下来，我们的细胞已经进入了一个低维空间，我们可以使用scRNA-seq分析的一些常规的方法（基于图像聚类以及非线性降维）来进行可视化。
• 使用函数RunUMAP(), FindNeighbors(), 和 FindClusters()

pbmc <- RunUMAP(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindNeighbors(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindClusters(object = pbmc, verbose = FALSE, algorithm = 3)
DimPlot(object = pbmc, label = TRUE) + NoLegend()

图片.png

8.创建基因活性表达矩阵（gene activity matrix）：

• 我们可以通过评估与基因相关的染色质区域的可及性来对基因组中的基因活性进行定量，并根据scATAC-seq数据创建一个新的gene活性表格。
• 我们在这里仅采用一种简单的方法：sum统计与基因体和启动子区域重叠的fragement的个数（Cicero也可实现这一目的）
• 为了创建gene活性表格，我们提取每个gene的坐标并延长至其上游的2kb区域，然后，我们统计每个细胞中映射到上述区域的fragement数目，这一步可以使用FeatureMatrix() 函数来进行，但是这步已经在 GeneActivity()函数中自动运行了。

#GeneActivity() 函数自动执行创建基因活性表达矩阵这一步骤
gene.activities <- GeneActivity(pbmc) #得到活性矩阵
#将活性矩阵添加到pbmc对象中，作为一个新的assays，并进行标准化。
pbmc <- NormalizeData(
  object = pbmc,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(pbmc$nCount_RNA)
)

• 接下来便可以使用经典的标记基因（canonical marker genes）的活性值可视化来帮助我们解释ATAC-seq簇。但是这个估计比scRNA-seq数据的估计的噪音值更大，这是因为ATAC-seq所代表的基因活性值与基因的表达不总是相对应的。
• 然而仅基于监督分析（supervised analysis）来对细胞亚群进行分类是具有挑战性的。

DefaultAssay(pbmc) <- 'RNA'

FeaturePlot(
  object = pbmc,
  features = c('MS4A1', 'CD3D', 'LEF1', 'NKG7', 'TREM1', 'LYZ'),
  pt.size = 0.1,
  max.cutoff = 'q95',
  ncol = 3
)

图片.png

9.与scRNA-seq数据进行整合（Integrating with scRNA-seq data）：

• 为了帮助解释ATAC-seq数据，我们可以基于来自于同一生物学系统的scRNA-seq的数据来对细胞进行分类。
• 我们使用一种标签转移（label transfer）和多组学的数据整合（cross-modality integration）的方法。目的是在基于基于活性矩阵的scATAC-seq数据和基于基因表达的scRNA-seq数据中鉴定到共同的相关模式。这个过程是对每个scRNA-seq数据鉴定出的细胞分类标签进行打分。
• 在这里，我们加载了已经经过预处理的人类PBMC的数据集，这数据集也是有10×提供的。你可以下载这个数据集以及处理该数据集生成seurat对象的代码。，同样，你也可以下载经过处理后的Seurat对象。

#加载下载后的PMBC的scRNA-seq的Seurat对象
pbmc_rna <- readRDS("../vignette_data/pbmc_10k_v3.rds")

transfer.anchors <- FindTransferAnchors(  #找到转移标签
  reference = pbmc_rna,
  query = pbmc,
  reduction = 'cca'
)

predicted.labels <- TransferData(          #预测的细胞标签
  anchorset = transfer.anchors,
  refdata = pbmc_rna$celltype,
  weight.reduction = pbmc[['lsi']],
  dims = 2:30
)

pbmc <- AddMetaData(object = pbmc, metadata = predicted.labels)  #将预测的细胞标签加入pmbc的meta.data中


##绘图对两种数据进行可视化
plot1 <- DimPlot(
  object = pbmc_rna,
  group.by = 'celltype',
  label = TRUE,
  repel = TRUE) + NoLegend() + ggtitle('scRNA-seq')

plot2 <- DimPlot(
  object = pbmc,
  group.by = 'predicted.id',
  label = TRUE,
  repel = TRUE) + NoLegend() + ggtitle('scATAC-seq')

plot1 + plot2

图片.png

• 根据图可以看到scRNA-seq的分类与我们仅基于scATAC-seq的数据几乎一致，我们可以很方便的注释scATAC-seq的细胞簇
• 但在图中可以看见细胞簇14被标注为CD4记忆细胞，并且这个小细胞簇具有较低的QC值，它可能代表低质量细胞。后续将该细胞类型去除。

pbmc <- subset(pbmc, idents = 14, invert = TRUE)
pbmc <- RenameIdents(
  object = pbmc,
  '0' = 'CD14 Mono',
  '1' = 'CD4 Memory',
  '2' = 'CD8 Effector',
  '3' = 'CD4 Naive',
  '4' = 'CD14 Mono',
  '5' = 'DN T',
  '6' = 'CD8 Naive',
  '7' = 'NK CD56Dim',
  '8' = 'pre-B',
  '9' = 'CD16 Mono',
  '10' = 'pro-B',
  '11' = 'DC',
  '12' = 'NK CD56bright',
  '13' = 'pDC'
)

10.鉴定不同细胞簇间差异可及性的peaks（Find differentially accessible peaks between clusters）：

• 我们采用differential accessibility (DA) test来鉴定不同细胞簇中可及性不同的区域。
• 我们对differential accessibility (DA) 利用逻辑回归的方法，并考虑fragment的数目的影响，来减少不同的测序深度对于结果的影响。
• 对于scATAC-seq这样的稀疏数据，我们发现通常需要将FindMarkers()函数中的min.pct从默认值（0.1）降低，因为默认值是专门为scRNA-seq数据设置的。
• 在这里，我们仅对Naive CD4 cells 和 CD14 monocytes进行比较。其他的所有细胞类型都可以进行这样的比较，并且可以使用Seurat中提供的函数进行可视化。（小提琴图、feature 图、dot点图、热图等等）

DefaultAssay(pbmc) <- 'peaks'  #将默认的assays转换回peaks

da_peaks <- FindMarkers(
  object = pbmc,
  ident.1 = "CD4 Naive",
  ident.2 = "CD14 Mono",
  test.use = 'LR',
  latent.vars = 'peak_region_fragments'
)

head(da_peaks)##查看差异peaks


plot1 <- VlnPlot(
  object = pbmc,
  features = rownames(da_peaks)[1],#对第一显著的差异peaks画小提琴图
  pt.size = 0.1,
  idents = c("CD4 Naive","CD14 Mono")
)
plot2 <- FeaturePlot(
  object = pbmc,
  features = rownames(da_peaks)[1],#对第一显著的peaks画feature图
  pt.size = 0.1
)

plot1 | plot2

图片.png

• 除了使用DA test的方法，我们还可以通过比较两组细胞的Fold Change accessibility 来鉴定DA区域。
• 这个方法可能比运行DA test的方法快很多，但这个方法没有做统计检验，所以无法避免一些潜在因素的影响，如细胞间的测序深度差异造成的影响。
• 但使用这个方法可以快速对数据进行探索，FoldChange()函数可以进行上述操作。

fc <- FoldChange(pbmc, ident.1 = "CD4 Naive", ident.2 = "CD14 Mono")
head(fc)

• 单纯的peak位置很难单独进行解释，我们使用ClosestFeature()函数查找这些peaks的最近的基因。如果我们查看基因列表，会发现naive T细胞中的开放峰与BCL11B和GATA3（T细胞分化的关键调控因子）等基因接近，而单核细胞中开放的峰与CEBPB（单核细胞分化的关键调控因子）基因接近。
• 在发现这些gene之后，还可以进行gene的富集分析，可以使用 GOstats R包进行富集分析。

open_cd4naive <- rownames(da_peaks[da_peaks$avg_log2FC > 3, ])
open_cd14mono <- rownames(da_peaks[da_peaks$avg_log2FC < -3, ])

closest_genes_cd4naive <- ClosestFeature(pbmc, regions = open_cd4naive)
closest_genes_cd14mono <- ClosestFeature(pbmc, regions = open_cd14mono)

head(closest_genes_cd4naive)
head(closest_genes_cd14mono)