1. 結(jié)果處理

現(xiàn)在我們已經(jīng)處理了 Greenleaf ATACseq 雙端數(shù)據(jù)，我們可以開始處理比對(duì)。

首先，我們將確定 ATACseq 數(shù)據(jù)的預(yù)期片段長(zhǎng)度分布。我們使用 GenomicAlignments 包讀取新對(duì)齊的數(shù)據(jù)。

這里我們只想要正確配對(duì)的讀取，因此我們將使用 ScanBamParam() 和 scanBamFlag() 函數(shù)來控制將讀入 R 的內(nèi)容。

我們將 scanBamFlag() 函數(shù)參數(shù) isProperPair 設(shè)置為 TRUE，以便僅讀取在我們預(yù)設(shè)的最大片段長(zhǎng)度 (2000bpp) 內(nèi)對(duì)齊配對(duì)的讀數(shù)

library(GenomicAlignments)
flags = scanBamFlag(isProperPair = TRUE)

我們現(xiàn)在可以將這些標(biāo)志與 ScanBamParam() 函數(shù)一起使用，以僅讀入正確配對(duì)的讀數(shù)。

我們還使用 what 參數(shù)指定要讀入 R 的信息。重要的是，我們指定插入大小信息 - isize。為了減少內(nèi)存占用，我們通過指定 GRanges 對(duì)象的參數(shù)來只讀取來自 20 號(hào)染色體的信息。

myParam = ScanBamParam(flag = flags, what = c("qname", "mapq", "isize"), which = GRanges("chr20",
    IRanges(1, 63025520)))
myParam

現(xiàn)在我們已經(jīng)設(shè)置了 ScanBamParam 對(duì)象，我們可以使用 readGAlignmentPairs() 函數(shù)以類似于我們使用 readGAlignments() 函數(shù)讀取單端 ChIP-seq 數(shù)據(jù)的方式讀取我們的雙端 ATACseq 數(shù)據(jù)。結(jié)果是一個(gè) GAlignmentPairs 對(duì)象。

atacReads <- readGAlignmentPairs(sortedBAM, param = myParam)
class(atacReads)

2. GAlignmentPairs

GAlignmentPairs 對(duì)象包含有關(guān)我們配對(duì)讀取的信息。它將每次讀取的信息成對(duì)存儲(chǔ)在并行的 GAlignments 對(duì)象中。

atacReads[1:2, ]

我們使用 first() 和 second() 訪問器函數(shù)訪問 GAlignments 對(duì)象，以分別獲取有關(guān)第一次或第二次讀取的信息。

read1 <- first(atacReads)
read2 <- second(atacReads)
read2[1, ]

3. MapQ 分?jǐn)?shù)

我們可以做的第一件事是獲取 read1 和 read2 的 MapQ 分?jǐn)?shù)分布。我們可以使用 mcols() 函數(shù)為每次讀取訪問它，以訪問每次讀取的 GAalignments 對(duì)象的 mapq 槽。

read1MapQ <- mcols(read1)$mapq
read2MapQ <- mcols(read2)$mapq
read1MapQ[1:2]

4. MapQ 頻率

然后我們可以使用 table() 函數(shù)來匯總每個(gè)成對(duì)讀取的分?jǐn)?shù)頻率。

read1MapQFreqs <- table(read1MapQ)
read2MapQFreqs <- table(read2MapQ)
read1MapQFreqs

read2MapQFreqs

5. 可視化

最后，我們可以使用 ggplot2 繪制成對(duì)讀取的每個(gè) MapQ 分布。

library(ggplot2)
toPlot <- data.frame(MapQ = c(names(read1MapQFreqs), names(read2MapQFreqs)), Frequency = c(read1MapQFreqs,
    read2MapQFreqs), Read = c(rep("Read1", length(read1MapQFreqs)), rep("Read2",
    length(read2MapQFreqs))))
toPlot$MapQ <- factor(toPlot$MapQ, levels = unique(sort(as.numeric(toPlot$MapQ))))
ggplot(toPlot, aes(x = MapQ, y = Frequency, fill = MapQ)) + geom_bar(stat = "identity") +
    facet_grid(~Read)

6. 插入大小

現(xiàn)在我們已經(jīng)將配對(duì)的對(duì)齊數(shù)據(jù)讀入 R，我們可以從附加到每個(gè)讀取對(duì)的 GAlignments 對(duì)象的 elementMetadata() 中檢索插入大小。由于正確配對(duì)的讀取將具有相同的插入大小長(zhǎng)度，因此我們從 read1 中提取插入大小。

atacReads_read1 <- first(atacReads)
insertSizes <- abs(elementMetadata(atacReads_read1)$isize)
head(insertSizes)

7. 可視化

ATACseq 應(yīng)該代表對(duì)應(yīng)于無核小體、單核小體和多核小體部分的片段長(zhǎng)度的混合。我們可以使用 table() 函數(shù)來檢索每個(gè)片段長(zhǎng)度出現(xiàn)的向量。

fragLenSizes <- table(insertSizes)
fragLenSizes[1:5]

現(xiàn)在我們可以使用新獲得的 20 號(hào)染色體插入長(zhǎng)度來繪制所有片段長(zhǎng)度的分布。

library(ggplot2)
toPlot <- data.frame(InsertSize = as.numeric(names(fragLenSizes)), Count = as.numeric(fragLenSizes))
fragLenPlot <- ggplot(toPlot, aes(x = InsertSize, y = Count)) + geom_line()
fragLenPlot + theme_bw()

我們可以對(duì)計(jì)數(shù)應(yīng)用 log2 轉(zhuǎn)換以闡明核小體模式。

fragLenPlot + scale_y_continuous(trans = "log2") + theme_bw()

我們現(xiàn)在可以像 Greenleaf 研究中那樣注釋我們的核小體游離 (< 100bp)、單核小體 (180bp-247bp) 和雙核小體 (315-437)。

fragLenPlot + scale_y_continuous(trans = "log2") + geom_vline(xintercept = c(180,
    247), colour = "red") + geom_vline(xintercept = c(315, 437), colour = "darkblue") +
    geom_vline(xintercept = c(100), colour = "darkgreen") + theme_bw()

歡迎Star -> 學(xué)習(xí)目錄

更多教程 -> 轉(zhuǎn)錄組測(cè)序分析教程合集

更多教程 -> 單細(xì)胞系列教程：合集

本文由mdnice多平臺(tái)發(fā)布