单细胞｜线粒体基因型和DNA可及性联合分析

引言

本研究分析了两组包含DNA可及性和线粒体突变信息的数据集，数据来自结直肠癌（CRC）患者的样本。

为了展示如何联合分析线粒体DNA变异和染色质可及性，将通过一个实例来说明，该实例涉及对来自原发性结直肠腺癌的细胞进行分析。这些细胞样本中包含了恶性上皮细胞和肿瘤内部浸润的免疫细胞。

导入 DNA 可及性数据

首先导入单细胞染色质可及性测序（scATAC-seq）数据，并按照scATAC-seq数据的标准分析流程，构建一个Seurat对象。

library(Signac)
library(Seurat)
library(ggplot2)
library(patchwork)
library(EnsDb.Hsapiens.v75)

# load counts and metadata from cellranger-atac
counts <- Read10X_h5(filename = "CRC_v12-mtMask_mgatk.filtered_peak_bc_matrix.h5")
metadata <- read.csv(
  file = "CRC_v12-mtMask_mgatk.singlecell.csv",
  header = TRUE,
  row.names = 1
)

# load gene annotations from Ensembl
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)

# change to UCSC style since the data was mapped to hg19
seqlevels(annotations) <- paste0('chr', seqlevels(annotations))
genome(annotations) <- "hg19"

# create object
crc_assay <- CreateChromatinAssay(
  counts = counts,
  sep = c(":", "-"),
  annotation = annotations,
  min.cells = 10,
  fragments = 'CRC_v12-mtMask_mgatk.fragments.tsv.gz'
)
crc <- CreateSeuratObject(
  counts = crc_assay,
  assay = 'peaks',
  meta.data = metadata
)

crc[["peaks"]]

## ChromatinAssay data with 81787 features for 3535 cells
## Variable features: 0 
## Genome: 
## Annotation present: TRUE 
## Motifs present: FALSE 
## Fragment files: 1

质控流程

能够对单细胞ATAC-seq（scATAC-seq）数据进行标准的质量评估，并依据这些评估指标剔除质量不佳的细胞。

# Augment QC metrics that were computed by cellranger-atac
crc$pct_reads_in_peaks <- crc$peak_region_fragments / crc$passed_filters * 100
crc$pct_reads_in_DNase <- crc$DNase_sensitive_region_fragments / crc$passed_filters * 100
crc$blacklist_ratio <- crc$blacklist_region_fragments / crc$peak_region_fragments

# compute TSS enrichment score and nucleosome banding pattern
crc <- TSSEnrichment(crc)
crc <- NucleosomeSignal(crc)

# visualize QC metrics for each cell
VlnPlot(crc, c("TSS.enrichment", "nCount_peaks", "nucleosome_signal", "pct_reads_in_peaks", "pct_reads_in_DNase", "blacklist_ratio"), pt.size = 0, ncol = 3)

# remove low-quality cells
crc <- subset(
  x = crc,
  subset = nCount_peaks > 1000 &
    nCount_peaks < 50000 &
    pct_reads_in_DNase > 40 &
    blacklist_ratio < 0.05 &
    TSS.enrichment > 3 & 
    nucleosome_signal < 4
)
crc

## An object of class Seurat 
## 81787 features across 1861 samples within 1 assay 
## Active assay: peaks (81787 features, 0 variable features)
##  2 layers present: counts, data

线粒体变异数据导入

接下来，将导入这些细胞的线粒体DNA变异数据，这些数据是通过mgatk工具进行量化的。Signac包中的ReadMGATK()函数使得mgatk的输出能够直接以适合于Signac后续分析的格式读入R语言环境中。在此过程中，将这些数据加载进来，并将其作为一个新检测添加到Seurat对象中。

# load mgatk output
mito.data <- ReadMGATK(dir = "crc/")

# create an assay
mito <- CreateAssayObject(counts = mito.data$counts)

# Subset to cell present in the scATAC-seq assat
mito <- subset(mito, cells = Cells(crc))

# add assay and metadata to the seurat object
crc[["mito"]] <- mito
crc <- AddMetaData(crc, metadata = mito.data$depth[Cells(mito), ], col.name = "mtDNA_depth")

