BH16.12/MachineLearning
提供：TogoWiki
JSTシソーラス - MeSH - 遺伝子 - 様々な特徴 を学習して、遺伝子と表現型（病気など）のアノテーション（関係）を見つける機械学習 (AI) をつくる
- 参加者：金城・片山
- サポート：渡辺・櫛田
目次 |
特徴ベクトル
- JSTシソーラスのRDFから特徴ベクトルを作成
- ランダムウォークでタームごとの URI 周辺のグラフパターンを学習
- ターム毎に特徴ベクトルを生成
- 遺伝子アノテーションの特徴ベクトルを生成
- UniProtやTogoGenomeのエントリからMeSHを含むリンクを抽出
- MeSHとJSTのタームの対応を学習
データセット
- JSTシソーラス
- https://stirdf.jglobal.jst.go.jp/page/lodchallenge_2016.html から利用申請
- 使いたいトリプル
- rdfs:seeAlso
- dct:subject
- skos:related
- skos:broader
- skos:narrower
- jst-csv2nt.rb で jst-mesh.nt を生成 (62MB; gzip 872KB)
- JSTシソーラスに skos:Concept のついているタームは 244832 個
- このうち、rdfs:seeAlso で MeSH への対応がついているものが 15945 個
- スキーマ図 を参照し、これの dct:subject, skos:related, broader, narrower のトリプルを取得
- これにMeSHのRDFを加える？ https://id.nlm.nih.gov/mesh/
- MeSH RDF trivia
- 超最新版のURI https://id.nlm.nih.gov/mesh/D017538
- 最新版のURL http://id.nlm.nih.gov/mesh/D017538
- 2016年版のURL http://id.nlm.nih.gov/mesh/2016/D017538 <- "/2016" が余分に入る。
- JSTのデータは2016年版で「最新版」の形式のリンクが入っているが、本日(2016/12/12)時点での最新版は2017年版なので、整合性が取れない…要Hack
- MeSH RDF trivia
- 遺伝子アノテーション
- TogoGenome, DDBJ, PDBj, UniProt などの RDF から、遺伝子の特徴に関わるトリプルと MeSH へのリンクを収集
UniProt
UniProt から MeSH に辿るには、下記の predicate を辿れば OK そう。disease と mesh の対応は diseases.rdf に含まれている。
<uniprot/#> :annotation <SHA-384/#> :disease <diseases/#> rdfs:seeAlso <mesh/#>
例としては
prefix : <http://purl.uniprot.org/core/>
select *
where {
<http://purl.uniprot.org/uniprot/Q9UDR5> :annotation ?a .
?a :disease ?d .
?d rdfs:seeAlso ?m .
FILTER strstarts(str(?m), "https://id.nlm.nih.gov/mesh/")
}
要注意なのは、現状の TogoGenome のエンドポイントのデータは少し古い（BH16.12/TogoGenome） ため
- 最新版の MeSH URI https://id.nlm.nih.gov/mesh/D020167
- 旧版の MeSH URI http://purl.uniprot.org/mesh/D020167
のうち、旧版の URI を使ったものになる（いずれにしても JST シソーラスのために URI を書き換える必要があるのと、2016 年前半時点のデータを使いたいのである意味ちょうどよかった？）。
prefix : <http://purl.uniprot.org/core/>
construct {
?u :annotation ?a .
?a :disease ?d .
?d rdfs:seeAlso ?mesh2016 .
}
where {
VALUES ?u { <http://purl.uniprot.org/uniprot/Q9UDR5> }
?u :annotation ?a .
?a :disease ?d .
?d rdfs:seeAlso ?m .
FILTER strstarts(str(?m), "http://purl.uniprot.org/mesh/")
BIND (iri(concat("https://id.nlm.nih.gov/mesh/", "2016/", strafter(str(?m), "mesh/"))) AS ?mesh2016)
}
MeSH 以外にどのような情報をひろうかを検討 → タンパク個別の情報は拾いたいが意味的に関連しないものにやたらとつながりがちな provenance/evidence 系の情報は拾わないことにするため、
- :attribution
- :citation
は無視することにし、あまり意味のない
- :alternativeName
- :conflictingSequence
- :created
- :modified
- :replaces
- :reviewed
- :version
も捨てて、
<uniprot/Q9UDR5> (a :Protein) :organism <taxonomy/9606> :mnemonic "AASS_HUMAN" :recommendedName/:fullName "Alpha-aminoadipic semialdehyde synthase, mitochondrial" (a :Structured_Name) :annotation/:disease <diseases/1773> (a :Disease_Annotation) :annotation/rdfs:comment "The disease is caused by mutations affecting the gene represented in this entry." (a :Disease_Annotation) :annotation/rdfs:comment "N6-acetyllysine; alternate" (a :Modified_Residue_Annotation) :annotation/rdfs:comment "In the N-terminal section; belongs to the AlaDH/PNT family." (a :Similarity_Annotation) :annotation/rdfs:comment "Saccharopine dehydrogenase" (a :Region_Annotation) :annotation/rdfs:comment "Induced by starvation." (a :Induction_Annotation) :annotation/rdfs:comment "Expressed in all 16 tissues examined with highest expression in the liver." (a :Tissue_Specificity_Annotation) :annotation/rdfs:comment "Mitochondrion" (a :Transit_Peptide_Annotation) :annotation/rdfs:comment "N(6)-(L-1,3-dicarboxypropyl)-L-lysine + NADP(+) + H(2)O = L-lysine + 2-oxoglutarate + NADPH." (a :Catalytic_Activity_Annotation) :annotation/rdfs:comment "Lysine-ketoglutarate reductase" (a :Region_Annotation) :annotation/rdfs:comment "Homodimer." (a :Subunit_Annotation) :annotation/rdfs:comment "Bifunctional enzyme that catalyzes the first two steps in lysine degradation. The N-terminal and the C-terminal contain lysine-ketoglutarate reductase and saccharopine dehydrogenase activity, respectively." (a :Function_Annotation) :annotation/rdfs:comment "In the C-terminal section; belongs to the saccharopine dehydrogenase family." (a :Similarity_Annotation) :annotation/rdfs:comment "Alpha-aminoadipic semialdehyde synthase, mitochondrial" (a :Chain_Annotation) :annotation/:locatedIn/:cellularComponent <locations/173> (a :Subcellular_Location_Annotation) :annotation/rdfs:seeAlso <unipathway/427.868.419.835> (a :Pathway_Annotation) :classifiedWith <http://purl.obolibrary.org/obo/GO_0004753> :classifiedWith <keywords/1185> :domain/:enzyme <enzyme/1.5.1.8> (a :Part) :encodedBy/skos:prefLabel "AASS" (a :Gene) :existence :Evidence_at_Protein_Level_Existence :isolatedFrom <tissues/564> :proteome <proteomes/UP000005640#Chromosome%207> :sequence <isoforms/Q9UDR5-1> :sequence/:mass "102132" (a :Simple_Sequence) :sequence/rdf:value ""MLQVHRTGLGRLGV ...." rdfs:seeAlso <http://identifiers.org/ncbigene/10157> rdfs:seeAlso/:database <database/DisGeNET> (a :Resource) rdfs:seeAlso <http://identifiers.org/reactome/R-HSA-71064> rdfs:seeAlso/:database <database/Reactome> (a Resource) rdfs:seeAlso/rdfs:comment "Lysine catabolism" (a Resource) :
このあたりを <SHA-384/#> も含めて拾うことに。
ただ、http://www.uniprot.org/uniprot/Q9UDR5 からダウンロードした RDF/XML に含まれる SHA-384 の URI をブラウザで開くと DESCRIBEが実行されるが、
select *
where {
?s ?p <http://purl.uniprot.org/SHA-384/ACB87331D17D2460F1FA8D686155044D4053D0FC575AEB433F5D74F2CFAF5A84B1449AE0068947CE4675EE8BB71B0A26>
}
limit 100
select *
where {
?s ?p <http://purl.uniprot.org/SHA-384/E2AA67B236D056422CB4A02B42561C44A901BADA4F19813CF481F77D6DFB1E7B193931927AE9E39A972E74F293E44890>
}
limit 100
があり、ヒットしない URI は SPARQL で使っても当然出てこない。最新のエントリに含まれる SHA-384 の ID と UniProt の公開エンドポイントに入っている RDF データは同期されていない（古い）のだろうか？
PREFIX : <http://purl.uniprot.org/core/>
CONSTRUCT {
?uniprot :organism ?organism .
?uniprot :annotation ?annotation .
?annotation :disease ?disease .
?disease rdfs:seeAlso ?mesh2016 .
?uniprot :mnemonic ?mnemonic .
?uniprot :recommendedName ?recommendedName .
?recommendedName :fullName ?fullName .
?annotation rdfs:comment ?comment .
?annotation :locatedIn ?locatedIn .
?locatedIn :cellularComponent ?cellularComponent .
?annotation rdfs:seeAlso ?pathway .
?uniprot :classifiedWith ?classifiedWith .
?uniprot :domain ?domain .
?domain :enzyme ?enzyme .
?uniprot :encodedBy ?encodedBy .
?encodedBy skos:prefLabel ?prefLabel .
?uniprot :existence ?existence .
?uniprot :isolatedFrom ?isolatedFrom .
?uniprot :proteome ?proteome .
?uniprot :sequence ?sequence .
?sequence :mass ?mass .
?uniprot rdfs:seeAlso ?seeAlso .
?seeAlso :database ?database .
?seeAlso rdfs:comment ?comment .
}
WHERE {
VALUES ?organism { <http://purl.uniprot.org/taxonomy/9606> }
?uniprot :organism ?organism .
### 以下は全て OPTIONAL に ###
?uniprot :annotation ?annotation .
# MeSH との対応
?annotation :disease ?disease .
?disease rdfs:seeAlso ?mesh .
FILTER strstarts(str(?mesh), "http://purl.uniprot.org/mesh/")
BIND (iri(concat("https://id.nlm.nih.gov/mesh/", "2016/", strafter(str(?mesh), "mesh/"))) AS ?mesh2016)
# その他
?uniprot :mnemonic ?mnemonic .
?uniprot :recommendedName ?recommendedName .
?recommendedName :fullName ?fullName .
?annotation rdfs:comment ?comment .
?annotation :locatedIn ?locatedIn .
?locatedIn :cellularComponent ?cellularComponent .
?annotation rdfs:seeAlso ?pathway .
?uniprot :classifiedWith ?classifiedWith .
?uniprot :domain ?domain .
?domain :enzyme ?enzyme .
?uniprot :encodedBy ?encodedBy .
?encodedBy skos:prefLabel ?prefLabel .
?uniprot :existence ?existence .
?uniprot :isolatedFrom ?isolatedFrom .
?uniprot :proteome ?proteome .
?uniprot :sequence ?sequence .
?sequence :mass ?mass .
?uniprot rdfs:seeAlso ?seeAlso .
?seeAlso :database ?database .
?seeAlso rdfs:comment ?comment .
}
ということで
PREFIX : <http://purl.uniprot.org/core/>
CONSTRUCT {
?uniprot :organism ?organism .
?uniprot :annotation ?annotation .
?annotation :disease ?disease .
?disease rdfs:seeAlso ?mesh2016 .
?uniprot :mnemonic ?mnemonic .
?uniprot :recommendedName ?recommendedName .
?recommendedName :fullName ?fullName .
?annotation rdfs:comment ?comment .
?annotation :locatedIn ?locatedIn .
?locatedIn :cellularComponent ?cellularComponent .
?annotation rdfs:seeAlso ?pathway .
?uniprot :classifiedWith ?classifiedWith .
?uniprot :domain ?domain .
?domain :enzyme ?enzyme .
?uniprot :encodedBy ?encodedBy .
?encodedBy skos:prefLabel ?prefLabel .
?uniprot :existence ?existence .
?uniprot :isolatedFrom ?isolatedFrom .
?uniprot :proteome ?proteome .
?uniprot :sequence ?sequence .
?sequence :mass ?mass .
?uniprot rdfs:seeAlso ?seeAlso .
?seeAlso :database ?database .
?seeAlso rdfs:comment ?comment .
}
WHERE {
VALUES ?organism { <http://purl.uniprot.org/taxonomy/9606> }
?uniprot :organism ?organism .
?uniprot :mnemonic ?mnemonic .
?uniprot :recommendedName ?recommendedName .
?recommendedName :fullName ?fullName .
# 以下は全て OPTIONAL に
OPTIONAL {
?uniprot :annotation ?annotation .
OPTIONAL {
# MeSH との対応
?annotation :disease ?disease .
?disease rdfs:seeAlso ?mesh .
FILTER strstarts(str(?mesh), "http://purl.uniprot.org/mesh/")
BIND (iri(concat("https://id.nlm.nih.gov/mesh/", "2016/", strafter(str(?mesh), "mesh/"))) AS ?mesh2016)
}
# その他の annotation
OPTIONAL {
?annotation rdfs:comment ?comment .
}
OPTIONAL {
?annotation :locatedIn ?locatedIn .
?locatedIn :cellularComponent ?cellularComponent .
}
OPTIONAL {
?annotation rdfs:seeAlso ?pathway .
}
}
OPTIONAL {
?uniprot :classifiedWith ?classifiedWith .
}
OPTIONAL {
?uniprot :domain ?domain .
?domain :enzyme ?enzyme .
}
OPTIONAL {
?uniprot :encodedBy ?encodedBy .
?encodedBy skos:prefLabel ?prefLabel .
}
OPTIONAL {
?uniprot :existence ?existence .
}
OPTIONAL {
?uniprot :isolatedFrom ?isolatedFrom .
}
OPTIONAL {
?uniprot :proteome ?proteome .
}
OPTIONAL {
?uniprot :sequence ?sequence .
OPTIONAL {
?sequence :mass ?mass .
}
}
OPTIONAL {
?uniprot rdfs:seeAlso ?seeAlso .
OPTIONAL {
?seeAlso :database ?database .
}
OPTIONAL {
?seeAlso rdfs:comment ?comment .
}
}
}
これを SPARQL エンドポイントで実行しても LIMIT 10000 くらいでないとレスポンスが帰ってこないので、construct-optional.rq として保存し、isql から実行することに。
#!/usr/bin/env ruby
head = <<HEAD
sparql
DEFINE output:format "NT"
HEAD
body = File.read("construct-optional.rq")
offset = 1
limit = 100000
input = "input.rq"
output = "output.nt"
before = 0
after = 1
while before < after do
tail = "OFFSET #{offset} LIMIT #{limit};"
File.open(input, "w") do |query|
query.puts head
query.puts body
query.puts tail
end
before = File.size(output)
system("isql VERBOSE=OFF BANNER=OFF PROMPT=OFF ECHO=OFF BLOBS=ON < #{input} >> #{output}")
after = File.size(output)
$stderr.puts "OFFSET #{offset} LIMIT #{limit} before: #{before} after: #{after}"
offset += limit
end
機械学習
DeepWalk
- Robert's implementation: https://github.com/bio-ontology-research-group/walking-rdf-and-owl
- DeepWalk: paper and implementation.
- Robert の前処理プログラム RDFWrapper がうまく動かない。Java のライブラリの問題？→しょうがないので、自前で前処理プログラムを作り、originalのDeepWalkプログラムを使うことにする。
- OCaml RDF libraryの都合により、ntriples -> turtleに変換。（rapper を使ったが、文字コードの問題で一部トリプルが省かれるみたい。）
- jst-mesh2016.ntとmesh2016.ntをマージしたグラフのデータを作った。ノード数 5,401,382
- DeepWalkを走らせてベクトルを作り始めた(2016/12/13 3PMごろ）。→まだ終わらない(17:52) →翌朝には終わっていた。
Prediction
あるタンパク質が、ある病気に関連するか否かを予測する（？）
- 必要なデータ
- - positive/negative sets
- - training/testing sets
- UniProt のヒトタンパク質に関連づけられているMeSH(病気）トップ２０
全部で973 Mesh terms のうち357個は１つのUniProtエントリにしか関連づけられていない。
- 151 https://id.nlm.nih.gov/mesh/D008607 "Intellectual Disability"
- 106 https://id.nlm.nih.gov/mesh/D006319 "Hearing Loss, Sensorineural"
- 83 https://id.nlm.nih.gov/mesh/D012174
- 59 https://id.nlm.nih.gov/mesh/D007153
- 57 https://id.nlm.nih.gov/mesh/D008831
- 53 https://id.nlm.nih.gov/mesh/D038901
- 51 https://id.nlm.nih.gov/mesh/D015419
- 50 https://id.nlm.nih.gov/mesh/D002607
- 49 https://id.nlm.nih.gov/mesh/D017237
- 49 https://id.nlm.nih.gov/mesh/D010009
- 49 https://id.nlm.nih.gov/mesh/D000015
- 46 https://id.nlm.nih.gov/mesh/D028361
- 45 https://id.nlm.nih.gov/mesh/D005124
- 45 https://id.nlm.nih.gov/mesh/D002386
- 43 https://id.nlm.nih.gov/mesh/D052177
- 42 https://id.nlm.nih.gov/mesh/D004392
- 41 https://id.nlm.nih.gov/mesh/D018981
- 39 https://id.nlm.nih.gov/mesh/D002524
- 39 https://id.nlm.nih.gov/mesh/D002311
- 37 https://id.nlm.nih.gov/mesh/D004476
- 機械学習のライブラリ
- TensorFlow
- Chainer
- SVM?
Chainerによる実装
- 使ったデータ
- - UniProt-human-subset.nt （片山さん製作 12/15 18:30現在まだ一部のみ）
- - mesh2016.nt (MeSHからとってきた）
- - go.nt (from UniProt)
- - pathways.nt (from UniProt)
- - enzyme.nt (from UniProt)
- - jst-mesh2016.nt (JST thesaurus)
これらを一つにまとめて、DeepWalkを走らせる。
- DeepWalk
% ./src/rdfwrapper -map all.map all.nt > all.edgelist # 自前のプログラムで前処理 % deepwalk --input all.edgelist --format edgelist --workers 4 --output all.vec --max-memory-data-size 3000000000
チュートリアルのコードをちょっとだけ書き換えたもの。
#!/usr/bin/env python
from __future__ import print_function
import argparse
import numpy as np
import chainer
import chainer.functions as F
import chainer.links as L
from chainer.datasets import tuple_dataset
from chainer import training
from chainer.training import extensions
# Network definition
class MLP(chainer.Chain):
def __init__(self, n_units, n_out):
super(MLP, self).__init__(
# the size of the inputs to each layer will be inferred
l1=L.Linear(None, n_units), # n_in -> n_units
l2=L.Linear(None, n_units), # n_units -> n_units
l3=L.Linear(None, n_out), # n_units -> n_out
)
def __call__(self, x):
h1 = F.relu(self.l1(x))
h2 = F.relu(self.l2(h1))
y = F.relu(self.l3(h2))
# print("hoge {}".format(y.data))
return y
def main():
parser = argparse.ArgumentParser(description='Chainer example: MNIST')
parser.add_argument('--batchsize', '-b', type=int, default=100,
help='Number of images in each mini-batch')
parser.add_argument('--epoch', '-e', type=int, default=20,
help='Number of sweeps over the dataset to train')
parser.add_argument('--gpu', '-g', type=int, default=-1,
help='GPU ID (negative value indicates CPU)')
parser.add_argument('--out', '-o', default='result',
help='Directory to output the result')
parser.add_argument('--resume', '-r', default='',
help='Resume the training from snapshot')
parser.add_argument('--unit', '-u', type=int, default=100,
help='Number of units')
args = parser.parse_args()
print('GPU: {}'.format(args.gpu))
print('# unit: {}'.format(args.unit))
print('# Minibatch-size: {}'.format(args.batchsize))
print('# epoch: {}'.format(args.epoch))
print('')
# Set up a neural network to train
# Classifier reports softmax cross entropy loss and accuracy at every
# iteration, which will be used by the PrintReport extension below.
model = L.Classifier(MLP(args.unit, 2))
if args.gpu >= 0:
chainer.cuda.get_device(args.gpu).use() # Make a specified GPU current
model.to_gpu() # Copy the model to the GPU
# Setup an optimizer
optimizer = chainer.optimizers.Adam()
optimizer.setup(model)
# データの入力：ここを書き換えた
# (URI) (UniProt エントリのベクトル64次元）（Meshのベクトル64次元）（０または１）の形式の空白区切りファイルを読み込む。
train_data = np.loadtxt("train.vec",delimiter=" ",skiprows=1,dtype=np.float32)
test_data= np.loadtxt("test.vec",delimiter=" ",skiprows=1,dtype=np.float32)
train = tuple_dataset.TupleDataset(train_data[:,1:-1],np.int32(train_data[:,-1]))
test = tuple_dataset.TupleDataset(test_data[:,1:-1],np.int32(test_data[:,-1]))
train_iter = chainer.iterators.SerialIterator(train, args.batchsize)
test_iter = chainer.iterators.SerialIterator(test, args.batchsize,
repeat=False, shuffle=False)
# Set up a trainer
updater = training.StandardUpdater(train_iter, optimizer, device=args.gpu)
trainer = training.Trainer(updater, (args.epoch, 'epoch'), out=args.out)
# Evaluate the model with the test dataset for each epoch
trainer.extend(extensions.Evaluator(test_iter, model, device=args.gpu))
# Dump a computational graph from 'loss' variable at the first iteration
# The "main" refers to the target link of the "main" optimizer.
trainer.extend(extensions.dump_graph('main/loss'))
# Take a snapshot at each epoch
trainer.extend(extensions.snapshot(), trigger=(args.epoch, 'epoch'))
# Write a log of evaluation statistics for each epoch
trainer.extend(extensions.LogReport())
# Print selected entries of the log to stdout
# Here "main" refers to the target link of the "main" optimizer again, and
# "validation" refers to the default name of the Evaluator extension.
# Entries other than 'epoch' are reported by the Classifier link, called by
# either the updater or the evaluator.
trainer.extend(extensions.PrintReport(
['epoch', 'main/loss', 'validation/main/loss',
'main/accuracy', 'validation/main/accuracy', 'elapsed_time']))
# Print a progress bar to stdout
trainer.extend(extensions.ProgressBar())
if args.resume:
# Resume from a snapshot
chainer.serializers.load_npz(args.resume, trainer)
# Run the training
trainer.run()
if __name__ == '__main__':
main()
実行
動いた！
[akinjo@cat ~/work/bh16.12/machinelearning]$ python ./chainer/train_deepwalk.py --epoch 100 GPU: -1 # unit: 100 # Minibatch-size: 100 # epoch: 100 epoch main/loss validation/main/loss main/accuracy validation/main/accuracy elapsed_time 1 0.492939 0.386433 0.824571 0.834586 0.114298 2 0.367239 0.363791 0.834857 0.834586 0.227231 3 0.340371 0.348169 0.835143 0.834586 0.337782 4 0.32338 0.332003 0.830588 0.834586 0.444711 5 0.299779 0.326216 0.834857 0.834586 0.552648 6 0.288021 0.321101 0.832571 0.834586 0.661041 7 0.283967 0.329176 0.832647 0.834586 0.770286 8 0.270831 0.33147 0.835429 0.834586 0.877979 9 0.269528 0.326854 0.833143 0.834586 0.987836 10 0.258968 0.317992 0.833824 0.834586 1.09585 11 0.250275 0.323071 0.832286 0.834586 1.20431 12 0.239626 0.318876 0.839714 0.841058 1.31732 13 0.234115 0.335189 0.882647 0.838183 1.44343 14 0.208163 0.342146 0.906857 0.832892 1.55745 15 0.19282 0.350635 0.915143 0.846878 1.67199 ... 97 0.000272794 0.968715 1 0.843739 10.7034 98 0.000258798 0.975599 1 0.847072 10.8188 99 0.000256288 0.979051 1 0.848183 10.9319 100 0.000244974 0.983407 1 0.849295 11.0407
