转录因子结合位点(TFBS)分析与“TFBSTools”包
2022年4月26日
摘要
转录因子结合位点的分析与调控。包
TFBSTools 1.34.0
内容
1简介
真核调控区域的特征是基于一组已发现的转录因子结合位点(TFBSs)，这些转录因子结合位点可以表示为具有不同简并程度的序列模式。
这TFBSTools包被设计为TFBSs分析的计算框架。基于著名的perl模块TFBS(伦哈德和沃瑟曼2002)，我们在交互环境中扩展了类定义和增强了实现。到目前为止，这个包包含了一套集成RS4样式类，工具，JASPAR数据库接口函数。大多数方法可以分为三个连续的阶段。首先，为已知由特定转录因子结合的一组目标序列生成模式。其次，分析一组DNA序列，以确定与所述结合模式一致的序列位置。最后，在高级情况下，基于检测到的模式的多次出现，构建了调控区域的预测统计模型。
自JASPAR2016以来，下一代转录因子结合位点TFFM(Mathelier and Wasserman 2013)，第一次引入到JASPAR中。现在TFBSTools还支持的操作TFFM．TFFM基于隐马尔可夫模型(HMM)。与基本PWM相比，TFFM的最大优点是它可以模拟TFBSs内的位置相互依赖性和可变motif长度。一种新的图形表示的TFFM主题，捕捉位置相互依赖也被介绍。有关TFFM的详情，请参阅http://cisreg.cmmt.ubc.ca/TFFM/doc/．
TFBSTools旨在支持环境中的所有这些功能R，除了外部motif查找软件，如MEME(贝利和埃尔坎1994)．
2S4TFBSTools的班级
该包是围绕一些S4其类别为XMatrix，SiteSet课程是最重要的。本节将简要解释中定义的大多数TFBSTools．
2.1XMatrix及其子类
XMatrix是一个虚拟类，这意味着不能直接从它创建具体对象。子类PFMatrix用于存储一个原始位置频率矩阵(PFM)的所有相关信息。该对象与JASPAR数据库中的一条记录兼容。PWMatrix用于存储位置权重矩阵(PWM)。相比之下,PFMatrix，它有一个额外的插槽pseudocounts．ICMatrix用于存储信息内容矩阵(ICM)。相比之下,PWMatrix，它有一个额外的插槽施耐德．
下面的示例演示如何创建PFMatrix，这些矩阵和为这些类定义的一些相关方法之间的转换。
library(TFBSTools) ## PFMatrix构造;不是所有的槽都需要初始化。烤瓷< - PFMatrix (ID = " MA0004.1 ", name = " Arnt " matrixClass =“拉链”,链= " + ",bg = c (c = 0.25, = 0.25 G = 0.25, T = 0.25),标签=列表(家庭=“Helix-Loop-Helix物种= " 10090 ",tax_group =“脊椎动物”,medline =“7592839”,类型=“SELEX”,ACC =“P53762 pazar_tf_id =“TF0000003 TFBSshape_ID =“十一”,TFencyclopedia_ID = " 580 "), profileMatrix =矩阵(c (4 l, 19 l, l 0, 0, 0, 0 l, l, 16 0 l, 20 l, l 0, 0, 0, 0 l, l, 1 0 l, 20 l 0 l, 20 l, l 0, 0, 0, 0 l, 20 l 0 l), byrow = TRUE, nrow = 4, dimnames =列表(c(“A”、“c”、“G”," T "))))烤瓷# >的对象类PFMatrix # > ID: MA0004.1 # >名称:Arnt # >矩阵类:拉链# >链:+ # >标签:# > $家族# >[1]“Helix-Loop-Helix”# > # > $物种# >[1]“10090”# > # > $ tax_group # >[1]“脊椎动物”# > # > $ medline # >[1]“7592839”# > # > $类型# >[1]“SELEX”# > # > $ ACC # >[1]“P53762”# > # > $ pazar_tf_id # >[1]“TF0000003”# > # > $ TFBSshape_ID # >[1]“十一”# > # > $ TFencyclopedia_ID # >[1]“580”# > # >背景:# > C G T # > 0.25 - 0.25 0.25 - 0.25 # >矩阵:# > [1] [2] [3] [4] [5] [6] 19 # > 4 0 0 0 0 # > 20 C 16 0 0 0 0 # > 20 G 0 1 0 0 20 # > 20 T 0 0 0 0 0 # #强迫矩阵as.matrix (pfm) # > [1] [2] [3] [4] [5] [6] 19 # > 4 0 0 0 0 # > 20 C 16 0 0 0 0 # > 20 G 0 1 0 0 20 # > 20 T 0 0 0 0 0 # #访问烤瓷ID的插槽(pfm) # >[1]“MA0004.1”名称(pfm) # >[1]“Arnt”矩阵(pfm) # > [1] [2] [3] [4] [5] [6] 19 # > 4 0 0 0 0 # > 20 C 16 0 0 0 0 # > 20 G 0 1 0 0 20 # > 20 T 0 0 0 0 0 ncol (pfm) # >[1] 6长度(pfm) # >[1] 6 # #转化PFM toPWM, ICM PWM <- toPWM(PFM, type="log2probratio"， pseudocounts=0.8, bg=c(A=0.25, c =0.25, G=0.25, T=0.25)) ICM <- toICM(PFM, pseudocounts=sqrt(rowsum (PFM)[1])， schneider=FALSE, bg=c(A=0.25, c =0.25, G=0.25, T=0.25)) ##获取除链外所有相同信息的逆补矩阵。pwmRevComp <-逆补(pwm)2．2XMatrixList及其子类
XMatrixList是用来存储一套的吗XMatrix对象。基本上它是一个简单列表，便于操作整个集合XMatrix．具体的对象可以是PFMatrix，PWMatrix而且ICMatrix．
pfm2 <- pfm pfmList <- PFMatrixList(pfm1=pfm, pfm2=pfm2, use.names=TRUE) pfmList #> PFMatrixList长度为2 #>名称(2):pfm1 pfm2名称(pfmList) #> [1] "pfm1" "pfm2"2．3SiteSet, SiteSetList, SitePairSet和SitePairSetList
的SiteSet类是一个容器，用于存储核苷酸序列上的一组假定的转录因子结合位点(起始、结束、链、分数、模式等)PWMatrix等)从扫描核苷酸序列与相应的PWMatrix．同样的,SiteSetList存储一组SiteSet对象。
为了保存从成对对齐扫描返回的结果，SitePairSet而且SitePairSetList被提供。
后面一节将给出使用这些类的更详细示例。
2.4MotifSet
这MotifSet类用于存储生成的主题新创Motif发现软件等MEME(贝利和埃尔坎1994)．
2.5TFFM及其子类
TFMM是一个虚班和两个班TFFMFirst而且TFFMDetail派生自此虚拟类。相比之下,PFMatrix类,TFFM有两个额外的槽，存储发射分布参数和跃迁概率。TFFMFirstclass表示一阶TFFMs，而TFFMDetail表示更详细和描述性的tffm。
尽管我们提供了构造函数TFFM类,TFFM对象通常通过从Python模块读取XML文件而生成TFFM．
xmlFirst <- file.path(system. path)file("extdata"， package="TFBSTools")， "tffm_first_order.xml") tffmFirst <- readXMLTFFM(xmlFirst, type="First") tffm <- getPosProb(tffmFirst) xmlDetail <- file.path(system. path)file("extdata"， package="TFBSTools")， " tffm_details .xml") tffmDetail <- readXMLTFFM(xmlDetail，type="Detail") getPosProb(tffmDetail) #> 12 34 56 7 #> A 0.2114735 0.2838839 0.16377668 0.1871573 0.01873841 #> C 0.3347593 0.2854457 0.2441757 0.1840625 0.41600416 0.935778162 0.77794455 #> G 0.2278705 0.2736597 0.2133885 0.4208620 0.49585999 0.3786690 0.2079182 0.05131866 0.003256293 0.12982385 #> 89 10 11 12 13 14 #> A 0.0493708 0.07792452 0.128592905 0.003276879 0.04503578 0.2113065 #> C 0.3350808 0.432493640.1518607 0.078588841 0.067743421 0.51743026 0.4022825 #> g 0.1622687 0.40736410 0.3456973 0.78624245657 0.924717895 0.40148909 0.1895450 #> t 0.4532797 0.08221774 0.0371010 0.006572597 0.004261805 0.03604487 0.1968661 #> 15 #> a 0.3661852 #> c 0.2161336 #> g 0.2429433 #> t 0.17473783.数据库接口为JASPAR2014数据库
本节将演示如何操作JASPAR 2014数据库。JASPAR是转录因子dna结合偏好的集合，建模为矩阵。它们可以转化为PWMs，用于扫描基因组序列。JASPAR是唯一具有此范围的数据库，其中数据可以不受限制地使用(开源)。一个Bioconducto实验数据包JASPAR2014随JASPAR的每个发行版一起提供。
3．1搜索JASPAR2014数据库
这个搜索函数获取数据库中所有矩阵的矩阵数据，匹配由命名参数定义的标准，并返回PFMatrixList对象。更多搜索条件，请参阅帮助页getMatrixSet．
suppressMessages(library(JASPAR2014)) opts <- list() opts[["species"]] <- 9606 opts[["name"]] <- "RUNX1" opts[["type"]] <- "SELEX" opts[["all_versions"]] <- TRUE PFMatrixList <- getMatrixSet(JASPAR2014, opts) PFMatrixList #> PFMatrixList of length 1 #> names(1): MA0002.1 opts2 <- list() opts2[["type"]] <- "SELEX" PFMatrixList2 <- getMatrixSet(JASPAR2014, opts2) PFMatrixList2 #> PFMatrixList of length 111 #> names(111): MA0004.1 MA0006.1 MA0008.1 MA0009.1…Ma0588.1 ma0589.1 ma0590.13.2存储、删除和初始化JASPAR2014数据库
我们还提供了一些函数来初始化一个空的JASPAR2014风格的数据库，存储newPFMatrix或PFMatrixList或者根据ID删除一些记录。数据库的后端是SQLite．
db <- "myMatrixDb. "sqlite" initializeJASPARDB(db, version="2014") #> [1] "Success" data("MA0043") storeMatrix(db, MA0043) #> [1] "Success" deleteMatrixHavingID(db，"MA0043.1") file.remove(db) #> [1] TRUE4PFM, PWM和ICM方法
本节将介绍矩阵运算，包括从PFM到PWM和ICM的转换，轮廓矩阵比较，动态随机轮廓生成。
4．1PFM到PWM
该方法toPWM可以转换PFM到PWM(沃瑟曼和桑德林2004)．可选参数包括类型，pseudocounts，bg．这个包中的实现与Biostrings．
首先，toPWM允许输入矩阵有不同的列和，这意味着计数矩阵可以有不同数量的序列贡献到每一列。这种情况很少见，但存在于JASPAR SELEX数据中。
其次，我们可以指定定制pseudocounts．pseudocounts对于纠正少量计数或在对数转换前消除零值是必要的。在TFBS perl模块中，每列所包含的序列数的平方根。然而，它已被证明过于苛刻(Nishida, Frith, and Nakai 2009)．因此，使用默认值0.8。当然，它可以更改为其他自定义值，甚至可以为每列更改不同的值。
pwm < - toPWM(烤瓷,pseudocounts = 0.8) pwm # >的对象类PWMatrix # > ID: MA0004.1 # >名称:Arnt # >矩阵类:拉链# >链:+ # > pseudocounts: 0.8 # >标签:# > $家族# >[1]“Helix-Loop-Helix”# > # > $物种# >[1]“10090”# > # > $ tax_group # >[1]“脊椎动物”# > # > $ medline # >[1]“7592839”# > # > $类型# >[1]“SELEX”# > # > $ ACC # >[1]“P53762”# > # > $ pazar_tf_id # >[1]“TF0000003”# > # > $ TFBSshape_ID # >[1]“十一”# > # > $ TFencyclopedia_ID # >[1]“580”# > # >背景:> A -0.3081223 > A -4.700440 -4.700440 -4.700440 -4.700440 #> G -4.7004397 -2.115477 -4.700440 1.957772 -4.700440 1.957772 #> T -4.7004397 -4.700440 -4.700440 -4.700440 1.957772 #> T -4.7004397 -4.700440 -4.700440 -4.700440 1.957772 -4.700440 > T -4.7004397 -4.700440 -4.700440 -4.700440 1.957772 -4.7004404.2PFM到ICM
该方法toICM能否将PFM转换为ICM(Schneider et al. 1986)．除了相似的pseudocounts，bg，你也可以选择做施耐德修正。
信息内容矩阵的列和在0(没有基础偏好)和2(只使用了1个基础)之间。通常这些信息用于绘制序列日志。
PFM如何转换为ICM:我们有PFM矩阵\ \ (x)，基本背景频率\ \ (bg)，\ (pseudocounts \)为修正。
\[Z[j] = \sum_{i=1}^{4} x[i,j]\]
\ [p (i, j) = {(x (i, j) + bg \[我]倍pseudocounts [j]) \ / (Z [j] + \ sum_{我}bg \[我]倍pseudocounts [j]} \]
\ [D [j] = \ log_2 {4} + \ sum_ {i = 1} ^ {4} p (i, j) *日志{p (i, j)} \ \]
\[ICM[i,j] = p[i,j] \times D[j]\]
icm <- toICM(pfm, pseudocounts=0.8, schneider=TRUE) icm# >类ICMatrix的对象#> ID: MA0004.1 #>名称:Arnt #>矩阵类:zipp - type #> strand: + #> pseudocounts: 0.8 #> schneider修正:TRUE #>标签:# > $家族# >[1]“Helix-Loop-Helix”# > # > $物种# >[1]“10090”# > # > $ tax_group # >[1]“脊椎动物”# > # > $ medline # >[1]“7592839”# > # > $类型# >[1]“SELEX”# > # > $ ACC # >[1]“P53762”# > # > $ pazar_tf_id # >[1]“TF0000003”# > # > $ TFBSshape_ID # >[1]“十一”# > # > $ TFencyclopedia_ID # >[1]“580”# > # >背景:# > C G T # > 0.25 - 0.25 0.25 - 0.25 # >矩阵:#> [，1] [，2] [，3] [，4] [，4] [，5] [，6] #> a 0.203985791 1.30439694 0.01588159 0.01588159 0.01588159 0.01588159 0.01588159 #> c 0.786802337 0.01358747 1.60404024 0.01588159 0.01588159 1.60404024 #> t 0.009713609 0.01358747 0.015882481 0.01588159 1.60404024 #> t 0.009713609 0.01358747 0.01588159 0.01588159 1.60404024 0.01588159 0.01588159为了绘制序列标志，我们使用包seqlogo．在序列标识中，每个位置给出每个核苷酸获得的信息内容。核苷酸对应的字母越高，信息含量越大，在该位置获得该核苷酸的概率也就越大。
seqLogo (icm)
4.３将PFM对准自定义矩阵或IUPAC字符串
在某些情况下，评估现有轮廓矩阵(如JASPAR)与新发现矩阵的相似性是有益的(就像使用Genbank时使用BLAST进行序列数据比较一样)。
TFBSTools提供了使用改进的Needleman-Wunsch算法比较PFM对或带IUPAC字符串的PFM的工具(Sandelin et al. 2003)．
##一对一比较数据(MA0003.2)数据(MA0004.1) pfmSubject <- MA0003.2 pfmQuery <- MA0004.1 PFMSimilarity(pfmSubject, pfmQuery) #>评分relScore #> 7.294736 60.789466 #一到几个比较PFMSimilarity(pfmList, pfmQuery) #> $pfm1 #>评分relScore #> 12 100 #> #> $pfm2 #>评分relScore #> 12 100 # align IUPAC string IUPACString <- "ACGTMRWSYKVHDBN" PFMSimilarity(pfmList, pfmQuery) #>评分relScore #>，IUPACString) #> $pfm1 #> score relScore #> 8.81500 73.45833 #> $pfm2 #> score relScore #> 8.81500 73.458334.4PWM相似
在三次测量中测量两个PWM矩阵的相似性:归一化欧氏距离，皮尔森相关而且Kullback Leibler散度(Linhart, Halperin, and Shamir 2008)．给定两个概率型pwm，\ (P ^ 1 \)而且\ (P ^ 2 \),在那里\ (l \)就是长度。\ (P ^ j_{我,b} \)值在列吗我\ \ ()与基础\ (b \)在脉宽调制\ (j \)．归一化欧氏距离计算在
\ [D (P ^ 1, P ^ 2) = {1 \ / {\ l sqrt {2}}} \ cdot \ sum_ {i = 1} ^ {l} \√6 {\ sum_ {b \ {\ {A、C、G、T \}}} (P_{我,b} ^ 1-P_{我,b} ^ 2) ^ 2} \]
这个距离在0(完全相同)和1(完全不相似)之间。
中计算pearson相关系数
\ [r (P ^ 1, P ^ 2) = {1 \ / l} \ cdot \ sum_ {i = 1} l ^ {\ sum_ {b \ \ {A、C、G、T \}} (P_{我,b} ^ 1 - 0.25) (P_{我,b} ^ 2 - 0.25) \ / \√6 {\ sum_ {b \ \ {A、C、G、T \}} (P_{我,b} ^ 1 - 0.25) ^ 2 \ cdot \ sum_ {b \ \ {A、C、G、T \}} (P_{我,b} ^ 2 - 0.25) ^ 2}} \]
计算Kullback-Leibler散度
\ [KL (P ^ 1, P ^ 2) = {1 \ / {2 l}} \ cdot \ sum_ {i = 1} ^ l \ sum_ {b \ \ {A、C、G、T \}} (P_{我,b} ^ 1 \ log {P_{我,b} ^ 1 \ / P_{我,b} ^ 2} + P_{我,b} ^ 2 \ log {P_{我,b} ^ 2 \ / {P_{我,b} ^ 1}}) \]
data(MA0003.2) data(MA0004.1) pwm1 <- toPWM(MA0003.2, type="prob") pwm2 <- toPWM(MA0004.1, type="prob") PWMSimilarity(pwm1, pwm2, method="欧氏")#> [1]0.5134956 PWMSimilarity(pwm1, pwm2, method="Pearson") #> [1] 0.2828507 PWMSimilarity(pwm1, pwm2, method="KL") #> [1] 2.3858664．5动态随机剖面生成
在本节中，我们将演示随机轮廓矩阵生成矩阵排列和概率抽样的能力。在许多计算/模拟研究中，特别需要一组随机矩阵。一些案例包括估计TFBS与转录起始位点之间的距离，评估矩阵之间的比较(Bryne et al. 2008)．期望这些随机矩阵与所选配置文件具有相同的统计特性，例如核苷酸含量或信息含量。
排列法相对简单。它只是简单地打乱每个矩阵中受约束的列，或者是所有选定矩阵中的列。概率抽样比较复杂，可以分两步进行:
- 在JASPAR中对所有可用矩阵训练Dirichlet多项混合模型。
- 随机列从训练Dirichlet模型的后验分布基于选定的轮廓进行采样。
##矩阵排列permuteMatrix(pfmQuery) #>类PFMatrix对象#> ID: MA0004.1 #>名称:Arnt #>矩阵类:拉链类型#>链:+ #>标签:# > $评论# > [1 ] "-" #> #> $ 家庭# >[1]“Helix-Loop-Helix”# > # > $ medline # >[1]“7592839”# > # > $ pazar_tf_id # >[1]“TF0000003”# > # > $ tax_group # >[1]“脊椎动物”# > # > $ tfbs_shape_id # >[1]“十一”# > # > $ tfe_id # >[1]“580”# > # > $类型# >[1]“SELEX”# # > >收藏美元# >[1]“核心”# > # > $物种# >[1]“10090”# > # > $ acc # >[1]“P53762”# > # >背景:# > C G T # > 0.25 - 0.25 0.25 - 0.25 # >矩阵:#> [，1] [，2] [，3] [，4] [，5] [，6] #> A 0 19 0 4 0 0 0 #> C 0 0 0 0 0 16 0 20 #> G 20 1 0 0 0 20 0 #> T 0 0 20 0 0 0 0 permuteMatrix(pfmList, type="intra") #> PFMatrixList长度为2 #>名称(2):pfm1 pfm2 permuteMatrix(pfmList, type="inter") #> PFMatrixList长度为2 #>名称(2):pfm1 pfm2##狄利克雷模型训练数据(MA0003.2)数据(MA0004.1) pfmList <- PFMatrixList(pfm1=MA0003.2, pfm2=MA0004.1, use.names=TRUE) dmmParameters <- dmmEM(pfmList, K=6, alg="C")5TFFM方法
5.1TFFM的图形表示
基本pwm可以用上面所示的序列标识图形化地表示。考虑到二核苷酸依赖性，需要一种新的TFFM图形表示。
对于序列标识的上部，我们表示该位置的核苷酸概率\ (p \)每个可能的核苷酸的位置\ (p - 1 \)．因此，每一列表示TFBS中的一个位置，每一行表示在该位置发现的核苷酸概率。每一行都假设前一个隐藏状态释放了一个特定的核苷酸。对应位置的两列之间的交点\ (p \)和一行对应核苷酸\ (n \)给出每个核苷酸在相应位置的概率\ (p \)如果\ (n \)在位置上被看到了吗\ (p - 1 \)．表示序列标志的不透明度与TFFM可能使用的行的概率成比例。
##序列标志的一阶TFFM seqLogo(tffmFirst)
##序列标识详细的TFFM seqLogo(tffmDetail)
6扫描序列和脉宽调制模式的对齐
6.1searchSeq
searchSeq用PWM中表示的模式扫描核苷酸序列。参数strand控制搜索序列中的哪一条链。当它是_*_时，两个链都将被扫描。
一个SiteSet对象将返回，该对象可以导出为GFF3或GFF2格式。比赛得分的经验p值可以用精确的方法计算TFMPvalue或者抽样分数的分布。
library(Biostrings) data(MA0003.2) data(MA0004.1) pwmList <- PWMatrixList(MA0003.2=toPWM(MA0003.2)， MA0004.1=toPWM(MA0004.1)， use.names=TRUE) subject <- DNAString("GAATTCTCTCTTGTTGTAGTCTCTTGACAAAATG") siteset <- searchSeq(pwm, subject, seqname="seq1"， min.score="60%"， strand="*") sitesetList <- searchSeq(pwmList, subject, seqname="seq1"， min.score="60%"， min.score="60%"，strand="*") ##生成gff2或gff3风格的输出头(writeGFF3(siteset)) #> seqname源特征开始结束评分strand帧#> 1 seq1 TFBS TFBS 8 13 -1.888154 +。#> 2 seq1 TFBS TFBS 21 26 -1.888154 +。#> 3 seq1 TFBS TFBS 29 34 -3.908935 +。#> 4 seq1 TFBS TFBS 8 13 -1.961403 -。#> 5 seq1 TFBS TFBS 10 15 -3.908935 -。#> 6 seq1 TFBS TFBS 21 26 -1.961403 -。# b> 1 TF=Arnt;class= zipp - type;sequence= ctcttg# b> 2 TF=Arnt;class= zipp - type;sequence= aaaatg# > 4 TF=Arnt;class= zipp - type;sequence= caagag# > 5 TF=Arnt;class= zipp - type;sequence= aacaag# > 6 TF=Arnt;class= zipp - type;sequence=CAAGAG头(writeGFF3(sitesetList)) #> seqname源功能开始结束得分帧#> MA0003.2 seq1 TFBS TFBS 18 32 -16.437682 -。#> MA0004.1.1 seq1 TFBS TFBS 8 13 -1.888154 +。#> MA0004.1.2 seq1 TFBS TFBS 21 26 -1.888154 +。 #> MA0004.1.3 seq1 TFBS TFBS 29 34 -3.908935 + . #> MA0004.1.4 seq1 TFBS TFBS 8 13 -1.961403 - . #> MA0004.1.5 seq1 TFBS TFBS 10 15 -3.908935 - . #> attributes #> MA0003.2 TF=TFAP2A;class=Zipper-Type;sequence=TTTTGTCAAGAGACT #> MA0004.1.1 TF=Arnt;class=Zipper-Type;sequence=CTCTTG #> MA0004.1.2 TF=Arnt;class=Zipper-Type;sequence=CTCTTG #> MA0004.1.3 TF=Arnt;class=Zipper-Type;sequence=AAAATG #> MA0004.1.4 TF=Arnt;class=Zipper-Type;sequence=CAAGAG #> MA0004.1.5 TF=Arnt;class=Zipper-Type;sequence=AACAAG head(writeGFF2(siteset)) #> seqname source feature start end score strand frame #> 1 seq1 TFBS TFBS 8 13 -1.888154 + . #> 2 seq1 TFBS TFBS 21 26 -1.888154 + . #> 3 seq1 TFBS TFBS 29 34 -3.908935 + . #> 4 seq1 TFBS TFBS 8 13 -1.961403 - . #> 5 seq1 TFBS TFBS 10 15 -3.908935 - . #> 6 seq1 TFBS TFBS 21 26 -1.961403 - . #> attributes #> 1 TF "Arnt"; class "Zipper-Type"; sequence "CTCTTG" #> 2 TF "Arnt"; class "Zipper-Type"; sequence "CTCTTG" #> 3 TF "Arnt"; class "Zipper-Type"; sequence "AAAATG" #> 4 TF "Arnt"; class "Zipper-Type"; sequence "CAAGAG" #> 5 TF "Arnt"; class "Zipper-Type"; sequence "AACAAG" #> 6 TF "Arnt"; class "Zipper-Type"; sequence "CAAGAG" ## get the relative scores relScore(siteset) #> [1] 0.6652185 0.6652185 0.6141340 0.6633668 0.6141340 0.6633668 relScore(sitesetList) #> $MA0003.2 #> [1] 0.6196884 #> #> $MA0004.1 #> [1] 0.6652185 0.6652185 0.6141340 0.6633668 0.6141340 0.6633668 ## calculate the empirical p-values of the scores pvalues(siteset, type="TFMPvalue") #> [1] 0.02734375 0.02734375 0.04638672 0.04052734 0.04638672 0.04052734 pvalues(siteset, type="sampling") #> [1] 0.0234 0.0234 0.0544 0.0359 0.0544 0.03596.2searchAln
searchAln扫描由PWM表示的模式的成对对齐。它只报告那些出现在两个序列的等效位置并且在两个序列中都超过指定阈值的命中，并且在指定的对齐区域之上找到。
library(Biostrings) data(MA0003.2) pwm <- toPWM(MA0003.2) aln1 <- DNAString("ACTTCACCAGCTCCCTGGCGGTAAGTTGATC- AAAGG- aaacgcaaagttcaag ") aln2 <- DNAString(" gtttcactacttcctttcgggtaagtaaatatatataaaaaatatataatttcatc ") sitePairSet <- searchn (pwm, aln1, aln2, seqname1="seq1"， seqname2="seq2"， min.score="50%"， cutoff=0.5, strand="*"，type="any") ##生成gff风格输出头(writeGFF3(sitePairSet)) #> seqname源功能开始结束评分链帧#> 1 seq1 TFBS TFBS 6 20 -9.515444 +。#> 2 seq1 TFBS TFBS 7 21 -13.348617 +。#> 3 seq1 TFBS TFBS 8 22 -13.182322 +。#> 4 seq1 TFBS TFBS 9 23 -3.729917 +。#> 5 seq1 TFBS TFBS 10 24 -7.677850 +。#> 6 seq1 TFBS TFBS 14 28 -20.774619 +。# # >属性> 1 TF = TFAP2A;类=拉链;序列= ACCAGCTCCCTGGCG # > 2 TF = TFAP2A;类=拉链;序列= CCAGCTCCCTGGCGG # > 3 TF = TFAP2A;类=拉链;序列= CAGCTCCCTGGCGGT # > 4 TF = TFAP2A;类=拉链;序列= AGCTCCCTGGCGGTA # > 5 TF = TFAP2A;类=拉链;序列= GCTCCCTGGCGGTAA # > 6 TF = TFAP2A;类=拉链;序列= CCTGGCGGTAAGTTG头(writeGFF2 (sitePairSet)) # > seqname源特性开始结束分数链帧# > 1 seq1 TFBS TFBS 20 -9.515444 + 6。#> 2 seq1 TFBS TFBS 7 21 -13.348617 +。#> 3 seq1 TFBS TFBS 8 22 -13.182322 +。 #> 4 seq1 TFBS TFBS 9 23 -3.729917 + . #> 5 seq1 TFBS TFBS 10 24 -7.677850 + . #> 6 seq1 TFBS TFBS 14 28 -20.774619 + . #> attributes #> 1 TF "TFAP2A"; class "Zipper-Type"; sequence "ACCAGCTCCCTGGCG" #> 2 TF "TFAP2A"; class "Zipper-Type"; sequence "CCAGCTCCCTGGCGG" #> 3 TF "TFAP2A"; class "Zipper-Type"; sequence "CAGCTCCCTGGCGGT" #> 4 TF "TFAP2A"; class "Zipper-Type"; sequence "AGCTCCCTGGCGGTA" #> 5 TF "TFAP2A"; class "Zipper-Type"; sequence "GCTCCCTGGCGGTAA" #> 6 TF "TFAP2A"; class "Zipper-Type"; sequence "CCTGGCGGTAAGTTG" ## search the Axt alignment library(CNEr) axtFilesHg19DanRer7 <- file.path(system.file("extdata", package="TFBSTools"), "hg19.danRer7.net.axt") axtHg19DanRer7 <- readAxt(axtFilesHg19DanRer7) #> The number of axt files 1 #> The number of axt alignments is 133 sitePairSet <- searchAln(pwm, axtHg19DanRer7, min.score="80%", windowSize=51L, cutoff=0.7, strand="*", type="any", conservation=NULL, mc.cores=1) GRangesTFBS <- toGRangesList(sitePairSet, axtHg19DanRer7) GRangesTFBS$targetTFBS #> GRanges object with 33 ranges and 5 metadata columns: #> seqnames ranges strand | matrix.ID matrix.strand abs.score #> | #> [1] chr11 31128026-31128040 + | MA0003.2 + 2.08909 #> [2] chr11 31128028-31128042 + | MA0003.2 - 1.40496 #> [3] chr11 31437161-31437175 + | MA0003.2 + 7.78600 #> [4] chr11 31437163-31437177 + | MA0003.2 - 7.79663 #> [5] chr11 31499628-31499642 + | MA0003.2 + 2.23903 #> ... ... ... ... . ... ... ... #> [29] chr11 32198154-32198168 + | MA0003.2 + 2.51181 #> [30] chr11 32221916-32221930 + | MA0003.2 + 4.67226 #> [31] chr11 32221918-32221932 + | MA0003.2 - 10.60956 #> [32] chr11 32456353-32456367 + | MA0003.2 + 4.20454 #> [33] chr11 32456355-32456369 + | MA0003.2 - 1.17697 #> rel.score sitesSeq #> #> [1] 0.813035 GGCCCCCCGTGGCGT #> [2] 0.805896 CCCCCCGTGGCGTTT #> [3] 0.872489 ACTGGCCTAGGGACA #> [4] 0.872600 TGGCCTAGGGACAGG #> [5] 0.814600 TCTGACCCAGAGGCA #> ... ... ... #> [29] 0.817447 ACCGTCTTTAGGGGA #> [30] 0.839993 TCTACCCTCGAGCAC #> [31] 0.901956 TACCCTCGAGCACTG #> [32] 0.835112 AAGGGCCCGTAGCGA #> [33] 0.803516 GGGCCCGTAGCGACA #> ------- #> seqinfo: 3 sequences from an unspecified genome; no seqlengths GRangesTFBS$queryTFBS #> GRanges object with 33 ranges and 5 metadata columns: #> seqnames ranges strand | matrix.ID matrix.strand abs.score #> | #> [1] chr25 15057430-15057444 + | MA0003.2 + 2.65652 #> [2] chr25 15057432-15057446 + | MA0003.2 - 5.00568 #> [3] chr25 15149515-15149529 + | MA0003.2 + 5.68610 #> [4] chr25 15149517-15149531 + | MA0003.2 - 6.43485 #> [5] chr25 15163682-15163696 + | MA0003.2 + 1.49257 #> ... ... ... ... . ... ... ... #> [29] chr25 15359169-15359183 + | MA0003.2 + 3.13474 #> [30] chr25 15472022-15472036 + | MA0003.2 + 1.34057 #> [31] chr25 15472024-15472038 + | MA0003.2 - 5.67419 #> [32] chr25 15570786-15570800 + | MA0003.2 + 1.60462 #> [33] chr25 15570789-15570803 + | MA0003.2 - 5.26177 #> rel.score sitesSeq #> #> [1] 0.818957 GGCCCCCTGCGGCTC #> [2] 0.843473 CCCCCTGCGGCTCTC #> [3] 0.850574 TCTCGCCTAAGGAAA #> [4] 0.858388 TCGCCTAAGGAAAAA #> [5] 0.806810 TGCCACCCTTGGGCA #> ... ... ... #> [29] 0.823948 ACCATCTTTAGGGGA #> [30] 0.805224 TCTCCCCTCCAGCTC #> [31] 0.850450 TCCCCTCCAGCTCTG #> [32] 0.807979 CCGTACCTGCAGGCC #> [33] 0.846146 TACCTGCAGGCCCCT #> ------- #> seqinfo: 3 sequences from an unspecified genome; no seqlengths 6.3searchPairBSgenome
searchPairBSgenome是用来进行基因组系统发育足迹的。给定两个BSgenome一个从一个基因组转移到另一个基因组的链文件，searchPairBSgenome鉴定在两个基因组中都保守的转录因子结合位点。
library(rtracklayer) library(JASPAR2014) library(BSgenome.Hsapiens.UCSC.hg19) library(BSgenome.Mmusculus.UCSC.mm10) pfm <- getMatrixByID(JASPAR2014, ID="MA0004.1") pwm <- toPWM(pfm) chain <- import.chain("Downloads/hg19ToMm10.over.chain") sitePairSet <-hg19 BSgenome.Mmusculus.UCSC。mm10 chr1 =“chr1 chr2 =“chr1 min.score =“90%”,链= " + ",链=链)7使用新创Motif发现软件
在本节中，我们将介绍外部主题发现程序的包装器函数。到目前为止,MEME是支持的。# #MEMErunMEME需要一个DNAStringSet或一套字符作为输入，并返回MotifSet对象。
motifSet <- runMEME(file.path(system. path))file("extdata"， package="TFBSTools")， "crp0.s")， binary="meme"， arguments=list("-nmotifs"=3)) ##获取站点序列和周围序列sitesSeq(motifSet, type="all") ##只获取站点序列sitesSeq(motifSet, type="none") consensusMatrix(motifSet)8会话信息
这是的输出sessionInfo ()在编译本文件的系统上:
#> R version 4.2.0 RC (2022-04-19 r82224) #>平台:x86_64-pc-linux-gnu (64-bit) #>运行在:Ubuntu 20.04.4 LTS #> #>矩阵产品:默认#> BLAS: /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas。所以#> LAPACK: /home/biocbuild/bbs-3.15-bioc/R/lib/libRlapack。so #> #> locale: #> [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# > [3] LC_TIME=en_GB LC_COLLATE= c# > [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 #> [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# > [9] LC_ADDRESS=C LC_TELEPHONE= c# > [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C #> #>附加基础包:#> [1]stats4 stats graphics grDevices utils datasets methods #> [8] base #> #>其他附加包:# b> [1] CNEr_1.32.0 JASPAR2014_1.31.0 Biostrings_2.64.0 #> [4] GenomeInfoDb_1.32.0 XVector_0.36.0 IRanges_2.30.0 #> [7] S4Vectors_0.34.0 BiocGenerics_0.42.0 TFBSTools_1.34.0 #> [10] BiocStyle_2.24.0 #> #>通过命名空间加载(并且未附加):#> [1] bitops_1.0-7 matrixStats_0.62.0 #> [5] dirichletmultiomial_1.38.0 bit64_0.5 #> [5] httr_1.4.2 tools_4.2.0 #> b[7] bslib_0.3.1 utf8_1.2.2 #> [9] R6_2.5.1 seqLogo_1.62.0 #> [11] DBI_1.1.2 colorspace_2.0-3 #> [15] compiler_4.2.0 cli_3.3.0 #> [17] Biobase_2.56.0 delayedarray_1.2.0 #> [19] rtracklayer_1.6.0 bookdown_0.26 #> [21] sass_0.4.1 caTools_1.18.2 #> [25] string_1 .4.0 digest_0.6.29 #> [27] Rsamtools_2.12.0 #R.utils_2.11.0 # > [29] rmarkdown_2.14 pkgconfig_2.0.3 # > [31] htmltools_0.5.2 MatrixGenerics_1.8.0 # > [33] highr_0.9 BSgenome_1.64.0 # > [35] fastmap_1.1.0 rlang_1.0.2 # > [37] RSQLite_2.2.12 jquerylib_0.1.4 # > [39] BiocIO_1.6.0 generics_0.1.2 # > [41] jsonlite_1.8.0 BiocParallel_1.30.0 # > [43] gtools_3.9.2 R.oo_1.24.0 # > [45] dplyr_1.0.8 rcurl_1.98 - 1.6 # > [47] magrittr_2.0.3 GO.db_3.15.0 # > [49] GenomeInfoDbData_1.2.8 Matrix_1.4-1 # > [51] Rcpp_1.0.8.3 munsell_0.5.0 # > [53] fansi_1.0.3R.methodsS3_1.8.1 # > [55] lifecycle_1.0.1 stringi_1.7.6 # > [57] yaml_2.3.5 SummarizedExperiment_1.26.0 # > [59] zlibbioc_1.42.0 plyr_1.8.7 # > [61] grid_4.2.0 blob_1.2.3 # > [63] parallel_4.2.0 crayon_1.5.1 # > [65] lattice_0.20-45 annotate_1.74.0 # > [67] KEGGREST_1.36.0 hms_1.1.1 # > [69] magick_2.7.3 knitr_1.38 # > [71] pillar_1.7.0 GenomicRanges_1.48.0 # > [73] rjson_0.2.21 reshape2_1.4.4 # > [75] TFMPvalue_0.0.8 xml_3.99 - 0.9 # > [77] glue_1.6.2 evaluate_0.15 # > [79] BiocManager_1.30.17 png_0.1-7 # >[81] vctrs_0.4.1 tzdb_0.3.0 #> [83] poweRlaw_0.70.6 gtable_0.3.0 #> [85] purrr_0.3.4 assertthat_0.2.1 #> [87] cachem_1.0.6 ggplot2_3.3.5 #> [89] xfun_0.30 xtable_1.8-4 #> [91] pracma_2.3.8 restfulr_0.0.13 #> [93] tibble_3.1.6 AnnotationDbi_1.58.0 #> [95] GenomicAlignments_1.32.0 memoise_2.0.1 #> [97] ellipsis_0.3.2参考文献
贝利，T L, C埃尔坎。1994。“用期望最大化拟合混合模型以发现生物聚合物中的基序。”Proc Int Conf Intell Syst Mol Biol2:几个。
Bryne, Jan Christian, Eivind Valen, Man-Hung Eric Tang, Troels Marstrand, Ole Winther, Isabelle da Piedade, Anders Krogh, Boris Lenhard和Albin Sandelin, 2008。“JASPAR，转录因子结合谱的开放获取数据库:2008年更新中的新内容和工具。”核酸测定。36(数据库问题):D102-106。https://doi.org/10.1093/nar/gkm955．
鲍里斯·伦哈德和怀斯·W·沃瑟曼，2002年。TFBS:转录因子结合位点分析的计算框架生物信息学18(8): 1135-6。
林哈特，查姆，尤尼特·哈尔佩林和罗恩·沙米尔，2008年。“转录因子和microRNA Motif发现:Amadeus平台和后生动物目标集纲要。”基因组Res。18(7): 1180-9。https://doi.org/10.1101/gr.076117.108．
麦瑟利尔，安东尼，怀斯·W·沃瑟曼，2013。“下一代转录因子结合位点预测”公共科学图书馆第一版。医学杂志。9 (9): e1003214。https://doi.org/10.1371/journal.pcbi.1003214．
Nishida, Keishin, Martin C Frith和Kenta Nakai, 2009。“转录因子结合位点的伪计数。”核酸测定。37(3): 939-44。https://doi.org/10.1093/nar/gkn1019．
桑德林，阿尔宾，安妮特Höglund，鲍里斯·伦哈德和怀斯·W·沃瑟曼。2003。生物连接基因簇的酵母调控序列的综合分析功能。中国。基因组学3(3): 125-34。https://doi.org/10.1007/s10142-003-0086-6．
施耐德，T D, G D Stormo, L Gold和A Ehrenfeucht, 1986。核苷酸序列上结合位点的信息含量J. Mol.生物学。188(3): 415-31。
Wasserman, Wyeth W, Albin Sandelin, 2004。应用生物信息学鉴定调控元件自然出版集团5(4): 276-87。