2.1安装

运行bioassayR的R软件可从CRAN (http://cran.at.r-project.org/)．的bioassayR包可以从R中使用BiocManager包中。

如果(!requireNamespace("BiocManager"， quiet =TRUE)) install.packages("BiocManager") BiocManager::install("bioassayR") #安装包

2．2加载包和文档

library(bioassayR) #加载包

vignette(“bioassayR”)#从R打开本手册

2．3快速教程

本示例将引导您创建一个新的空数据库，添加示例小分子生物活性数据，并对这些数据执行查询。如果您只对查询预构建的PubChem BioAssay数据库感兴趣，请跳到后面的“预构建数据库示例:调查已知药物的活性”一节。

这个例子包括一个抗生素发现实验的真实实验数据。这些数据是一个“验证性生物测定”，其中57个小分子被筛选对抗来自甲羟戊酸激酶蛋白链球菌引起的肺炎(SP)的细菌。甲戊酸激酶抑制剂是一种可能对这种臭名昭著的细菌有效的抗生素药物。这些数据由Thomas S. Leyh博士作为检测标识符(辅助)1000发表在NCBI PubChem BioAssay数据库中。

首先，创建一个新数据库。本手册使用临时文件，但您可以替换tempfile如果希望为以后保存结果数据库，则使用您选择的文件名调用函数。

library(bioassayR) myDatabaseFilename <- tempfile() mydb <- newBioassayDB(myDatabaseFilename, indexed=FALSE)

接下来，指定计划加载的数据的来源和版本。这是一个必需的步骤，可以更容易地跟踪数据的起源。如果数据没有版本控制，可以随意使用当前日期作为版本号。

addDataSource(mydb, description="PubChem BioAssay"， version="unknown")

添加数据源后，可以创建或导入data.frame其中包含您的实验中每个分子的活性分数。这data.frame必须包含三列，其中包括每个化合物的cid(唯一化合物标识符)，二进制活动分数(1=活跃，0=不活跃，NA=不确定)和数字活动分数。查阅生物测定有关格式化的详细信息，请参见手册页data.frame．的bioassayR包中包含一个活动分数数据帧示例，可以如下方式访问:

Data (samplebioassay) samplebioassay[1:10，] #打印前10个分数

## cid活动评分## 1 730195 0 0 ## 2 16749973 1 80 ## 3 16749974 1 80 ## 4 16749975 1 80 ## 5 16749976 1 80 ## 6 16749977 1 80 ## 7 16749978 1 80 ## 8 16749979 1 80 ## 9 16749980 1 80 ## 10 16749981 1 80

所有生物活性数据都被加载到数据库中，或从数据库中检索生物测定对象，其中包含有关分析实验设计、分子靶标和活性分数的详细信息。可以按照如下方式创建包含活动分数的生物测定对象。源id值必须与先前加载的值完全匹配addDataSource．用于检测的分子靶标是可选的，并且可以为单个检测指定无限数量的靶标作为传递给靶标选项的载体。目标类型字段应该是一个长度相等的向量，以相同的顺序描述每个目标的类型。

myAssay <- new("bioassay"，aid="1000"， source_id="PubChem bioassay"， assay_type="验证性"，有机体="未知"，scoring="活性等级"，targets="116516899"， target_types="蛋白"，scores=samplebioassay) myAssay

##类:生物测定##援助:1000 ## source_id: PubChem bioassay ## assay_type:验证##生物:未知##评分:活动排名##目标:116516899 ## target_types:蛋白质##总得分:57

的生物测定对象可以加载到数据库中loadBioassay函数。通过对不同的数据重复此步骤，可以将大量不同的分析结果加载到数据库中。

loadBioassay (mydb myAssay)

建议使用NCBI GI编号作为任何生物分子分析靶标的标签，因为这是在包中提供的预构建数据库中使用的。在某些情况下，还可以存储标识符翻译，其中包含来自其他生物数据库的相应id。例如，在本例中，甲戊酸激酶靶蛋白GI 116516899在UniProt数据库的Q8DR51中有一个对应的标识符。可以按如下方式存储以供将来参考。通过这种方式可以存储任意类别的任意数量的翻译。它还可以用于存储目标的注释数据。例如，如果目标上有序列级集群数据，可以将它们存储在“sequencecluclusters”类别中，并将集群号存储为标识符。

loadIdMapping(mydb, target="116516899"， category="UniProt"， identifier="Q8DR51")

等一下!当我们知道这实际上是对一种蛋白质的筛选时，我们意外地将该检测标记为“未知生物”链球菌引起的肺炎．将分析加载到数据库后，稍后可以使用getAssay函数。通过将此功能与删除分析的能力相结合dropBioassay函数)可以通过以下方式编辑数据库:(1)从数据库中取出化验结果，(2)从数据库中删除化验结果，(3)修改取出的对象，以及(4)重新加载化验结果。例如，我们可以为我们的分析更新物种注释如下:

