ArchR官网教程学习笔记14:ArchR的Footprinti

2020-12-02 本文已影响0人生信start_site

转录因子(TF)footprinting分析可以预测TF在特定位点的精确结合位置。这是因为直接与转录因子结合的DNA碱基实际上受到保护，不会发生转位，而与转录因子结合紧邻的DNA碱基是可以接近的。

理想情况下，TF footprinting是在单个位点进行，以确定TF的精确结合位置。然而，在实践中，这需要非常高的测序深度，通常比大多数用户从bulk或单细胞ATAC-seq获得的深度要高得多。为了解决这个问题，我们可以用多个预测TF结合实例结合Tn5插入位置。例如，我们可以取所有含有CTCF motif的峰，并在整个基因组中为CTCF生成一个TF足迹集合。

这个footprint的准确性依赖于为感兴趣的TF生成一个可靠的预测结合位点的列表。ArchR使用addMotifAnnotations()函数通过搜索任何匹配motif的DNA序列的峰区域来实现这一点。根据感兴趣的motif的简并性，这可能是充分的，也可能是不充分的。这些motif注释以每个峰为基础(0 =没有motif, 1 = motif present)的二进制表示形式添加到ArchRProject中。一旦你有了这些motif注释，ArchR使用getfootprint()函数执行footprinting分析，该函数接受一个ArchRProject对象和一个包含motifs位置的GenomicRanges对象作为输入。这些位置可以通过getPositions()函数从ArchRProject访问。然后可以使用plotFootprints()函数绘制这些footprints。

最重要的是，ArchR的 footprinting 分析考虑了已知的Tn5插入序列偏差。为了做到这一点，ArchR使用六聚体位置频率矩阵和Tn5插入位点k-mer频率矩阵:

所有这些放在一起，这个工作流程生成考虑Tn5插入偏差的footprint图。

ArchR支持motif足迹和用户自定义特征的footprinting，这两个在本章讨论。

（一）Motif Footprinting

重要的是，由本教程数据生成的footprints并不像预期的那样干净，但这是因为本教程数据集较小。从更大的数据集产生的footprints预计将显示更少的variation。

在做footprinting分析的时候，我们首先要做的是获得相关motif的位置。为此，我们调用getPositions()函数。这个函数有一个名为name的可选参数，它可以接受peakAnnotation对象的名称，我们从中获得位置。如果name = NULL，那么ArchR将使用peakAnnotation槽中的第一个条目。在下面的示例中，我们没有指定名称，ArchR使用第一个条目，即我们的CIS-BP motifs。

> motifPositions <- getPositions(projHeme5)

这时创建了一个GRangesList对象，每个TF motif由一个分割的GRanges对象表示。

> motifPositions
GRangesList object of length 870:
$TFAP2B_1
GRanges object with 16660 ranges and 1 metadata column:
          seqnames              ranges strand |            score
             <Rle>           <IRanges>  <Rle> |        <numeric>
      [1]     chr1       873916-873927      + | 8.31937625544643
      [2]     chr1       873916-873927      - | 8.31937625544643
      [3]     chr1       875505-875516      + |   10.21484637226
      [4]     chr1       875505-875516      - | 9.05303812217238
      [5]     chr1       896671-896682      + | 9.95732790593787
      ...      ...                 ...    ... .              ...
  [16656]     chrX 153991101-153991112      + | 8.38589666203426
  [16657]     chrX 154299568-154299579      + | 8.89537817432743
  [16658]     chrX 154664929-154664940      - | 8.16160181551093
  [16659]     chrX 154807684-154807695      + | 9.57032220106342
  [16660]     chrX 154807684-154807695      - | 10.6068302166026
  -------
  seqinfo: 23 sequences from an unspecified genome; no seqlengths

...
<869 more elements>

我们可以把这个GRangesList子集化为几个我们感兴趣的TF motifs。因为SREBF1 TF在我们搜索“EBF1”时出现，所以我们使用%ni%辅助函数从下游分析中删除它。

> motifs <- c("GATA1", "CEBPA", "EBF1", "IRF4", "TBX21", "PAX5")
> markerMotifs <- unlist(lapply(motifs, function(x) grep(x, names(motifPositions), value = TRUE)))
> markerMotifs <- markerMotifs[markerMotifs %ni% "SREBF1_22"]
> markerMotifs
[1] "GATA1_383" "CEBPA_155" "EBF1_67"   "IRF4_632"  "TBX21_780"
[6] "PAX5_709"

为了准确地分析TF footprints，需要进行大量的读取。因此，对细胞进行分组以创建伪批量ATAC-seq文件，然后可用于TF footprinting。这些伪批量文件以为group覆盖率文件储存，这是我们在前一章中创建的用来执行peak calling的。如果你还没有在ArchRProject中添加group coverages，现在可以添加:

> projHeme5 <- addGroupCoverages(ArchRProj = projHeme5, groupBy = "Clusters2")

