1概述

igvR包提供了简单的编程访问R的基于web的javascript库igv.js。R的数据访问和操作数据,辅以igv丰富的交互式web浏览器界面显示和探索基因组和epigentic数据。

在这个描述我们现在有些刻意ChIP-seq研究——做作,它不是围绕一个实际的研究问题。相反,您可能会使用这个小插图展示了方法做视觉QC和探索性数据分析ChIP-seq实验。

我们开始与CTCF的编码数据绑定3号染色体上,hg19参考基因组,附近的GATA2基因。我们将

• 显示一个252 kb的基因组区域
• 加载和显示一个GenomicAlignments对象从一个索引bam文件覆盖地区
• 加载和显示一个狭窄的山峰bedGraph跟踪,由MACS2预先计算的
• 使用Biostrings:: matchPWM方法识别区域显示的序列中,至少有80%的忠诚,与CTCF的主题,从MotifDb Jaspar 2018年MA0138.1
• 浏览器弹出显示图案标志。
• 负载保护包phastCons100way.UCSC Bioconductor注释的数据序列。hg19,寻找芯片的交集山峰,主题匹配和保守序列,所有这些一起显示,功能性TF结合发生在这样的位置。

下面的两张图片说明进口的最终状态在你的浏览器下面的代码发现的结论。第一个视图跨越35 kb:

第二个跨度小于1200基地。

2我们需要加载库

库(igvR)库(MotifDb)图书馆(AnnotationHub)图书馆(BSgenome.Hsapiens.UCSC.hg19)库(phastCons100way.UCSC.hg19)

创建igvR实例,所有默认参数(portRange,安静,标题)。Javascript和HTML加载到浏览器中,进口。js是初始化,websocket连接之间R过程和构造web页面,在后续的命令和数据将旅行。

igv < - igvR () setBrowserWindowTitle (igv CTCF ChIP-seq) setGenome(进口,“hg19”)

在chr3显示1.4 mb或多或少地集中在GATA2基因

showGenomicRegion (igv chr3:128,079,020 - 128331275)

3创建一个农庄对象指定我们的地区的利益,获得一个列表(铬,开始、结束字符串)来自进口

loc < - getGenomicRegion(进口)<——(loc,农庄(seqnames =铬,范围= IRanges(开始、结束)))参数< - ScanBamParam (=, = scanBamWhat ()) bamFile < -“~ / github / ChIPseqMotifMatch /体积/ GSM749704 / GSM749704_hg19_wgEncodeUwTfbsGm12878CtcfStdAlnRep1。bam”file.exists (bamFile) x < - readGAlignments (bamFile、use.names = TRUE参数=参数)跟踪< - GenomicAlignmentTrack (ctcf bam, x, visibilityWindow = 10000000, trackHeight = 200) # 30000默认displayTrack(进口、跟踪)

4负载窄峰值之前被MACS2 bam的文件

我们以前跑MACS2使用码头工人,使,这个makefile,一个方法,你可能会发现有用的但要求您安装码头工人和MACS2形象。由此产生的窄峰文件是包含在包。

Macs2 narrowpeaks产生一个可能的总结和bam连环相撞的减少,尽管bam跟踪直方图可以说是包含更多的信息。

DATADIR = /用户/保罗/ github / igvR /小插曲/ macs2 BAM = GSM749704_hg19_wgEncodeUwTfbsGm12878CtcfStdAlnRep1。bam NAME = GSM749704_hg19_chr19运行:码头工人跑- v (DATADIR):美元/数据/ fooliu macs2 callpeak - t /数据/ (bam) - n(名字)——美元outdir数据/

这产生多个文件;我们只有在狭窄的山峰:感兴趣

文件名< -系统。文件(包=“igvR”、“extdata”、“GSM749704_hg19_chr19_peaks.narrowPeak”)资源。pk < - get(负载(系统。文件(包=“igvR”、“extdata”、“GSM749704_hg19_chr19_subset.RData”)))跟踪< - DataFrameQuantitativeTrack (“macs2峰”,资源描述。pk,颜色=“红色”,自动定量= TRUE) displayTrack(进口、跟踪)

5与CTCF主题对序列的当前显示区域

