動(dòng)動(dòng)發(fā)財(cái)?shù)男∈?，點(diǎn)個(gè)贊吧！
1. 包加載
我們可以使用 rGREAT 包中提供的 GREAT Bioconductor 接口。
library(rGREAT)
2. GO和功能測(cè)試
要提交作業(yè)，我們可以使用 Myc 峰的 GRanges 并使用 submitGreatJob 函數(shù)指定基因組。
此函數(shù)返回一個(gè) GreatJob 對(duì)象，其中包含對(duì)我們?cè)?GREAT 服務(wù)器上的結(jié)果的引用。要查看可用結(jié)果的類別，我們可以在 GreatJob 對(duì)象上使用 availableCategories 函數(shù)。
great_Job <- submitGreatJob(macsPeaks_GR, species = "mm10", version = "3.0.0", request_interval = 1)
availableCategories(great_Job)
可以使用 getEnrichmentTables 函數(shù)檢索結(jié)果表并指定我們希望查看的表。
在這里，我們檢索包含 2 個(gè)獨(dú)立數(shù)據(jù)庫(kù)結(jié)果的“Regulatory Motifs”基因集的結(jié)果表。
great_ResultTable = getEnrichmentTables(great_Job, category = "Regulatory Motifs")
names(great_ResultTable)
現(xiàn)在我們可以在“MSigDB 預(yù)測(cè)的啟動(dòng)子基序”基因集的 TSS 中使用 Myc 峰查看我們的基因的富集情況。
msigProMotifs <- great_ResultTable[["MSigDB Predicted Promoter Motifs"]]
msigProMotifs[1:4, ]
3. Motifs 分析
3.1. Motifs
轉(zhuǎn)錄因子 ChIPseq 的一個(gè)常見做法是研究峰下富集的基序?？梢栽?R/Bioconductor 中進(jìn)行從頭富集基序，但這可能非常耗時(shí)。在這里，我們將使用在線提供的 MEME-ChIP 套件來(lái)識(shí)別新的基序。
MEME-ChIP 需要一個(gè)包含峰下序列的 FASTA 文件作為輸入，因此我們使用 BSgenome 包提取它。
3.2. 序列提取
首先，我們需要為我們正在處理的基因組加載 BSgenome 對(duì)象，UCSC 為小鼠基因組構(gòu)建的 mm10，BSgenome.Mmusculus.UCSC.mm10。
library(BSgenome)
library(BSgenome.Mmusculus.UCSC.mm10)
BSgenome.Mmusculus.UCSC.mm10
我們現(xiàn)在有一個(gè) GRanges，以山頂為中心，每個(gè)山峰的最高信號(hào)點(diǎn)。
macsSummits_GR
一旦我們使峰重新居中，我們就可以將 getSeq 函數(shù)與調(diào)整大小的常見峰的 GRanges 和 mm10 的 BSgenome 對(duì)象一起使用。
getSeq 函數(shù)返回包含峰下序列的 DNAStringSet 對(duì)象。
peaksSequences <- getSeq(BSgenome.Mmusculus.UCSC.mm10, macsSummits_GR)
names(peaksSequences) <- paste0(seqnames(macsSummits_GR), ":", start(macsSummits_GR),
"-", end(macsSummits_GR))
peaksSequences[1:2, ]
3.3. 寫入 FASTA 文件
writeXStringSet 函數(shù)允許用戶將 DNA/RNA/AA（氨基酸）StringSet 對(duì)象寫入文件。默認(rèn)情況下，writeXStringSet 函數(shù)以 FASTA 格式寫入序列信息（根據(jù) MEME-ChIP 的要求）。
writeXStringSet(peaksSequences, file = "mycMel_rep1.fa")
3.4. MEME-ChIP
現(xiàn)在文件“mycMel_rep1.fa”包含適合 MEME-ChIP 中 Motif 分析的峰幾何中心周圍的序列。
在您自己的工作中，您通常會(huì)在本地安裝了 MEME 的筆記本電腦上運(yùn)行它，但今天我們會(huì)將生成的 FASTA 文件上傳到他們的門戶網(wǎng)站。按照此處的說明在本地安裝 MEME?？梢栽?a target="_blank">此處找到 MEME-ChIP 的結(jié)果文件
3.5. 結(jié)果解析
我們可以從 FIMO 輸出中檢索 MEME-ChIP 中識(shí)別的 Myc 基序的位置。
FIMO 將 Myc 基序位置報(bào)告為 GFF3 文件，我們應(yīng)該能夠在 IGV 中對(duì)其進(jìn)行可視化。遺憾的是，這個(gè) GFF 文件的命名約定只導(dǎo)致報(bào)告了一小部分圖案。
3.6. FIMO to R
幸運(yùn)的是，我們可以將 motif 的 GFF 文件解析為 R 并使用 rtracklayer 包中的導(dǎo)入函數(shù)解決這個(gè)問題。
library(rtracklayer)
motifGFF <- import("~/Downloads/fimo.gff")
3.7. 獲取有效 GFF3
我們可以給序列一些更合理的名稱并將 GFF 導(dǎo)出到文件以在 IGV 中可視化。
motifGFF$Name <- paste0(seqnames(motifGFF), ":", start(motifGFF), "-", end(motifGFF))
motifGFF$ID <- paste0(seqnames(motifGFF), ":", start(motifGFF), "-", end(motifGFF))
export.gff3(motifGFF, con = "~/Downloads/fimoUpdated.gff")
3.8. 掃描已知 motifs
我們之前看到我們可以使用一些 Biostrings 功能 matchPattern 來(lái)掃描序列。通常使用 ChIPseq，我們可能知道我們正在尋找的基序，或者我們可以使用來(lái)自數(shù)據(jù)庫(kù)（例如 JASPAR）的一組已知基序。
library(JASPAR2020)
JASPAR2020
3.9. 使用 TFBStools 從 JASPAR 獲取 motifs
我們可以使用 TFBSTools 包及其 getMatrixByName 函數(shù)訪問我們感興趣的motif的模型。
library(TFBSTools)
pfm <- getMatrixByName(JASPAR2020, name = "MYC")
pfm
3.10. 使用 motifmathr 進(jìn)行 motifs 掃描
有了這個(gè) PWM，我們可以使用 motifmathr 包來(lái)掃描我們的山峰以尋找 Myc motif并返回motif的位置。
我們需要提供我們的 PWM、要在內(nèi)部掃描的 GRanges 和要從中提取序列的 BSGenome 對(duì)象。我們還將輸出參數(shù)設(shè)置為這個(gè)實(shí)例的位置。
library(motifmatchr)
MycMotifs <- matchMotifs(pfm, macsSummits_GR, BSgenome.Mmusculus.UCSC.mm10, out = "positions")
MycMotifs
3.11. 導(dǎo)出匹配的 motifs
我們可以導(dǎo)出峰內(nèi)的 Myc 基序位置，以便稍后在 IGV 中使用或用于元圖可視化。
export.bed(MycMotifs[[1]], con = "MycMotifs.bed")
本文由mdnice多平臺(tái)發(fā)布