通过计算group coverages，我们现在可以使用getfootprint()函数计算之前选择的标记motif子集。尽管ArchR实现了一个高度优化的footprinting工作流，但建议对motif的子集而不是所有的motif执行footprinting分析。因此，我们通过position参数向footprint提供motifs的子集:

> seFoot <- getFootprints(
  ArchRProj = projHeme5, 
  positions = motifPositions[markerMotifs], 
  groupBy = "Clusters2"
)

一旦我们获得了这些footprints，我们就可以使用plotfootprint()函数来绘制它们。此函数可以同时以各种方式标准化footprints。下一节将讨论这种标准化和footprints的实际绘制。

（二）Footprints对于Tn5偏差的标准化

使用ATAC-seq数据进行TF足迹分析的一个主要挑战是Tn5转座酶的插入序列偏差，这可能导致TF footprint的错误分类。为了考虑Tn5插入偏差，ArchR识别围绕每个Tn5插入位点的k-mer(用户定义长度，默认长度6)序列。为了进行分析，ArchR鉴定每个伪批量的单碱基分辨率Tn5插入位点，调整这些1-bp位点到k-bp 窗口 (- k / 2 + (k / 2 - 1) bp从插入)。然后使用Biostrings 包创建一个使用oligonucleotidefrequency(w = k, simplify.as ="collapse") 的k-mer频率表。然后，ArchR使用与BSgenome相关的基因组文件和相同的函数计算预期的全基因组范围的k-mers。为了计算伪批量footprint的插入偏差，ArchR创建了一个k-mer频率矩阵，它表示为从motif中心穿过一个窗口+/- N bp(用户定义，默认为250 bp)的所有可能的k-mers。然后，在每个motif位点上迭代，ArchR将定位的k-mers填充到k-mer频率矩阵中。然后在全基因组范围内计算每个motif的位置。使用样本的k-mer频率表，ArchR可以通过将k-mer位置频率表乘以观察/期望的Tn5 k-mer频率来计算预期的Tn5插入。

所有这些都是在plotfootprint()函数内部进行的。

（１）减去Tn5偏差

一种标准化方法是从footprinting信号中减去Tn5偏差。这个标准化是通过在调用plotFootprints()时设置normMethod = "Subtract"来执行的。

> plotFootprints(
  seFoot = seFoot,
  ArchRProj = projHeme5, 
  normMethod = "Subtract",
  plotName = "Footprints-Subtract-Bias",
  addDOC = FALSE,
  smoothWindow = 5
)

默认情况下，这些图将保存在ArchRProject的outputDirectory中。如果你要求绘制所有的motif并返回一个ggplot对象，这个ggplot对象将会非常大。下面是一个已经减去偏差的motif footprints分析结果：

（２）没有标准化Tn5偏差的Footprinting

虽然我们强烈建议对Tn5序列插入偏差的footprints进行规范化，但也可以通过在plotFootprints()函数中设置normMethod = "None"来执行不标准化的footprinting。

> plotFootprints(
  seFoot = seFoot,
  ArchRProj = projHeme5, 
  normMethod = "None",
  plotName = "Footprints-No-Normalization",
  addDOC = FALSE,
  smoothWindow = 5
)

以这种方法得到的结果会是这样：

这里比较奇怪的是，我的这一步并没有出图（所以上面的图用的是官网的图片），而输出信息都是和官网是一致的。我查阅了google和百度，用没有发现有和我有相同的问题的人。这一步等后面的学习结束后，还会再返回来研究一下的。

（３）特征footprinting分析

除了motif的footprinting外，ArchR还允许对任何用户自定义的特征集进行footprinting。为了演示这一功能，我们将使用plotFootprints()函数创建一个TSS插入profile文件(之前在数据质量控制一节中介绍过)。TSS插入文件只是footprinting的一个特殊子案例。

正如前一节所讨论的，使用从伪批量重复派生的group coverage文件来执行footprinting。我们在前一章中创建这些是为了执行peak calling。如果你还没有在ArchRProject中添加group coverages，现在可以添加：

> projHeme5 <- addGroupCoverages(ArchRProj = projHeme5, groupBy = "Clusters2")

我们在创建TSS插入配置文件时没有对Tn5偏差进行标准化。与我们之前的分析主要不同的是，我们指定了flank = 2000来扩展每个TSS的footprints 2000 bp。

> seTSS <- getFootprints(
  ArchRProj = projHeme5, 
  positions = GRangesList(TSS = getTSS(projHeme5)), 
  groupBy = "Clusters2",
  flank = 2000
)

绘制每个细胞群的TSS插入profiles，使用plotFootprints()：

> plotFootprints(
  seFoot = seTSS,
  ArchRProj = projHeme5, 
  normMethod = "None",
  plotName = "TSS-No-Normalization",
  addDOC = FALSE,
  flank = 2000,
  flankNorm = 100
)

同样的，这一步的图也是没有出来，但是奇怪的是并没有报错。后面有时间会来复盘一下这一章节。