< - dna (loc, getSeq (BSgenome.Hsapiens.UCSC。hg19、铬、开始、结束)烤瓷。ctcf < - MotifDb::查询(MotifDb c (“ctcf”、“智人”,“jaspar2018”), notStrings =“ctcfl”) motif.name < -名字(pfm.ctcf)[1]烤瓷< - pfm.ctcf [[1]]

6使用从Biostrings matchPWM包匹配两个链上的图案

我们用宽容匹配阈值反映了我们的观点,TF绑定通常是相当宽松的,不准确。

点击。转发< - matchPWM(烤瓷,as.character (dna)。分数= TRUE, min.score =“80%”)。反向< - matchPWM (reverseComplement (pfm) as.character (dna)。分数= TRUE, min.score =“80%”)资源。< - as.data.frame向前台(范围(hits.forward))。反向< - as.data.frame(范围(hits.reverse))(资源。美元远期美元得分< - mcols (hits.forward)分数(资源。反向评分< - mcols美元美元(hits.reverse)分数(资源。匹配< - rbind(资源描述。向前,tbl.reverse)资源。铬< - loc铬资源美元相匹配。美元开始< -资源匹配。比赛开始+ loc美元开始资源。比赛结束美元< -台。比赛结束+ loc美元开始

7设置主题标志查看

它可能是有用的视图上下文中的标志匹配的主题。igvR支持这个,(目前)准备:

1. 告诉igvR你希望使用这个功能
2. 添加一个可辨认的“MotifDb::”协议字符串中的主题名称注释

启用了这个选项,名称(或标签)领域的一个注释床data.frame前缀MotifDb: / /,然后单击主题元素CTCF_MA0139.1轨道,相应的主题将显示在一个弹出窗口。这是特别有用,当你显示跟踪多个主题。

enableMotifLogoPopups(进口)(资源。< - $名称匹配paste0 (“MotifDb::”, motif.name)资源。< -资源匹配。匹配[c(“铬”、“开始”、“结束”,“名字”,“分数”)]暗(tbl.matches)

8呈现的主题匹配两种方法,有两条途径:一个注释,一个分数

在igvR我们区分注释跟踪,定量跟踪,其他跟踪(往下VCF和bam)。注释床的跟踪是一个表示格式。定量追踪代表假发的格式。

简而言之:想象你与这五个data.frame列:

1. 铬
2. 开始
3. 结束
4. 的名字
5. 分数

构造一个DataFrameAnnotationTrack列1 - 4。

构造一个DataFrameQuantitativeTrack列1、2、3、5。

如果名称(四)列包含MotifDb条目的名称(例如,“Hsapiens-jaspar2018-CTCF-MA0139.1”)然后预谋“MotifDb: / /”来支持主题图标弹出。

跟踪< - DataFrameAnnotationTrack (“CTCF-MA0139.1”,资源描述。匹配[1:4],颜色=“随机”,trackHeight = 25) displayTrack(进口、跟踪)跟踪< - DataFrameQuantitativeTrack (“MA0139.1分数”,资源描述。匹配(c(1、2、3、5)],颜色=“随机”,自动定量= FALSE, min = 8, max = 12) displayTrack(进口、跟踪)

9是否有理由相信这些主题匹配功能吗?

主题匹配的DNA序列的频率,尤其是当(如上所述),一个轻松的主题匹配阈值使用。然而,有三个因素,增加比赛的可能性识别功能结合位点:

1. 他们发现在或接近高信任度ChIP-seq网站
2. 跨物种DNA序列高度保守的和时间
3. 序列存在于一个已知的监管区域与目标相关联的基因(这里是GATA2)

第一个因素可以从浏览器中直接检查。第二,让我们添加phastCons100way基因组,Bioconductor注释包代表在100种经过数百万年的保守序列。

获得细粒度保护的评估中,我们使用箱子短于图案。我们添加一个保护跟踪显示每个底的保护两块。保护是在0和1之间。

loc < - getGenomicRegion(进口)开始< - (loc, seq(= 5)开始、结束)< -结束开始+ 5数< -长度(开始)资源。块< - data.frame(铬=代表(loc铬美元数),开始=开始,结束=结束,stringsAsFactors = FALSE)暗(tbl.blocks)资源。cons100 < - as.data.frame (gscores (phastCons100way.UCSC.hg19农庄(tbl.blocks)), stringsAsFactors = FALSE)资源。cons100铬< - as.character (tbl.cons100 seqnames美元)资源。cons100 < -台。cons100 [c(“铬”、“开始”、“结束”,“违约”)]跟踪< - DataFrameQuantitativeTrack (“cons100”,资源描述。cons100,自动定量= TRUE,颜色=“暗绿色”)displayTrack(进口、跟踪)

10添加甲基化跟踪

H3K4me3组蛋白甲基化标记往往与启动子的转录起始相关区域。我们将添加三个跟踪H3K4me3 AnnotationHub的甲基化Gm12878(lymphoblastoid)细胞。

啊< - AnnotationHub()啊。人类< -子集(物种啊,= =“智人”)组蛋白。歌曲< - AnnotationHub::查询(啊。人类,c (“H3K4me3”、“Gm12878”,“峰”,“窄”)< -组蛋白)# 3跟踪描述。跟踪美元描述标题< -组蛋白。跟踪美元标题颜色< -代表(terrain.colors(6), 4)颜色。指数< - 0 roi <——(getGenomicRegion(进口),农庄(seqnames =铬,IRanges(开始=开始,结束=))),(我在seq_len(长度(histone.tracks))){名称< -名称(histone.tracks)(我)的颜色。指数< -颜色。指数+ 1 gr < -组蛋白。跟踪[[名字]]ov < - findOverlaps (gr, roi)组蛋白< - gr (queryHits (ov)]。组蛋白< GRangesQuantitativeTrack(标题[我],组蛋白,“pValue”,颜色=颜色(颜色。指数),trackHeight = 50,自动定量= TRUE) displayTrack (igv track.histones) Sys.sleep(5)} #跟踪

showGenomicRegion (igv chr3:128,216,201 - 128217324)

11一些初步的结论

大ChIP-seq峰值位于中间的大基因内区GATA2-AS1是集中在一个好的搭配的CTCF Jaspar 2018主题,和大部分的序列匹配与高度保守的DNA。这三个特征的同现表明这个位置可能是CTCF的功能结合位点。

12会话信息

sessionInfo ()' ' ' '