关于程序员:ChIPseq-分析Differential-Peaks15

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

然而，辨认特定于细胞系或条件的峰并不能捕捉表观遗传事件的全副变动。

为了辨认表观遗传事件的差别，咱们能够尝试量化 IP 样本中非冗余峰组中片段丰度的变动。

咱们首先必须建设一组区域，在这些区域中量化 IP ed 片段。

一种成熟的技术是产生一组非冗余峰，这些峰呈现在至多一个被评估的试验条件的大多数中。

在这里，咱们确定了在 Mel 或 Ch12 细胞系的两个反复中呈现的峰。

HC_Peaks <- 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")])) >= 2]
HC_Peaks

export.bed(HC_Peaks, "HC_Peaks.bed")

咱们将从对齐的 BAM 文件中计数以量化 IP 片段。

正如咱们之前所见，咱们能够应用 BamFileList() 函数来指定要计数的 BAM，重要的是，为了管制内存，咱们应用 yield() 参数指定一次要保留在内存中的读取次数。

library(Rsamtools)

bams <- c("~/Projects/Results/chipseq/testRun/BAMs/Sorted_Myc_Ch12_1.bam", "~/Projects/Results/chipseq/testRun/BAMs/Sorted_Myc_Ch12_2.bam",
    "~/Projects/Results/chipseq/testRun/BAMs/Sorted_Myc_Mel_1.bam", "~/Projects/Results/chipseq/testRun/BAMs/Sorted_Myc_Mel_2.bam")
bamFL <- BamFileList(bams, yieldSize = 5e+06)
bamFL

咱们能够应用 summarizeOverlaps 函数计算与峰重叠的片段数。因为 ChIPseq 是无链的，咱们将 ignore.strand 参数设置为 TRUE。

返回的对象是一个相熟的 RangedSummarizedExperiment，其中蕴含咱们的非冗余峰的 GRanges 以及咱们 BAM 文件在这些区域中的计数。

library(GenomicAlignments)
myMycCounts <- summarizeOverlaps(HC_Peaks, reads = bamFL, ignore.strand = TRUE)
class(myMycCounts)
save(myMycCounts, file = "data/MycCounts.RData")

为了评估跨细胞系的 ChIPseq 信号变动，咱们将应用 DESeq2 包。

DESeq2 包蕴含一个工作流程，用于评估复制条件之间片段 / 读取丰度的部分变动。此工作流程包含标准化、方差预计、离群值移除 / 替换以及实用于高通量测序数据（即整数计数）的显着性测试。

要应用 DESeq2 工作流程，咱们必须首先创立一个感兴趣条件的 data.frame，并将行名设置为咱们的 BAM 文件名。

metaDataFrame <- data.frame(CellLine = c("Ch12", "Ch12", "Mel", "Mel"))
rownames(metaDataFrame) <- colnames(myMycCounts)
metaDataFrame

咱们能够应用 DESeqDataSetFromMatrix() 函数来创立 DESeq2 对象。

咱们必须将咱们的计数矩阵提供给 countData 参数，咱们的元数据 data.frame 提供给 colData 参数，并且咱们在 rowRanges 的可选参数中蕴含咱们能够依附的非冗余峰值集。

最初，咱们在咱们心愿测试设计参数的元数据 data.frame 中提供列的名称。

library(DESeq2)
deseqMyc <- DESeqDataSetFromMatrix(countData = assay(myMycCounts), colData = metaDataFrame,
    design = ~CellLine, rowRanges = HC_Peaks)

咱们当初能够应用 DESeq() 函数在咱们的 DESeq2 对象上运行 DESeq2 工作流程。

deseqMyc <- DESeq(deseqMyc)

咱们的 DESeq2 对象已更新，以蕴含有用的统计信息，例如咱们的标准化值和每个非冗余峰调用中的信号方差。

deseqMyc

咱们能够应用 results() 函数提取差别区域的信息。

咱们向 results() 函数提供 DESeq2 对象、对比照参数感兴趣的比拟以及返回格局参数的输入类型。

与比照参数的比拟作为长度为 3 的向量提供，包含感兴趣的元数据列和要测试的组。

咱们能够应用 order() 函数按 pvalue 对后果进行排序，以按最显着的变动进行排名。

MelMinusCh12 <- results(deseqMyc, contrast = c("CellLine", "Mel", "Ch12"), format = "GRanges")
MelMinusCh12 <- MelMinusCh12[order(MelMinusCh12$pvalue), ]
class(MelMinusCh12)

GRanges 对象蕴含无关在 DESeq2 中进行的比拟的信息。

最有用的是它蕴含 IP 信号的差别，如 log2FoldChange 中的 log2 倍变动、pvalue 列中变动的重要性以及调整后的 p 值以解决 padj 列中的多重校对。

MelMinusCh12[1,]

咱们当初能够通过过滤 log2FoldChange 和 padj（针对多重校对调整的 p 值）小于 0.05，将咱们的非冗余峰过滤为 Mel 或 Ch12 细胞系中信号显著更多的峰。

MelMinusCh12Filt <- MelMinusCh12[!is.na(MelMinusCh12$pvalue) | !is.na(MelMinusCh12$padj)]
UpinMel <- MelMinusCh12[MelMinusCh12$padj < 0.05 & MelMinusCh12$log2FoldChange >
    0]
DowninMel <- MelMinusCh12[MelMinusCh12$padj < 0.05 & MelMinusCh12$log2FoldChange <
    0]
export.bed(UpinMel, "UpinMel.bed")
export.bed(DowninMel, "DowninMel.bed")

最初，咱们能够应用 tracktables 包让咱们在 IGV 中审查站点更容易一些。

tracktables 包的 makebedtable() 函数承受一个 GRanges 对象并编写一个蕴含 IGV 链接的 HTML 报告。

library(tracktables)
myReport <- makebedtable(MelMinusCh12Filt, "MelMinusCh12.html", basedirectory = getwd())

browseURL(myReport)

本文由 mdnice 多平台公布

关于程序员:ChIPseq-分析Differential-Peaks15

1. 寻找差别区域

2. 计数区域

3. DESeq2