tempAssay <- getAssay(mydb， "1000") # get assay from database dropBioassay(mydb， "1000") # delete assay from database organism(tempAssay) <- "Streptococcus pneumonia" # update organism loadBioassay(mydb, tempAssay)

建议在加载所有数据后对数据库进行索引。这大大加快了对数据库的访问速度，但如果在加载之前执行索引，也会降低数据的加载速度。

addBioassayIndex (mydb)

##创建索引:注意，对于大型数据库，这可能需要很长时间

建立索引后，可以查询数据库。下面是一些查询示例。首先查看由提供的数据库摘要bioassayR：

mydb

##类:BioassayDB ## assays: 1 ##来源:PubChem BioAssay ##源版本:未知##可写:是

接下来，您可以通过cid查询数据库中给定化合物的活动目标。在这种情况下，由于只加载了一个试验，只能找到一个靶标。实验加载更多的分析更有趣的结果!当使用预构建的PubChem BioAssay数据库时，这些目标作为NCBI蛋白标识符返回。

activeTargets (mydb, 16749979)

## fraction_active total_screens ## 116516899 11

如果目标翻译是在前面的步骤中加载的，则可以使用translateTargetId函数如下。这只接受一个目标，并将返回所选类别中所有对应标识符的向量。在某些情况下，您可能希望对该结果进行子集化，以便在数据库包含大量针对某些目标的标识符时仅获得单个标识符。

在这里，我们请求GI 116516899的UniProt标识符，就像前面存储的一样loadIdMapping函数。

translateTargetId(mydb, target="116516899"， category="UniProt")

##[1]“q8dr51”

最后，在分析后断开与数据库的连接可以减少数据损坏的机会。如果您正在使用只读模式下的预构建数据库(如预构建数据库示例部分所示)，则可以选择跳过此步骤，因为只有可写数据库容易因断开连接失败而损坏。

disconnectBioassayDB (mydb)

3．1加载用户提供的PubChem生物测定数据

本节演示从用户提供的数据创建新的生物活性数据库的过程。作为一个例子，我们将演示从NCBI PubChem BioAssay生物活性数据存储库下载一个分析，并将其加载到一个新的数据库的过程(Wang et al. 2012)．

首先，从PubChem BioAssay获取感兴趣的实验的两个文件:一个XML文件，其中包含如何执行实验的细节，一个CSV(逗号分隔值)文件，其中包含实际的活动分数。为了本例的目的，我们将使用来自试验1000的数据，这是针对甲羟戊酸激酶蛋白的57个小分子的验证性试验(滴定试验)。在“快速教程”中提供了关于此分析的更多详细信息，其中使用了相同的数据。这些文件可以从PubChem BioAssay下载http://pubchem.ncbi.nlm.nih.gov/或者从这个包中包含的示例数据存储库中加载，如下所示:

library(bioassayR) extdata_dir <- system. bat。file("extdata"， package="bioassayR") assayDescriptionFile <- file。path(extdata_dir， "exampleAssay.xml") activityScoresFile <- file. path(extdata_dir， "exampleAssay.xml")路径(extdata_dir,“exampleScores.csv”)

接下来，创建一个新的空数据库，用于加载这些数据。本例使用Rtempfile函数在随机位置创建数据库。如果希望保留生成的数据库，请替换myDatabaseFilename使用所需的路径和文件名。

myDatabaseFilename <- tempfile() mydb <- newBioassayDB(myDatabaseFilename, indexed=F)

我们还将向该数据库添加一个数据源，指定这里的数据反映PubChem BioAssay提供的分析。

addDataSource(mydb, description="PubChem BioAssay"， version="unknown")

PubChem BioAssay提供的XML文件包含了大量关于如何执行检测、分子目标和结果评分方法的详细信息。可以使用parsePubChemBioassay函数如下。的parsePubChemBioassay函数还需要一个. CSV文件，其中包含每个实验的活性分数，以PubChem BioAssay提供的标准CSV格式。对于来自PubChem BioAssay以外来源的数据，您可能需要编写自己的代码来解析分析细节—或者手动输入它们。

