BH16.12/MachineLearning

提供：TogoWiki

（版間での差分）
移動： 案内, 検索
(→実行)
(→実行)
 
565行： 565行：
2          0.367239    0.363791              0.834857      0.834586                  0.227231       
2          0.367239    0.363791              0.834857      0.834586                  0.227231       
3          0.340371    0.348169              0.835143      0.834586                  0.337782       
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       
97          0.000272794  0.968715              1              0.843739                  10.7034       

2016年12月16日 (金) 02:21時点における最新版

BH16.12

JSTシソーラス - MeSH - 遺伝子 - 様々な特徴 を学習して、遺伝子と表現型（病気など）のアノテーション（関係）を見つける機械学習 (AI) をつくる

  • 参加者：金城・片山
  • サポート：渡辺・櫛田

参考：Robert's paper

目次

特徴ベクトル

  • JSTシソーラスのRDFから特徴ベクトルを作成
    • ランダムウォークでタームごとの URI 周辺のグラフパターンを学習
    • ターム毎に特徴ベクトルを生成
  • 遺伝子アノテーションの特徴ベクトルを生成
    • UniProtやTogoGenomeのエントリからMeSHを含むリンクを抽出
    • MeSHとJSTのタームの対応を学習

データセット

  • 遺伝子アノテーション
    • 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） ため

のうち、旧版の 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

  1. Robert の前処理プログラム RDFWrapper がうまく動かない。Java のライブラリの問題？→しょうがないので、自前で前処理プログラムを作り、originalのDeepWalkプログラムを使うことにする。
  2. OCaml RDF libraryの都合により、ntriples -> turtleに変換。（rapper を使ったが、文字コードの問題で一部トリプルが省かれるみたい。）
  3. jst-mesh2016.ntとmesh2016.ntをマージしたグラフのデータを作った。ノード数 5,401,382
  4. DeepWalkを走らせてベクトルを作り始めた(2016/12/13 3PMごろ）。→まだ終わらない(17:52)　→翌朝には終わっていた。

Prediction

あるタンパク質が、ある病気に関連するか否かを予測する（？）

  • 必要なデータ
- positive/negative sets
- training/testing sets
  • UniProt のヒトタンパク質に関連づけられているMeSH(病気）トップ２０

全部で973 Mesh terms のうち357個は１つのUniProtエントリにしか関連づけられていない。

  1. 151 https://id.nlm.nih.gov/mesh/D008607 "Intellectual Disability"
  2. 106 https://id.nlm.nih.gov/mesh/D006319 "Hearing Loss, Sensorineural"
  3. 83 https://id.nlm.nih.gov/mesh/D012174
  4. 59 https://id.nlm.nih.gov/mesh/D007153
  5. 57 https://id.nlm.nih.gov/mesh/D008831
  6. 53 https://id.nlm.nih.gov/mesh/D038901
  7. 51 https://id.nlm.nih.gov/mesh/D015419
  8. 50 https://id.nlm.nih.gov/mesh/D002607
  9. 49 https://id.nlm.nih.gov/mesh/D017237
  10. 49 https://id.nlm.nih.gov/mesh/D010009
  11. 49 https://id.nlm.nih.gov/mesh/D000015
  12. 46 https://id.nlm.nih.gov/mesh/D028361
  13. 45 https://id.nlm.nih.gov/mesh/D005124
  14. 45 https://id.nlm.nih.gov/mesh/D002386
  15. 43 https://id.nlm.nih.gov/mesh/D052177
  16. 42 https://id.nlm.nih.gov/mesh/D004392
  17. 41 https://id.nlm.nih.gov/mesh/D018981
  18. 39 https://id.nlm.nih.gov/mesh/D002524
  19. 39 https://id.nlm.nih.gov/mesh/D002311
  20. 37 https://id.nlm.nih.gov/mesh/D004476


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      
...
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
「/mw/BH16.12/MachineLearning」より作成