可以查看每个细胞的线粒体测序深度，并根据线粒体测序深度进一步对细胞进行子集化。

VlnPlot(crc, "mtDNA_depth", pt.size = 0.1) + scale_y_log10()

# filter cells based on mitochondrial depth
crc <- subset(crc, mtDNA_depth >= 10)
crc

## An object of class Seurat 
## 214339 features across 1359 samples within 2 assays 
## Active assay: peaks (81787 features, 0 variable features)
##  2 layers present: counts, data
##  1 other assay present: mito

降维和聚类

接下来，将利用单细胞ATAC-seq（scATAC-seq）数据进行常规的降维和聚类分析，目的是识别不同的细胞群组。

crc <- RunTFIDF(crc)
crc <- FindTopFeatures(crc, min.cutoff = 10)
crc <- RunSVD(crc)
crc <- RunUMAP(crc, reduction = "lsi", dims = 2:50)
crc <- FindNeighbors(crc, reduction = "lsi", dims = 2:50)
crc <- FindClusters(crc, resolution = 0.7, algorithm = 3)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 1359
## Number of edges: 63301
## 
## Running smart local moving algorithm...
## Maximum modularity in 10 random starts: 0.7329
## Number of communities: 6
## Elapsed time: 0 seconds

DimPlot(crc, label = TRUE) + NoLegend()

基因得分计算

为了更好地理解这些细胞群的特征并确定它们的细胞类型，将计算基因的活性得分，方法是将基因本体区域和启动子区域的DNA可及性进行累加。

# compute gene accessibility
gene.activities <- GeneActivity(crc)

# add to the Seurat object as a new assay
crc[['RNA']] <- CreateAssayObject(counts = gene.activities)

crc <- NormalizeData(
  object = crc,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(crc$nCount_RNA)
)

展示特定基因活性得分

在CRC数据集中，识别了以下不同细胞类型的标志性基因：

• EPCAM：上皮细胞的标志
• TREM1：髓系细胞的标志
• PTPRC（也称为CD45）：泛免疫细胞的标志
• IL1RL1：嗜碱性粒细胞的标志
• GATA3：T细胞的标志

DefaultAssay(crc) <- 'RNA'

FeaturePlot(
  object = crc,
  features = c('TREM1', 'EPCAM', "PTPRC", "IL1RL1","GATA3", "KIT"),
  pt.size = 0.1,
  max.cutoff = 'q95',
  ncol = 2
)

使用这些基因评分值，可以分配簇标识：

crc <- RenameIdents(
  object = crc,
  '0' = 'Epithelial',
  '1' = 'Epithelial',
  '2' = 'Basophil',
  '3' = 'Myeloid_1',
  '4' = 'Myeloid_2',
  '5' = 'Tcell'
)

在髓系细胞群中，发现有一个亚群表现出较低的峰值区域片段比例，整体线粒体测序覆盖度也不高，并且其核小体的排列模式与其他亚群有所区别。

p1 <- FeatureScatter(crc, "mtDNA_depth", "pct_reads_in_peaks") + ggtitle("") + scale_x_log10()
p2 <- FeatureScatter(crc, "mtDNA_depth", "nucleosome_signal") + ggtitle("") + scale_x_log10()

p1 + p2 + plot_layout(guides = 'collect')

观察到，大多数FRIP值较低的细胞属于髓系1群。这些细胞很可能是肿瘤内的粒细胞，它们在染色质可及性方面的富集程度相对较低。与此相似，这种细胞类型的核染色质包装方式也较为特殊，导致其mtDNA覆盖度略低于髓系2群。

识别有意义的mtDNA变异

接下来，可以确定线粒体基因组中在不同细胞间存在差异的位点，并根据这些变异在细胞中的出现频率将细胞分为不同的克隆型。Signac在识别高质量变异时，采用了原始mtscATAC-seq研究中建立的原则。

variable.sites <- IdentifyVariants(crc, assay = "mito", refallele = mito.data$refallele)
VariantPlot(variants = variable.sites)

