ChIP-seq 分析：Consensus Peaks（14）

动动发财的小手，点个赞吧！

1. 数据读入

首要，咱们需要将来自 MACS2 的峰值调用读取到 R 中。

咱们将检查的 Myc peak 调用坐落 peaks 目录中，因此咱们在这里运用 dir() 函数列出与咱们预期的文件形式匹配的一切文件。

peakFiles <- dir("data/peaks/", pattern = "*.peaks", full.names = TRUE)
peakFiles

咱们能够循环遍历咱们的制表符分隔文件（伪装成 .xls 函数）并运用循环将它们作为 data.frames 列表导入到 R 中。

macsPeaks_DF <- list()
for (i in 1:length(peakFiles)) {
    macsPeaks_DF[[i]] <- read.delim(peakFiles[i], comment.char = "#")
}
length(macsPeaks_DF)

现在有了咱们的 data.frames 峰值调用列表，咱们循环遍历列表并为每个峰值调用创立一个 GRanges。

请记住，您也能够运用 rtracklayer 的导入功能来履行此操作。

library(GenomicRanges)
macsPeaks_GR <- list()
for (i in 1:length(macsPeaks_DF)) {
    peakDFtemp <- macsPeaks_DF[[i]]
    macsPeaks_GR[[i]] <- GRanges(seqnames = peakDFtemp[, "chr"], IRanges(peakDFtemp[,
        "start"], peakDFtemp[, "end"]))
}
macsPeaks_GR[[1]]

咱们将要为咱们的峰值呼叫分配一组合理的称号。

咱们能够将 gsub() 和 basename() 函数与咱们的文件名一起运用来创立一些样本称号。

basename() 函数承受文件途径（例如咱们的 bam 文件的途径）并仅回来文件名（删除目录途径）。

gsub() 函数承受要替换的文本、替换文本和要替换的字符向量。

fileNames <- basename(peakFiles)
fileNames

sampleNames <- gsub("_peaks.xls", "", fileNames)
sampleNames

现在咱们有一个作为 GRanges 目标的峰值调用的命名列表。

咱们能够运用 GRangesList() 函数将 GRanges 目标列表转换为 GRangesList。

macsPeaks_GRL <- GRangesList(macsPeaks_GR)
names(macsPeaks_GRL) <- sampleNames
class(macsPeaks_GRL)

names(macsPeaks_GRL)

2. GRangesList 目标

GRangesList 目标的行为与咱们的规范列表相同。在这里，咱们运用 lengths() 函数来获取每个重复中的峰数。

lengths(macsPeaks_GRL)

GRangesList 目标的一个主要长处是咱们能够将许多 GRanges 访问器和运算符函数直接运用于咱们的 GRangesList。

这意味着假如咱们希望通过通用办法更改咱们的 GRanges，则无需运用并转换回 GRangesList。

library(rtracklayer)
macsPeaks_GRLCentred <- resize(macsPeaks_GRL, 10, fix = "center")
width(macsPeaks_GRLCentred)

现在咱们有了 GRangesList，咱们能够提取 Mel 仿制的峰值调用。

Mel_1_Peaks <- macsPeaks_GRL$Mel_1
Mel_2_Peaks <- macsPeaks_GRL$Mel_2
length(Mel_1_Peaks)  # ## [1] 13777
length(Mel_2_Peaks)  # ## [1] 13512

3. 寻找 unique peaks

咱们能够运用 %over% 运算符提取仅有的峰值调用以仿制 1 或 2。

Mel_1_Unique <- Mel_1_Peaks[!Mel_1_Peaks %over% Mel_2_Peaks]
Mel_2_Unique <- Mel_2_Peaks[!Mel_2_Peaks %over% Mel_1_Peaks]
length(Mel_1_Unique)  # ## [1] 4668
length(Mel_2_Unique)  # ## [1] 4263
export.bed(Mel_1_Unique, "Mel_1_Unique.bed")
export.bed(Mel_2_Unique, "Mel_2_Unique.bed")

4. 寻找 common peaks

同样，咱们能够提取仿制 1 或 2 常见的峰值调用。

然而，一起的数字不同。这是由于一个样本中的 2 个峰调用能够与另一个重复中的 1 个峰调用堆叠。

Mel_1_Common <- Mel_1_Peaks[Mel_1_Peaks %over% Mel_2_Peaks]
Mel_2_Common <- Mel_2_Peaks[Mel_2_Peaks %over% Mel_1_Peaks]
length(Mel_1_Common)  # ## [1] 9109
length(Mel_2_Common)  # ## [1] 9249
export.bed(Mel_1_Common, "Mel_1_Common.bed")
export.bed(Mel_2_Common, "Mel_2_Common.bed")

尽管堆叠，但这些峰并不相同。那么咱们如何确定几个样本的一起共识峰。

