1. 切割位点

ATACseq 应该在较小的保护区（如转录因子结合位点）周围生成较短的片段（我们的无核小体区域）。

因此，我们可以在不同组织/细胞类型/样本中寻找围绕感兴趣基序的切割位点堆积。

为了从我们的 BAM 文件中生成切割位点，我们首先将读取大小调整为 1bp，并根据链进行 4/-5 bp 的偏移，以调整插入 Tn5 转座酶的预期偏移。

在这里，我们将识别通过任意截止的 CTCF 基序，然后使用 soGGi 在它们周围绘制切割位点。我们将跳回我们的 Greenleaf 数据集来执行此操作。

2. 查找 motifs

我们需要确定 CTCF 基序在基因组中的位置，因此首先我们需要知道 CTCF 基序是什么样的。

motifDB 包包含来自公共数据库（例如 JASPAR）的有关 Motif 的信息。在这里，我们使用带有我们感兴趣的主题 (CTCF) 的 query() 函数来提取 CTCF 主题。

library(MotifDb)
library(Biostrings)
library(BSgenome.Hsapiens.UCSC.hg19)
CTCF <- query(MotifDb, c("CTCF"))
CTCF

我们可以为 CTCF 提取一个点权重矩阵，它指定了 DNA 碱基出现在 CTCF 基序中的可能性。在这里，我们从 Human JASPAR Core 数据库中提取 CTCF 的主题。

names(CTCF)

ctcfMotif <- CTCF[[1]]
ctcfMotif[, 1:4]

3. PWMs 可视化

我们可以使用 seqLogo 包和 seqLogo 函数可视化主题中 DNA 碱基的频率。

library(seqLogo)
seqLogo(ctcfMotif)

4. PWMs 搜索

我们现在可以将 matchPWM() 函数与我们新获得的 CTCF PWM 一起使用。在这里，我们将使用 BSgenome 库中为人类 BSgenome.Hsapiens.UCSC.hg19 提供的序列搜索 Chr20 上的序列。结果是一个 Views 对象，类似于 IRanges 对象。

myRes <- matchPWM(ctcfMotif, BSgenome.Hsapiens.UCSC.hg19[["chr20"]])
myRes

我们需要将 Views 对象转换为 GRanges，以便我们可以在 soGGi 中使用它们来绘制切割站点。

toCompare <- GRanges("chr20", ranges(myRes))
toCompare

5. 切割位点分析

要绘制切割位点，我们希望只考虑读取的 5' 端，并且需要调整已知的 5' 读取偏移量到实际 T5 切割位点。

这将涉及捕获读数的 5' 端并将正链和负链上的读数分别移动 4bp 或 -5bp。

首先，我们读入我们的无核小体区域 BAM 文件并提取读取对。

BAM <- "~/Downloads/ATAC_Workshop/ATAC_Data/ATAC_BAM/Sorted_ATAC_50K_2_openRegions.bam"
atacReads_Open <- readGAlignmentPairs(BAM)
read1 <- first(atacReads_Open)
read2 <- second(atacReads_Open)
read2[1, ]

现在我们可以根据链将两个读取对的 5' 端移动 4bp 或 -5bp。这从两个读数中产生了我们所有切割位点的 GRanges。

Firsts <- resize(granges(read1), fix = "start", 1)
First_Pos_toCut <- shift(granges(Firsts[strand(read1) == "+"]), 4)
First_Neg_toCut <- shift(granges(Firsts[strand(read1) == "-"]), -5)

Seconds <- resize(granges(read2), fix = "start", 1)
Second_Pos_toCut <- shift(granges(Seconds[strand(read2) == "+"]), 4)
Second_Neg_toCut <- shift(granges(Seconds[strand(read2) == "-"]), -5)

test_toCut <- c(First_Pos_toCut, First_Neg_toCut, Second_Pos_toCut, Second_Neg_toCut)
test_toCut[1:2, ]

现在我们可以使用 coverage() 函数使用切割位点位置的 GRanges 生成整个基因组切割位点的 RLElist。

cutsCoverage <- coverage(test_toCut)
cutsCoverage20 <- cutsCoverage["chr20"]
cutsCoverage20[[1]]

我们可以使用带有 soGGi 的 RLElist 围绕我们发现的 CTCF 图案生成切割位点图。

我们将格式更改为 rlelist，并将 distanceAround 参数更改为 500bp。

CTCF_Cuts_open <- regionPlot(cutsCoverage20, testRanges = toCompare, style = "point",
    format = "rlelist", distanceAround = 500)

现在我们可以使用 plotRegion() 函数绘制切割点。

plotRegion(CTCF_Cuts_open, outliers = 0.001) + ggtitle("NucFree Cuts Centred on CTCF") +
    theme_bw()