上图清晰地展示了一组变异，这些变异具有较高的变异等位基因比率（VMR）和链一致性。理论上，高链一致性减少了等位基因频率受测序错误影响的可能性（这种错误通常只发生在一个链上，而不是另一个链上，这是由于前导核苷酸的内容和mtDNA基因分型中的一个常见错误）。另一方面，具有高VMR的变异更可能是克隆变异，因为这些变异的替代等位基因倾向于在特定细胞中聚集，而不是均匀分布在所有细胞中，后者则可能表明其他一些技术问题。

注意到，那些具有极低VMR和极高链一致性的变异是这个样本的同源变异。虽然这些变异在某些情况下可能具有研究价值（例如，在供体样本的去复用过程中），但对于推断亚克隆结构，它们并没有太大的用处。

根据这些标准，可以筛选出一组在不同细胞中表现出差异的、具有信息价值的线粒体变异。

# Establish a filtered data frame of variants based on this processing
high.conf <- subset(
  variable.sites, subset = n_cells_conf_detected >= 5 &
    strand_correlation >= 0.65 &
    vmr > 0.01
)

high.conf[,c(1,2,5)]

##          position nucleotide      mean
## 1227G>A      1227        G>A 0.0090107
## 6081G>A      6081        G>A 0.0031485
## 9804G>A      9804        G>A 0.0035833
## 12889G>A    12889        G>A 0.0217851
## 9728C>T      9728        C>T 0.0150412
## 16147C>T    16147        C>T 0.6582835
## 824T>C        824        T>C 0.0050820
## 2285T>C      2285        T>C 0.0034535
## 16093T>C    16093        T>C 0.0093126

一些关键点值得注意。首先，在12个变异中，有10个在群体中的出现频率不到1%。但是，16147C>T的频率大约为62%。将会认识到这是一个标记上皮细胞的克隆变异。此外，所有检测到的变异都是转换类型（A - G 或 C - T），而不是颠换类型（A - T 或 C - G）。这与对线粒体基因组中这些突变产生方式的认知相吻合。

计算变异等位基因频率

目前拥有存储在mito检测中的每个链的信息，这使得可以评估链的一致性。现在，已经确定了一组高可信度的信息变异，可以使用AlleleFreq()函数，为这些变异创建一个新的检测，包含每个细胞的链折叠等位基因频率计数。

crc <- AlleleFreq(
  object = crc,
  variants = high.conf$variant,
  assay = "mito"
)
crc[["alleles"]]

## Assay data with 9 features for 1359 cells
## First 9 features:
##  1227G>A, 6081G>A, 9804G>A, 12889G>A, 9728C>T, 16147C>T, 824T>C,
## 2285T>C, 16093T>C

展示变异分布

随着等位基因频率被存储为额外的检测数据，可以利用Seurat的标准功能来展示这些频率如何在不同细胞间分布。在这里，通过FeaturePlot()和DoHeatmap()展示了变异的一个子集，以便直观地理解它们的分布情况。

DefaultAssay(crc) <- "alleles"
alleles.view <- c("12889G>A", "16147C>T", "9728C>T", "9804G>A")
FeaturePlot(
  object = crc,
  features = alleles.view,
  order = TRUE,
  cols = c("grey", "darkred"),
  ncol = 4
) & NoLegend()

DoHeatmap(crc, features = rownames(crc), slot = "data", disp.max = 1) +
  scale_fill_viridis_c()

在这些选定的变异中，可以观察到一些有趣的分布模式。变异16147C>T几乎存在于所有的上皮细胞中，而且基本上只出现在上皮细胞内（少数例外情况也对应着UMAP和聚类结果不一致的案例）。这个变异的等位基因频率达到了100%，这强烈表明肿瘤的原始细胞在发生突变后进行了克隆性扩增。接着，至少发现了3个变异——1227G>A、12889G>A和9728C>T，它们主要特定于定义亚克隆的上皮细胞中。另外，包括3244G>A、9804G>A和824T>C在内的其他变异，特别在免疫细胞群体中被发现，这暗示这些变异可能源自一个共同的造血祖细胞（很可能位于骨髓中）。