myAssay <- parsePubChemBioassay("1000"， activityScoresFile, assayDescriptionFile

##类:生物测定##援助:1000 ## source_id: PubChem bioassay ## assay_type:验证##生物:NA ##评分:IC50 ##目标:116516899 ## target_types:蛋白质##总评分:57

接下来，将从XML和CSV文件解析的结果数据加载到数据库中。这将在数据库中为化验本身和分子靶标创建记录。

loadBioassay (mydb myAssay)

要加载额外的分析，重复上述步骤。加载完所有数据后，可以通过向数据库添加索引来显著提高后续查询性能。

addBioassayIndex (mydb)

##创建索引:注意，对于大型数据库，这可能需要很长时间

在建立索引之后，对数据库执行测试查询，以确认正确加载了数据。

activeAgainst (mydb,“116516899”)

# # fraction_active total_assays # # 16749973 1 1 # # 16749973 # # 16749975 1 1 # # 16749976 1 1 # # 16749976 # # 16749978 1 1 1 # # 16749979 1 1 # # 16749979 # # 16749981 1 1 # # 16749981 1 1 # # 16749983 1 1 # # 16749983 # # 16749985 1 1 # # 16749985 1 1 # # 16749987 1 1 # # 16749987 # # 16749989 1 1 # # 16749989 1 1 # # 16749991 1 1 # # 16749991 # # 16749993 1 1 # # 16749993 1 1 # # 16749995 1 1 # # 16749995 # # 16749997 1 1 # # 16749997 1 1 # # 16749999 1 1 # # 16749999 # # 1 # # 16750001 1 116750002 1 1 ## 16750003 1 1 ## 16750004 1 1 ## 16750005 1 1 ## 16750006 1 1 ## 16750007 1 1 ## 16750016 1 1 1

最后，断开与数据库的连接以防止数据损坏。

disconnectBioassayDB (mydb)

3.2预建数据库示例:调查已知药物的活性

一个预先构建的数据库，包含大量公共领域的生物活性数据，来自PubChem BioAssay数据库，可以从http://chemmine.ucr.edu/bioassayr．虽然建议下载完整的数据库，但可以使用数据库的一个小子集运行此示例bioassayR包用于测试目的。的实用功能bioassayR用于鉴定类药物小分子的生物活性模式。在本例中，我们观察药物乙酰水杨酸(又名阿司匹林)的结合活性模式，并将这些结合数据与文献中注释的靶标进行比较DrugBank药品数据库(Wishart et al. 2008)．

药物库数据库是一个有价值的资源，包含大量关于人类药物活性的数据，包括已知的分子靶点。在本练习中，首先通过搜索此名称来查看注释的乙酰水杨酸分子靶标http://drugbank.ca．这将为我们在预构建的PubChem BioAssay数据库中找到的生物活性数据进行比较提供一个参考点。请注意，DrugBank还包含该化合物的PubChem CID，您可以使用它来查询bioassayR PubChem BioAssay数据库。

要开始，首先连接到数据库。的变量sampleDatabasePath可以用您下载的完整测试数据库的文件名替换，如果您想使用该文件名而不是此软件包中包含的小型示例版本。

library(bioassayR) extdata_dir <- system. bat。sampleDatabasePath <- file. file("extdata"， package="bioassayR")path(extdata_dir， "sampleDatabase.sqlite") pubChemDatabase <- connectBioassayDB(sampleDatabasePath)

接下来，使用activeTargets函数用于查找加载的数据库中乙酰水杨酸显示活性的所有蛋白质目标。这些目标id是由PubChem BioAssay提供的NCBI蛋白标识符(Tatusova et al. 2014)．在哪些情况下，这些结果与DrugBank中标注的目标一致或不一致?

drugTargets <- activeTargets(pubChemDatabase， "2244"

## fraction_active total_screens ## 116241312 1.0 1 ## 125987835 1.0 1 ## 166897622 1.0 1 ## 226694183 1.0 1 ## 317373262 1.0 6 ## 3914304 1.0 14 ## 3915797 0.8 10 ## 6686268 1.0 1 ## 754286265 1.0 22 ## 84028191 1.0 1

现在我们想连接到UniProt (Universal Protein Resource)数据库，以获取这些目标的注释细节
(贝特曼等，2015)．的biomaRtBioconductor软件包是一个很好的工具，但最好的工作UniProt标识符，而不是我们目前拥有的NCBI蛋白质标识符，因此我们必须首先翻译标识符(Durinck et al. 2009, 2005)．

我们将使用translateTargetId函数bioassayR获取每个NCBI蛋白标识符(GI)对应的UniProt标识符。这些标识符翻译是从UniProt获得的，并预先加载到数据库中。的translateTargetId只接受一个查询，并返回一个或多个UniProt标识符。这里我们称它为酸式焦磷酸钠哪个自动调用函数多次，存储的每个蛋白质一次drugTargets．在许多情况下，单个查询GI将转换为多个UniProt标识符。在本例中，我们只保留第一个，因为我们在这里寻找的注释细节可能对所有这些注释都是相同的。

#在每个目标标识符上运行translateTargetId uniProtIds <- lapply(row.names(drugTargets)， translateTargetId, database=pubChemDatabase, category="UniProt") #如果任何目标有多个UniProt ID，只保留第一个uniProtIds <- sapply(uniProtIds, function(x) x[1])

接下来，我们通过biomaRt连接到Ensembl以获得每个目标的描述智人基因。有关更多信息，请参阅biomaRt文档。检索这些数据后，我们调用匹配函数，以确保它们与查询数据的顺序相同。

library(biomaRt) ensembl <- useEnsembl(biomaRt ="ensembl"，dataset=" hapiens_gene_ensembl ") proteinDetails <- getBM(attributes=c("description"，"uniprotswissprot"，"external_gene_name")， filters=c("uniprotswissprot")， mart=ensembl, values=uniProtIds) proteinDetails <- proteinDetails[match(uniProtIds, proteinDetails$uniprotswissprot)，]

现在我们可以查看注释数据了。NAs表示蛋白质没有发现在智人合群，可能来自其他种。

proteinDetails

##描述## 4细胞色素P450，家族3，亚家族A，多肽4[来源:HGNC符号;Acc:HGNC:2637] ## 2细胞色素P450，家族1，亚家族A，多肽2[来源:HGNC符号;Acc:HGNC:2596] ## 1整合素，β 3(血小板糖蛋白IIIa，抗原CD61)[来源:HGNC符号;Acc:HGNC:6156] ## NA  ## 3整合素，α 2b (IIb/IIIa复合体的血小板糖蛋白IIb，抗原CD41)[来源:HGNC符号;Acc:HGNC:6138] ## 7前列腺素过氧化物合酶1(前列腺素G/H合酶和环加氧酶)[来源:HGNC符号;Acc:HGNC:9604] ## NA.1  ## 8前列腺素过氧化物合酶2(前列腺素G/H合酶和环加氧酶)[来源:HGNC符号;Acc:HGNC:9605] ## 6细胞色素P450，家族2，亚家族C，多肽9[来源:HGNC符号;Acc:HGNC:2623] ## NA.2  ## 5细胞色素P450，家族2，亚家族D，polypeptide 6[来源:HGNC Symbol;Acc:HGNC:2625] ## uniprot_swissprot external_gene_name ## 4 P08684 CYP3A4 ## 2 P05177 CYP1A2 ## 1 P05106 ITGB3 ## NA   ## 3 P08514 ITGA2B ## 7 P23219 PTGS1 ## NA.1   ## 8 P35354 PTGS2 ## 6 P11712 CYP2C9 ## NA.2   ## 5 P10635 CYP2D6

最后，让我们再次查看活动目标列表，以及旁边的注释。注意，这些只在长度和顺序上匹配，因为在上面的代码中，我们删除了每个目标蛋白的所有UniProt ID，只保留一个，然后用匹配函数以正确的顺序获取它们。

drugTargets <- drugTargets[! !na(proteinDetails[， 1])，] proteinDetails <- proteinDetails[!na(proteinDetails[， 1])，] cbind(proteinDetails, drugTargets)

##描述## 4细胞色素P450，家族3，亚家族A，多肽4[来源:HGNC符号;Acc:HGNC:2637] ## 2细胞色素P450，家族1，亚家族A，多肽2[来源:HGNC符号;Acc:HGNC:2596] ## 1整合素，β 3(血小板糖蛋白IIIa，抗原CD61)[来源:HGNC符号;Acc:HGNC:6156] ## 3整合素，α 2b (IIb/IIIa复合体的血小板糖蛋白IIb，抗原CD41)[来源:HGNC符号;Acc:HGNC:6138] ## 7前列腺素过氧化物合酶1(前列腺素G/H合酶和环加氧酶)[来源:HGNC符号;Acc:HGNC:9604] ## 8前列腺素过氧化物合酶2(前列腺素G/H合酶和环加氧酶)[来源:HGNC符号;Acc:HGNC:9605] ## 6细胞色素P450家族2亚家族C多肽9[来源:HGNC符号;Acc:HGNC:2623] ## 5细胞色素P450家族2亚家族Dpolypeptide 6[来源:HGNC Symbol;Acc:HGNC:2625] ## uniprot_swissprot external_gene_name fraction_active total_screens ## 4 P08684 CYP3A4 1.0 2 ## 2 P05177 CYP1A2 1.0 1 ## 1 P05106 ITGB3 1.0 1 ## 3 P08514 ITGA2B 1.0 1 ## 7 P23219 PTGS1 1.0 6 ## 8 P35354 PTGS2 0.8 10 ## 6 P11712 CYP2C9 1.0 1 ## 5 P10635 CYP2D6 1.0 1 ## 1

3．3确定目标选择性化合物

在前面的例子中，发现乙酰水杨酸对许多蛋白质具有结合活性，包括COX-1环加氧酶(NCBI Protein ID 166897622)。理论上COX-1活性是该分子作为非甾体抗炎药(NSAID)发挥作用的关键机制。在本例中，我们将寻找选择性结合COX-1的其他小分子，假设这些小分子可能值得作为潜在的非甾体抗炎药物进行进一步研究。这个例子展示了bioassayR可用于识别选择性结合到感兴趣的目标的小分子，并协助发现小分子药物和分子探针。

首先，我们将从连接到数据库开始。还是变量sampleDatabasePath可以替换为完整的PubChem BioAssay数据库的文件名(可从http://chemmine.ucr.edu/bioassayr)，如果你想使用该版本，而不是此软件包附带的小型示例版本。

library(bioassayR) extdata_dir <- system. bat。sampleDatabasePath <- file. file("extdata"， package="bioassayR")path(extdata_dir， "sampleDatabase.sqlite") pubChemDatabase <- connectBioassayDB(sampleDatabasePath)

的activeAgainst函数可用于显示数据库中对COX-1具有活性的所有小分子，如下所示。每一行代表一个小分子酸。标记为“总化验”的列显示了针对感兴趣的目标筛选每个小分子的总次数。标记为“活跃分数”的列显示了这些分数中标注为活跃的0到1之间的数字的部分。该功能允许用户将来自不同来源的不同结合分析作为重复，以帮助区分潜在的假结合结果与证明可重复性的结合结果。

activeCompounds <- activeAgainst(pubChemDatabase， "166897622") activeCompounds[1:10，] #查看前10个化合物

## fraction_active total_assays ## 237 1 1 ## 2244 1 1 ## 2662 1 4 ## 3033 1 1 ## 3194 1 1 ## 3672 13 ## 3715 1 6 ## 133021 1 2 ## 156391 1 1 ## 247704 1 1 1

仅观察与感兴趣靶标结合的化合物不足以识别候选药物，因为这些化合物中的一部分可能是无选择性靶标化合物(TUCs)，它们不加区别地与大量不同的蛋白质靶标结合。R函数selectiveAgainst向用户提供化合物列表，这些化合物对感兴趣的靶标显示活性(至少在一种测定中)，同时对其他靶标显示有限活性。

的maxCompounds选项限制返回结果的最大数目，而minimumTargets选择限制返回的化合物是针对指定的最小不同目标进行筛选的化合物。结果被格式化为data.frame其中每个行名代表一个不同的复合。第一列显示该化合物对不同目标显示活性的数量，第二列显示它对筛选目标的总数。

selectiveCompounds <- selectiveAgainst(pubChemDatabase， "166897622"， maxCompounds = 10, minimumTargets = 1) selectiveCompounds

## active_targets tested_targets ## 2662 11 ## 3033 11 ## 133021 11 ## 11314954 11 ## 13015959 11 ## 44563999 11 ## 44564000 11 ## 44564002 1 ## 44564003 11 1

在示例数据库中，这些化合物只显示了一个测试目标，因为加载的测试很少。我们鼓励用户使用PubChem BioAssay完整的数据库来尝试这个例子http://chemmine.ucr.edu/bioassayr为了获得更有趣和有信息的结果。

用户可以组合bioassayR与ChemmineR库，以获得这些目标选择性化合物的结构信息，然后执行进一步的分析-如结构聚类，可视化，和计算物理化学性质。

的ChemmineR软件库可以用来下载任何这些化合物的结构数据，并将这些结构可视化如下(Cao et al. 2008)．这个示例需要一个活动的internet连接，因为复合结构是从远程服务器获得的。

库(ChemmineR)结构<- getIds(as.numeric(row.names(selectiveCompounds)))

在这里，我们只看到了前四种化合物selectiveAgainst．参考附带的插图ChemmineR为了进一步可视化和分析这些结构的许多例子。

plot(structures[1:4]， print=FALSE) #将结构绘制到R图形设备

3.4按活性谱聚类化合物

本例演示了在几个不同的生物测定实验中通过相似的生物活性谱聚类小分子的示例。在许多情况下，在数据库中聚集所有的复合数据是cpu和内存密集型的，所以我们首先从数据库中提取这些数据的一个子集bioassaySet对象，然后将其转换为化合物与目标的活动矩阵，以便根据活动剖面的相似性进行后续聚类。这个函数getBioassaySetByCids提取给定化合物列表的活性数据。或者，可以使用该函数提取给定化验id列表的全部数据getAssays．

首先，连接到包含的示例数据库:

library(bioassayR) extdata_dir <- system. bat。sampleDatabasePath <- file. file("extdata"， package="bioassayR")path(extdata_dir， "sampleDatabase.sqlite") sampleDB <- connectBioassayDB(sampleDatabasePath)

接下来，从3个化合物中选择数据提取到bioassaySet对象进行后续分析。

compoundsOfInterest <- c("2244"， "2662"， "3715") selectedAssayData <- getBioassaySetByCids(sampleDB, compoundsOfInterest) selectedAssayData . c("2244"， "2662"， "3715"

##类:bioassaySet ##分析:792 ##化合物:3 ##目标:551 ##来源:PubChem BioAssay

这个函数perTargetMatrix将先前提取的活动数据转换为目标(行)与复合(列)的矩阵。通过用户指定的方式将击中同一目标的多个测定的数据汇总为单个值。如果你选择summarizeReplicates选项activesFirst，任何活跃的分数优先于不活跃的分数。如果你选择了这个选项模式最丰富的活性或非活性都存储在结果矩阵中。您还可以传递一个自定义函数来决定如何总结复制。在结果矩阵中，“2”表示活性化合物与目标组合，“1”表示非活性组合，“0”表示未经测试或不确定的组合。通过传递Inactive = FALSE选项，可以有选择地排除非活动结果。这里使用稀疏矩阵(忽略实际存储0值)来节省内存。在稀疏矩阵中，一个周期”。"项等于零的" 0 "项，但暗示不占用额外的内存空间。

在这里，我们创建了一个活动矩阵，选择包括非活动值，并根据统计模式总结复制:

myActivityMatrix <- perTargetMatrix(selectedAssayData, inactives=TRUE, summarizereplcopies = "mode") myActivityMatrix[1:15，] #打印前15行

## 15 x 3稀疏矩阵类“dgCMatrix”## 2244 2662 3715 ## 166897622 22 ## 125987835 2 .。2 . .2 . .2 . .2 . .2 . .2 . .2 . .## 1162413121 . . ## 160794 1 . . ## 21464101 1 . . ## 154146191 1 . . ## 68565074 1 . .

活动矩阵也可以选择使用原始数字分数或缩放和居中的z分数来创建，而不是离散的活动/不活动值。选项的手册页，以获得有关此选项的详细信息perTargetMatrix而且scaleBioassaySet功能。

接下来，我们将重新创建活性矩阵，其中在序列水平上非常相似的蛋白质目标(例如来自不同物种的同源物)被视为复制，并合并。

这个函数perTargetMatrix包含一个选项assayTargets这样你就可以为每个实验指定目标而不是从bioassaySet对象。这个函数assaySetTargets中每个化验的目标向量bioassaySet对象，其中每个元素的名称对应于它的分析标识符(aid)。这是必须传递给的格式perTargetMatrix为了指定哪些检测被视为要合并的副本，因此首先我们可以获得这些数据，然后用以相同方式格式化的自定义合并条件替换它们。

myAssayTargets <- assaySetTargets(selectedAssayData) myAssayTargets[1:5] #打印前5个目标

## 399401 399404 399411 443726 347221 ##“166897622”“166897622”“166897622”“166897622”“166897622”“166897622”

预构建的PubChem BioAssay数据库包括用kClust工具生成的序列级蛋白质目标聚类结果(选项s=2.93, E-value < 1*10^-4, c=0.8)(Hauser, Mayer, and Söding 2013)．每个集群都有一个唯一的编号，聚集在一起的目标被分配相同的集群编号。这些聚类结果存储在数据库中，作为类别“kClust”下的目标翻译。现在我们将访问这些翻译，并制作如下所示的复合与目标聚类矩阵。

#获取单个目标的kClust蛋白簇号translateTargetId(database = sampleDB, target = "166897622"， category = "kClust")

“1945”

- sapply(myAssayTargets, translateTargetId, database = sampleDB, category = "kClust") customMerge[1:5]

## 399401 399404 399411 443726 347221 ##“1945”“1945”“1945”“1945”“1945”“1945”“1945”

mergedActivityMatrix <- perTargetMatrix(selectedAssayData, inactives=TRUE, assayTargets=customMerge) mergedActivityMatrix[1:15，] #打印前15行

## 15 x 3稀疏矩阵类“dgCMatrix”## 2244 2662 3715 ## 1945 22 2 ## 1210 2 .。2 . .2 . .2 . .2 . .1 .;1 . .1 . .1 . .1 . . ## 5734 1 . . ## 1250 1 . . ## 3182 1 . . ## 5745 1 . .

请注意，合并后的矩阵更小，因为几个相似的蛋白质目标已经被分解成单个簇。

获取未合并矩阵的行数和列数(myActivityMatrix)

## [1] 418 3

#获取合并矩阵的行数和列数

## [1] 365

现在，为了能够使用二进制聚类方法，我们将矩阵简化为二进制矩阵，其中“1”表示活动，“0”表示非活动或未经测试的组合。我们提醒用户仔细考虑这在正在分析的具体实验中是否有意义。

binaryMatrix <- 1*(mergedActivityMatrix > 1) binaryMatrix[1:15，] #打印前15行

## 15 x 3稀疏矩阵类“dgCMatrix”## 2244 2662 3715 ## 1945 1 1 1 ## 1210 1 .。1 . .1 . .1 . .1 . .## 2260 0 . .0 . .## 1218 0 . .## 1484 0 . .## 4946 0 . . ## 5734 0 . . ## 1250 0 . . ## 3182 0 . . ## 5745 0 . .

如前所述，“0”和“。“稀疏矩阵中的项在数值上是等价的。

使用R内内置的欧几里得聚类函数进行聚类。这提供了一个树状图，根据它们的活性谱表明化合物之间的相似性。

转置矩阵<- t(binaryMatrix) distanceMatrix <- dist(转置矩阵)clusterResults <- hclust(distanceMatrix, method="average") plot(clusterResults)

比较活动概要和集群数据的第二种方法是将活动矩阵传递给ChemmineR化学信息学图书馆FPset(二进制指纹)对象。这将生物活性数据表示为二进制指纹(生物活性指纹)，这是每个化合物的二进制字符串，其中每个位表示对这些化合物显示活性的全部目标的活性(带1)或不活性(带0)。这样就可以对化合物之间的生物活性进行两两比较。看到ChemmineR文档在//www.andersvercelli.com/packages/ChemmineR/有关其他示例和详细信息。

library(ChemmineR) fpset <- bioactivityFingerprint(selectedAssayData) fpset

##一个具有3个分子的551位“生物活性”类型的“FPset”实例

通过将第一个化合物与所有化合物进行比较，使用FPset对象执行活动剖面相似性搜索。

fpSim(fpset[1]， fpset, method="Tanimoto")

## 2244 2662 3715 ## 1.0000000 0.1818182 0.1818182

计算所有对所有这些化合物的生物活性指纹相似矩阵。

simMA <- sapply(cid(fpset)， function(x) fpSim(fpset[x]， fpset, sorted=FALSE, method="Tanimoto")) simMA

## 2244 2662 3715 ## 2244 1.0000000 0.1818182 0.1818182 ## 2662 0.1818182 1.0000000 1.0000000 ## 3715 0.1818182 1.0000000 1.0000000

将相似度矩阵转换为距离矩阵并进行层次聚类。

clusterResults <- hclust(as.dist(1-simMA)， method="single")

一种可视化大量化合物的相对生物活性相似性的方法是使用多维标度图(又名MDS或主坐标分析)，其中每个化合物表示为一个点，生物活性距离与任意两点之间的距离成正比。注意，X轴和Y轴都表示生物活性距离。下面的例子也应用了一个小的位置抖动，这样代表具有相同生物活性的化合物的点就不会重叠。

库(ggplot2) #二维MDS转换plotdata <- cmdscale(as.dist(1- simma)， k=2) dat <- data.frame(xvar=plotdata[，1]， yvar=plotdata[，2]) #设置plot主题mytheme = theme(plot. 1)保证金=单位(c(。2。2。2,,。2),单位=“行”),轴。Text = element_blank()，轴。ticks = element_blank()， axis.ticks.length = unit(0， "lines")) # x和y变量的scatterplot minR <- min(range(dat$xvar)， range(dat$yvar)) - 0.1 maxR <- 0.2 + max(range(dat$xvar)， range(dat$yvar)) scatter <- ggplot(dat, aes(xvar, yvar)) + xlim(minR, maxR) + ylim(minR,maxR) + geom_point(shape=19, position = position_jitter(w = 0.05，h = 0.05)) + scale_colour_brewer(palette="Set1") + mytheme + theme(legend.position="none") + xlab("生物活性指纹距离")+ ylab("生物活性指纹距离")plot(scatter)

最后，断开与数据库的连接。

disconnectBioassayDB (sampleDB)

3．5分析和加载原始筛选数据

本例演示了分析和加载来自21,888种不同化合物的高通量筛选实验的数据的基础知识。

这个例子基于cellHTS2库(Boutros, L, and Huber 2006)．使用cellHTS2中包含的示例数据。这实际上是针对活细胞筛选dsRNA的数据，但由于数据格式相同，我们将其视为针对蛋白质目标的小分子结合数据。

首先读入cellHTS2提供的筛选数据。

- system. library(cellHTS2) library(bioassayR) dataPath <- system. library(cellHTS2)file("KcViab"， package="cellHTS2") x <- readPlateList("Platelist.txt"， name="KcViab"， path=dataPath) x <- configure(x, descripFile="Description.txt"， confFile="Plateconf.txt"， logFile="Screenlog.txt"， path=dataPath) xn <- normalizePlates(x, scale="multiplicative"， log=FALSE, method="median"， varianceAdjust="none")

接下来，对重复进行评分和总结。

xsc <- scorereplcopies (xn, sign="-"， method="zscore") xsc <- summarizereplcopies (xsc, summary="mean")

解析注释数据。

xsc <- annotate(xsc, geneIDFile="GeneIDs_Dm_HFA_1.1.txt"， path=dataPath)

应用s型变换生成二进制调用。

y <- scores2calls(xsc, z0=1.5, lambda=2) binaryCalls <- round(Data(y))

将二进制调用转换为一个活动表bioassayR可以解析。

scoreDataFrame <- cbind(geneAnno(y)， binaryCalls) rawScores <- as.vector(Data(xsc)) rawScores <- rawScores[wellAnno(y) == "sample"] scoreDataFrame <- scoreDataFrame[wellAnno(y) == "sample"，] activityTable <- cbind(cid=scoreDataFrame[，1]， activity=scoreDataFrame[，2]， score=rawScores) activityTable <- as.data.frame(activityTable) activityTable[1:10，]

## cid活动评分## 1 CG11371 1 2.147814893189 ## 2 CG31671 1 4.00105151825918 ## 3 CG11376 0 0.955550327380792 ## 4 CG11723 0 0.768879954745058 ## 5 CG11723 0 1.12105534792054 ## 6 CG7261 0 -0.126481089856157 ## 7 CG2674 0 0.706942495432899 ## 8 CG7263 0 0.920455054477106 ## 9 CG4822 0 0.2514405137389305878 ## 10 CG4265 0 0.414137389305878

创建一个新的(在本例中是临时的)bioassayR数据库来加载这些数据。

myDatabaseFilename <- tempfile() mydb <- newBioassayDB(myDatabaseFilename, indexed=F) addDataSource(mydb, description="other"， version="unknown")

为新分析创建一个分析对象。

myAssay <- new("bioassay"，aid="1"， source_id="other"， assay_type="confirmatory"， organism="unknown"， scoring="activity rank"， targets="2224444"， target_types="protein"， scores=activityTable)

将此分析对象加载到bioassayR数据库。

loadBioassay(mydb, myAssay) mydb

##类:BioassayDB ## assays: 1 ##源:其他##源版本:未知##可写:是

现在已经加载了这些数据，您可以使用它们来执行本文档中的任何其他分析示例。

最后，为了本例的目的，断开与示例数据库的连接。

disconnectBioassayDB (mydb)

3.6自定义SQL查询

虽然提供了许多预先构建的查询(请参阅其他示例和手册页)，但高级用户也可以构建自己的SQL查询。由于bioassayR使用SQLite数据库，您可以咨询http://www.sqlite.org有关构建SQL查询的详细信息。我们也推荐《使用SQLite》这本书。(Kreibich 2010)．

要开始，首先连接到数据库。如果你从作者的网站下载了完整的PubChem BioAssay数据库，变量sampleDatabasePath可以替换为您下载的数据库的文件名，如果您想使用该文件名而不是此软件包中包含的小型示例版本。

library(bioassayR) extdata_dir <- system. bat。sampleDatabasePath <- file. file("extdata"， package="bioassayR")path(extdata_dir， "sampleDatabase.sqlite") pubChemDatabase <- connectBioassayDB(sampleDatabasePath)

首先，你想看到数据库的结构如下所示:

queryBioassayDB(pubChemDatabase， "SELECT * FROM sqlite_master WHERE type='table'")

# # 1 # #类型名称tbl_name rootpage表活动活动2 # # 2表分析化验3 # # 3表源源4 # # 4表目标目标5 # # 5表targetTranslations targetTranslations 6 # # 1 # # sql创建表活动(活动援助整数,cid整数,整数,整数)得分# # 2创建表分析(source_id整数,整数的援助,assay_type文本、生物文本,得分文本)# # 3创建表来源(ASC source_id整数主键、描述文本,## 5创建表targetTranslations(目标文本，类别文本，标识符文本)

例如，你可以找到一个给定的化合物参与的前10个化验如下:

queryBioassayDB(pubChemDatabase， "SELECT DISTINCT(aid) FROM activity WHERE cid = '2244' LIMIT 10")

##援助## 1 399401 ## 2 399404 ## 3 399411 ## 4 443726 ## 5 410 ## 6 411 ## 7 422 ## 8 429 ## 9 436 ## 10 445

下面的例子打印了一个指定化验的活性分数:

queryBioassayDB(pubChemDatabase， "SELECT * FROM activity WHERE aid = '393818'")

## aid cid活动评分## 1 393818 3672 1 NA ## 2 393818 2662 1 NA ## 3 393818 25258347 NA NA ## 4 393818 25258348 1 NA ## 5 393818 25258349 1 NA ## 6 393818 25258350 1 NA

一个自然的加入自动合并共享公共行的表，从而更容易解析分布在多个表中的数据。在这里，我们将活性表(原始分数)与测定表(测定注释细节)和蛋白质靶标表合并:

queryBioassayDB(pubChemDatabase， "SELECT * FROM activity NATURAL JOIN assay NATURAL JOIN targets WHERE cid = '2244' LIMIT 10")

# #援助cid活动得分source_id assay_type有机体得分目标target_type # # 1 399401 2244 NA NA 1确认Bos_taurus IC50 166897622蛋白# # 2 399404 2244 NA NA 1确认Bos_taurus IC50 166897622蛋白# # 3 399411 2244 NA NA 1确认Bos_taurus IC50 2244蛋白质# # 4 443726 166897622 NA 1确认Bos_taurus IC50 2244蛋白质# # 5 410 166897622 0 0 1确认Homo_sapiens < NA > 2244蛋白质# # 6 411 73915100 0 0 1确认Photinus_pyralis效力160794蛋白7 422 2244 0 -6 1筛选Homo_sapiens  21464101蛋白8 429 2244 0 -2 1筛选Homo_sapiens  154146191蛋白9 429 2244 0 -2 1筛选Homo_sapiens  4261762蛋白10 436 2244 0 5 1筛选Homo_sapiens  68565074蛋白

最后，在分析后断开与数据库的连接可以减少数据损坏的机会。如果您正在使用只读模式下的预构建数据库(如预构建数据库示例部分所示)，则可以选择跳过此步骤，因为只有可写数据库容易因断开连接失败而损坏。

disconnectBioassayDB (pubChemDatabase)