5. 界说consensus, redundant 集

为了解决这个问题，ChIPseq 中的一个常见操作是在一切样本中界说一组非冗余峰。

为此，咱们首要将一切重复的一切峰（这里是 Mel 和 Ch12）汇集到一组冗余的堆叠峰中。

allPeaksSet_Overlapping <- unlist(macsPeaks_GRL)
allPeaksSet_Overlapping

然后咱们能够运用 reduce() 函数将咱们的峰折叠成非冗余的、不同的峰，代表任何样本中存在的峰。

allPeaksSet_nR <- reduce(allPeaksSet_Overlapping)
allPeaksSet_nR

export.bed(allPeaksSet_nR, "allPeaksSet_nR.bed")

6. 界说 common peaks

运用咱们新界说的非冗余峰集，咱们现在能够运用 %over% 运算符和逻辑表达式从该集中识别咱们的重复中存在哪些峰。

commonPeaks <- allPeaksSet_nR[allPeaksSet_nR %over% Mel_1_Peaks & allPeaksSet_nR %over%
    Mel_2_Peaks]
commonPeaks

export.bed(commonPeaks, "commonPeaks.bed")

7. 界说 unique peaks

同样，咱们能够确定哪些峰仅呈现在一次重复中。

mel1_Only <- allPeaksSet_nR[allPeaksSet_nR %over% Mel_1_Peaks & !allPeaksSet_nR %over%
    Mel_2_Peaks]
mel2_Only <- allPeaksSet_nR[!allPeaksSet_nR %over% Mel_1_Peaks & allPeaksSet_nR %over%
    Mel_2_Peaks]
length(mel1_Only)  # ## [1] 4445
length(mel2_Only)  # ## [1] 4185
export.bed(mel1_Only, "mel1_Only.bed")
export.bed(mel2_Only, "mel2_Only.bed")

8. 杂乱堆叠

当处理大量峰时，咱们通常会界说一个逻辑矩阵来描绘咱们的非冗余峰呈现在哪些样本中。

首要，咱们运用循环为每个样本中呈现的非冗余峰生成一个逻辑向量。

overlap <- list()
for (i in 1:length(macsPeaks_GRL)) {
    overlap[[i]] <- allPeaksSet_nR %over% macsPeaks_GRL[[i]]
}
overlap[[1]][1:2]

咱们现在能够运用 to do.call 和 cbind 函数将咱们的堆叠列表列绑定到咱们的峰值呈现矩阵中。

overlapMatrix <- do.call(cbind, overlap)
colnames(overlapMatrix) <- names(macsPeaks_GRL)
overlapMatrix[1:2, ]

咱们能够运用 mcols() 访问器将矩阵添加回非冗余峰的 GRanges() 的元数据列中。

现在咱们有了非冗余峰以及每个样本中这些峰的呈现，咱们能够轻松识别重复和条件/细胞系特有或共有的峰。

mcols(allPeaksSet_nR) <- overlapMatrix
allPeaksSet_nR[1:2, ]

limma 包常用于分析 RNAseq 和微阵列数据，并包括许多有用的功能。

一个十分有用的函数是 vennDiagram 函数，它答应咱们制作逻辑矩阵的堆叠，就像咱们创立的那样。

library(limma)
vennDiagram(mcols(allPeaksSet_nR))

limma 包的 vennCounts 函数答应咱们检索维恩图中显示的计数作为 data.frame。

vennCounts(mcols(allPeaksSet_nR))

9. 高置信度峰

运用咱们的非冗余峰集和峰呈现矩阵，咱们能够在条件下界说仿制峰。在这里，咱们界说了在两个 Ch12 重复中呈现的峰值。

由于逻辑矩阵等效于 1 或 0 矩阵（1 = TRUE 和 0 = FALSE），咱们能够运用 rowSums 函数在至少 2 个 Ch12 重复中提取峰。

ch12_HC_Peaks <- allPeaksSet_nR[rowSums(as.data.frame(mcols(allPeaksSet_nR)[, c("ch12_1",
    "ch12_2")])) >= 2]
export.bed(ch12_HC_Peaks, "ch12_HC_Peaks.bed")
ch12_HC_Peaks[1:2, ]

10. 高置信度的仅有峰

同样，咱们能够界说在 Ch12 中仿制但在 Mel 样本中不存在的峰。

ch12_HC_UniquePeaks <- allPeaksSet_nR[rowSums(as.data.frame(mcols(allPeaksSet_nR)[,
    c("ch12_1", "ch12_2")])) >= 2 & rowSums(as.data.frame(mcols(allPeaksSet_nR)[,
    c("Mel_1", "Mel_2")])) == 0]
export.bed(ch12_HC_UniquePeaks, "ch12_HC_UniquePeaks.bed")
ch12_HC_UniquePeaks[1, ]

本文由mdnice多渠道发布