关于程序员:ChIPseq-分析Consensus-Peaks14

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

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 多平台公布