共计 2841 个字符,预计需要花费 8 分钟才能阅读完成。
动动发财的小手,点个赞吧!
1. 包加载
咱们能够应用 rGREAT 包中提供的 GREAT Bioconductor 接口。
library(rGREAT)
2. GO 和功能测试
要提交作业,咱们能够应用 Myc 峰的 GRanges 并应用 submitGreatJob 函数指定基因组。
此函数返回一个 GreatJob 对象,其中蕴含对咱们在 GREAT 服务器上的后果的援用。要查看可用后果的类别,咱们能够在 GreatJob 对象上应用 availableCategories 函数。
great_Job <- submitGreatJob(macsPeaks_GR, species = "mm10", version = "3.0.0", request_interval = 1)
availableCategories(great_Job)
能够应用 getEnrichmentTables 函数检索后果表并指定咱们心愿查看的表。
在这里,咱们检索蕴含 2 个独立数据库后果的“Regulatory Motifs”基因集的后果表。
great_ResultTable = getEnrichmentTables(great_Job, category = "Regulatory Motifs")
names(great_ResultTable)
当初咱们能够在“MSigDB 预测的启动子基序”基因集的 TSS 中应用 Myc 峰查看咱们的基因的富集状况。
msigProMotifs <- great_ResultTable[["MSigDB Predicted Promoter Motifs"]]
msigProMotifs[1:4,]
3. Motifs 剖析
3.1. Motifs
转录因子 ChIPseq 的一个常见做法是钻研峰下富集的基序。能够在 R/Bioconductor 中进行从头富集基序,但这可能十分耗时。在这里,咱们将应用在线提供的 MEME-ChIP 套件来辨认新的基序。
MEME-ChIP 须要一个蕴含峰下序列的 FASTA 文件作为输出,因而咱们应用 BSgenome 包提取它。
3.2. 序列提取
首先,咱们须要为咱们正在解决的基因组加载 BSgenome 对象,UCSC 为小鼠基因组构建的 mm10,BSgenome.Mmusculus.UCSC.mm10。
library(BSgenome)
library(BSgenome.Mmusculus.UCSC.mm10)
BSgenome.Mmusculus.UCSC.mm10
咱们当初有一个 GRanges,以山顶为核心,每个山峰的最高信号点。
macsSummits_GR
一旦咱们使峰从新居中,咱们就能够将 getSeq 函数与调整大小的常见峰的 GRanges 和 mm10 的 BSgenome 对象一起应用。
getSeq 函数返回蕴含峰下序列的 DNAStringSet 对象。
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 函数容许用户将 DNA/RNA/AA(氨基酸)StringSet 对象写入文件。默认状况下,writeXStringSet 函数以 FASTA 格局写入序列信息(依据 MEME-ChIP 的要求)。
writeXStringSet(peaksSequences, file = "mycMel_rep1.fa")
3.4. MEME-ChIP
当初文件“mycMel_rep1.fa”蕴含适宜 MEME-ChIP 中 Motif 剖析的峰几何核心四周的序列。
在您本人的工作中,您通常会在本地装置了 MEME 的笔记本电脑上运行它,但明天咱们会将生成的 FASTA 文件上传到他们的门户网站。依照此处的阐明在本地装置 MEME。能够在此处找到 MEME-ChIP 的后果文件
3.5. 后果解析
咱们能够从 FIMO 输入中检索 MEME-ChIP 中辨认的 Myc 基序的地位。
FIMO 将 Myc 基序地位报告为 GFF3 文件,咱们应该可能在 IGV 中对其进行可视化。遗憾的是,这个 GFF 文件的命名约定只导致报告了一小部分图案。
3.6. FIMO to R
侥幸的是,咱们能够将 motif 的 GFF 文件解析为 R 并应用 rtracklayer 包中的导入函数解决这个问题。
library(rtracklayer)
motifGFF <- import("~/Downloads/fimo.gff")
3.7. 获取无效 GFF3
咱们能够给序列一些更正当的名称并将 GFF 导出到文件以在 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 来扫描序列。通常应用 ChIPseq,咱们可能晓得咱们正在寻找的基序,或者咱们能够应用来自数据库(例如 JASPAR)的一组已知基序。
library(JASPAR2020)
JASPAR2020
3.9. 应用 TFBStools 从 JASPAR 获取 motifs
咱们能够应用 TFBSTools 包及其 getMatrixByName 函数拜访咱们感兴趣的 motif 的模型。
library(TFBSTools)
pfm <- getMatrixByName(JASPAR2020, name = "MYC")
pfm
3.10. 应用 motifmathr 进行 motifs 扫描
有了这个 PWM,咱们能够应用 motifmathr 包来扫描咱们的山峰以寻找 Myc motif 并返回 motif 的地位。
咱们须要提供咱们的 PWM、要在外部扫描的 GRanges 和要从中提取序列的 BSGenome 对象。咱们还将输入参数设置为这个实例的地位。
library(motifmatchr)
MycMotifs <- matchMotifs(pfm, macsSummits_GR, BSgenome.Mmusculus.UCSC.mm10, out = "positions")
MycMotifs
3.11. 导出匹配的 motifs
咱们能够导出峰内的 Myc 基序地位,以便稍后在 IGV 中应用或用于元图可视化。
export.bed(MycMotifs[[1]], con = "MycMotifs.bed")
本文由 mdnice 多平台公布