从scRNA-seq和scATAC-seq与区室图直接推断单个细胞中的高阶染色质结构
本·约翰逊
2022年4月26日
区域图谱:从scRNA-seq和scATAC-seq中直接推断单细胞高阶染色质
Compartmap扩展了Fortin和Hansen 2015年提出的方法，Genome Biology (https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0741-y)进行高阶染色质的直接推断单个细胞scRNA-seq和scATAC-seq。最初，Fortin和Hansen证明染色质构象可以从(sc)ATAC-seq、亚硫酸氢盐测序、DNase-seq和甲基化阵列推断出来，类似于HiC在组水平上提供的结果。因此，除了上述分析所提供的基本信息外，染色质状态也可以推断出来。
在这里，我们提出了一种从scRNA-seq和scATAC-seq推断群体和单细胞水平高阶染色质状态的方法。为了实现这一目标，我们使用来自scRNA-seq和scATAC-seq的染色体或全基因组信息，采用James-Stein估计器(JSE)对全局或目标平均值进行估计。此外，由于单细胞数据的稀疏性，我们采用自举程序来量化与推断的隔间的状态和边界相关的不确定性。然后，与正交分析类型相比，分隔图的输出可以直接可视化，并且/或嵌入UMAP或t-SNE之类的东西。此外，为了探索从分隔图推断出的高阶相互作用域，我们使用随机矩阵理论(RMT)方法来解析“格状”图案，类似于在Hi-C和scHi-C中观察到的图案。
快速开始示例数据
输入
隔间地图的预期输入是aRangedSummarizedExperiment对象。这些可以使用内置函数来构建importBigWig ()如果从BigWigs(推荐用于scRNA-seq)开始，或者从特性级对象(如a)开始SingleCellExperiment与rowRanges插槽中填充了每个特性的grange(参见下面的示例)。
在组和单细胞水平上推断高阶染色质结构域
####组级推断#####处理示例K562 scRNA-seq数据的chr14，并在1Mb分辨率下推断高阶染色质k562_compartments < -scCompartments(k562_scrna_chr14空空的=“chr14”，res =1 e6，组=真正的，引导=假，基因组=“hg19”，分析=“rna”）####单细胞级推断####要在单细胞中推断高阶域，并通过自举过程量化符号一致性，您可以运行:以10个单元格为例k562_scrna_chr14。子< -k562_scrna_chr14 (,样本（colnames(k562_scrna_chr14),大小=10，取代=假)]k562_compartments。启动< -scCompartments(k562_scrna_chr14.sub空空的=“chr14”，res =1 e6，组=假，引导=真正的，num.bootstraps =10，基因组=“hg19”，分析=“rna”）#翻转域符号，如果符号一致性在80%的自举中不一致k562_compartments.boot。修复< -fixCompartments(k562_compartments.bootmin.conf =0.8）查看GRangesList对象中的第一个单元格k562_compartments.boot.fix [[1]]GRanges对象有89个范围和13个元数据列:## seqnames range strand | PC## | ## chr14 1900000 -19999999 * | 0.824163## chr14 20000000-20999999 * | 0.427980chr14 21000000-21999999 * | 0.272336## chr14 2200000 -22999999 * | 0.017647## chr14 2300000 -23999999 * | -0.194594## ... ... ... ... . ...chr14 10300000 -103999999 * | -0.0672494chr14 10400000 -104999999 * | 0.0600952chr14 10500000 -105999999 * | 0.3147846## chr14 10600000 -106999999 * | 0.5694740## chr14 10700000 -107349539 * | 0.8241634##隔间评分引导。打开boot.closed## <字符> <数字> <数字> <数字>开0.824163 82## chr14:20000000-20999999 open 0.427980 8 2## chr14:21000000-21999999 open 0.272336 8 2## chr14:2200000 -22999999 open 0.017647 4 6## chr14:2300000 -23999999关闭-0.194594 2 8## ... ... ... ... ...## chr14:10300000 -103999999关闭-0.0672494 7 3## chr14:10400000 -104999999 open 0.0600952 7 3## chr14:10500000 -105999999 open 0.3147846 8 2## chr14:10600000 -106999999 open 0.5694740 10 0## chr14:1070000 -107349539 open 0.8241634 10 0## conf.est conf.est. upperci conf.est. lowerci## <数字> <数字> <数字>## chr14:190000 -19999999 0.716740 0.954113 0.479368## chr14:20000000-20999999 0.716740 0.954113 0.479368## chr14:21000000-21999999 0.716740 0.954113 0.479368## chr14:2200000 -22999999 0.427753 0.688396 0.167111## chr14:2300000 -23999999 0.716740 0.954113 0.479368## ... ... ... ...## chr14:10300000 -103999999 0.355507 0.607675 0.103338## chr14:10400000 -104999999 0.644493 0.896662 0.392325## chr14:10500000 -105999999 0.716740 0.954113 0.479368## chr14:10600000 -106999999 0.861234 1.000000 0.679113## chr14:1070000 -107349539 0.861234 1.000000 0.679113##翻转。分数flip.conf.est## <逻辑> <数字> <数字>## chr14:190000 -19999999 FALSE 0.824163 0.716740## chr14:20000000-20999999 FALSE 0.427980 0.716740## chr14:21000000-21999999 FALSE 0.272336 0.716740## chr14:2200000 -22999999 FALSE 0.017647 0.427753## chr14:230000 -23999999 FALSE -0.194594 0.716740## ... ... ... ...## chr14:10300000 -103999999 FALSE -0.0672494 0.355507## chr14:10400000 -104999999 FALSE 0.0600952 0.644493## chr14:10500000 -105999999 FALSE 0.3147846 0.716740## chr14:10600000 -106999999 FALSE 0.5694740 0.861234## chr14:1070000 -107349539 FALSE 0.8241634 0.861234## flip.conf.est.upperCI flip.conf.est.lowerCI## <数字> <数字>## chr14:190000 -19999999 0.954113 0.479368## chr14:20000000-20999999 0.954113 0.479368## chr14:21000000-21999999 0.954113 0.479368## chr14:2200000 -22999999 0.688396 0.167111## chr14:2300000 -23999999 0.954113 0.479368## ... ... ...## chr14:10300000 -103999999 0.607675 0.103338## chr14:10400000 -104999999 0.896662 0.392325## chr14:10500000 -105999999 0.954113 0.479368## chr14:10600000 -106999999 1.000000 0.679113## chr14:1070000 -107349539 1.000000 0.679113## -------## seqinfo: 1个来自未指定基因组的序列;没有seqlengths可见推断的染色质域
在组级或单单元级处理数据后，可以使用plotAB功能在隔间地图。值得注意的是，我们可以包括置信区间和中位数，从符号一致性的自举过程中得到的染色体范围的置信估计。在50%的情况下，这表明在开放状态和封闭状态之间的估计是平均的。这可能是由于数据的稀疏性或异构性。解决这个问题的一种可能的方法是，如果初始设置较低(例如10)，则增加执行的引导次数。或者，它可能是一个值得为您的数据集研究的区域。
#用plotAB从上面的单元格1中绘制“固定”结果#包括置信区间和中位数置信估计plotAB(k562_compartments.boot.fix [[1]],空空的=“chr14”，什么=“flip.score”，与。ci =真正的，median.conf =真正的）
结构域屈折变化的提取
当染色质结构域从“开放”状态转变为“关闭”状态时，提取它们的拐点以寻找附近的CTCF位点等通常是有兴趣的。我们可以使用getDomainInflections功能在隔间地图。
#提取单细胞结构域屈折#域变化可以用来寻找附近的CTCF站点等。k562_cell_1_inflections < -getDomainInflections(k562_compartments.boot.fix [[1]],什么=“flip.score”，res =1 e6，空空的=“chr14”，基因组=“hg19”）#显示拐点k562_cell_1_inflections## GRanges对象有12个范围和0个元数据列:## seqnames表示字符串## ## [1] chr14 19000000 *## [2] chr14 22999999 *## [3] chr14 27000000 *## [4] chr14 33999999 *## [5] chr14 42000000 *## ... ... ... ...## [8] chr14 78999999 *## [9] chr14 90000000 *## [10] chr14 96999999 *## [11] chr14 104000000 *## [12] chr14 107349539 *## -------## seqinfo: 1个来自未指定基因组的序列;没有seqlengths导入bigWigs作为输入隔间地图
目前推荐的scRNA-seq的compartmap输入文件是单单元格bigWigs，尽管它可以使用基于特征/计数的对象(下一节将演示)。单细胞bigWigs可以通过几种工具生成，例如deeptools (https://deeptools.readthedocs.io/en/latest/)．要导入bigWigs，我们可以使用importBigWig功能在隔间地图。这将读入一个bigWig文件，并有选择地总结为任意大小的bin。在compartmap手稿中使用的箱子大小是1kb，这也是我们在这里所做的。
我们可以导入一个bigWig文件列表，并将它们合并到一个rangedsummarizeexperiment对象中#列举大人物的例子要人< -list.files（执行（“extdata”，包=“compartmap”)，full.names =真正的）#生成1kb的箱子数据（“hg19.gr”，包=“compartmap”）kb_bins < -tileGenome（seqlengths =seqlengths(hg19.gr) [“chr14”]，tilewidth =1000，cut.last.tile.in.chrom =真正的）#进口bigwigs_lst < -拉普兰人(要人,函数(x) {importBigWig(x,垃圾箱=kb_bins,总结=真正的，基因组=“hg19”）})#结合bigwig_rse < -do.call(cbind bigwigs_lst)colnames(bigwig_rse) < -c（“cell_1”，“cell_2”）从一个基于特性或计数的对象开始
在没有或不能从bigWigs对象(例如scATAC)开始的情况下，我们可以从特征级或计数对象开始(例如scATAC)。SingleCellExperiment)．必须要做的两件事是确保rowRanges而且colnames为每个特征和单元格/样本设置。我们将使用K562中的scATAC-seq作为对象外观的示例，但也将展示添加的一种方法rowRanges到一个SingleCellExperiment，这对a也是一样的SummarizedExperiment对象。
#加载scATAC-seq示例#这些数据是使用csaw包进行预处理的数据（“k562_scatac_chr14”，包=“compartmap”）#显示整体对象结构＃请注意colname也在这里化验名称是countsk562_scatac_chr14类:rangedsummarizeexperiment## dim: 11404 279##元数据(6):spacing width…param final.ext## assays(1):计数## rownames: NULL## rowData名称(0):## colnames(279): cell_1 cell_2…cell_287 cell_288## colData名称(4):bam。文件总数显示rowRanges槽位为农庄对于每个特性
rowRanges(k562_scatac_chr14)包含11404个范围和0个元数据列的GRanges对象:## seqnames表示字符串## ## [1] chr14 19614351-19614500 *## [2] chr14 19614401-19614550 *## [3] chr14 19614451-19614600 *## [4] chr14 19614501-19614650 *## [5] chr14 19614551-19614700 *## ... ... ... ...## [11400] chr14 106938701-106938850 *## [11401] chr14 106938751-106938900 *## [11402] chr14 106938801-106938950 *## [11403] chr14 107280901-107281050 *## [11404] chr14 107280951-107281100 *## -------## seqinfo:来自未知基因组的1个序列但如果我们没有rowRanges比如SingleCellExperiment我们的工作，我们需要产生它们。因此，我们将展示一个如何添加的示例rowRanges从GTF文件到SingleCellExperiment．
#定义一个辅助函数来添加rowrange#修改自https://github.com/trichelab/velocessor/blob/master/R/import_plate_txis.R＃请注意你可以修改rtracklayer::import而不是子集到基因级别#如果使用特征级信息(例如scATAC的片段床文件)getRowRanges < -函数(gtf, asys) {gx < -子集(rtracklayer：：进口(gtf)类型= =“基因”）的名字(gx) < -gx＄gene_id农庄(gx) [rownames(容易)｝从SMART-seq3论文中导入HEK293T singlecel实验示例数据（“ss3_umi_sce”，包=“compartmap”）使用helper函数导入示例GTFgtf_rowranges < -getRowRanges（执行（“extdata / grch38_91_hsapiens.gtf.gz”，包=“compartmap”)，ss3_umi_sce)#添加row range到singlecelexperiment中rowRanges(ss3_umi_sce) < -gtf_rowranges#我们现在可以继续下一步并运行分隔图工作流运行分隔图工作流
一旦我们有了数据RangedSummarizedExperiment或者其他类型的SummarizedExperiment与rowRanges槽填满后，我们可以继续进行隔间图工作流程。我们将使用14号染色体手稿中显示的相同K562 scRNA-seq数据作为示例。
加载使用importBigWig导入的K562数据示例数据（“k562_scrna_raw”，包=“compartmap”）# TF-IDF变换信号k562_scrna_chr14_tfidf < -transformTFIDF（分析(k562_scrna_se_chr14))#将TF-IDF计数加回计数槽中的对象分析(k562_scrna_se_chr14“计数”) < -t(k562_scrna_chr14_tfidf)#计算组级染色质结构域k562_scrna_chr14_raw_domains < -scCompartments(k562_scrna_se_chr14空空的=“chr14”，res =1 e6，组=真正的，引导=真正的，num.bootstraps =10，基因组=“hg19”，分析=“rna”）#对于单单元格，运行以下命令# k562_scrna_chr14_raw_domains <- scCompartments(k562_scrna_se_chr14，# CHR = "chr14"，# res = 1e6，# group = FALSE，# bootstrap = TRUE，# num.bootstrap = 10，# genome = "hg19"，#化验= "rna")#“修复”不协调的标志一致性隔间k562_scrna_chr14_raw_domains。修复< -fixCompartments(k562_scrna_chr14_raw_domains)#剧情结果plotAB(k562_scrna_chr14_raw_domains.fix空空的=“chr14”，什么=“flip.score”，与。ci =真正的，median.conf =真正的）
高阶染色质相互作用图
我们可以从scRNA和scATAC中得出的另一个有趣的方面是通过使用随机矩阵理论方法对相关矩阵进行去噪来获得高阶相互作用域。这通常用Hi-C和scHi-C方法中显示的“格状”模式表示，其中较强的相关性(例如较大的红色强度)表示相对于较低的相关性的相互作用域。我们可以在这里做类似的事情。
使用随机矩阵理论的迭代去噪方法k562_scrna_chr14_rmt < -getDenoisedCorMatrix(k562_scrna_chr14res =1 e6，空空的=“chr14”，基因组=“hg19”，分析=“rna”，iter =2）##收缩容器与JSE。##计算相关性…# #做……##使用RMT去噪相关矩阵。迭代去噪。迭代:2#绘制交互地图plotCorMatrix(k562_scrna_chr14_rmtuppertri =真正的）