doc/user/topics/tutorial-trait-search.html.textile.liquid

   1 ---
   2 layout: default
   3 navsection: userguide
   4 title: "Querying the Metadata Database"
   5 ...
   6 {% comment %}
   7 Copyright (C) The Arvados Authors. All rights reserved.
   8
   9 SPDX-License-Identifier: CC-BY-SA-3.0
  10 {% endcomment %}
  11
  12 This tutorial introduces the Arvados Metadata Database.  The Metadata Database stores information about files in Keep.  This example will use the Python SDK to find public WGS (Whole Genome Sequencing) data for people who have reported a certain medical condition.
  13
  14 {% include 'tutorial_expectations' %}
  15
  16 In the tutorial examples, three angle brackets (&gt;&gt;&gt;) will be used to denote code to enter at the interactive Python prompt.
  17
  18 Start by running Python.
  19
  20 <notextile>
  21 <pre><code>~$ <span class="userinput">python</span>
  22 Python 2.7.3 (default, Jan  2 2013, 13:56:14)
  23 [GCC 4.7.2] on linux2
  24 Type "help", "copyright", "credits" or "license" for more information.
  25 &gt;&gt;&gt;
  26 </code></pre>
  27 </notextile>
  28
  29 If everything is set up correctly, you will be able to import the arvados SDK.
  30
  31 notextile. <pre><code>&gt;&gt;&gt; <span class="userinput">import arvados</span></pre></code>
  32
  33 This tutorial will also use the regular expression (re) python module:
  34
  35 <notextile>
  36 <pre><code>&gt;&gt;&gt; <span class="userinput">import re</span>
  37 </code></pre>
  38 </notextile>
  39
  40 h2. Finding traits
  41
  42 notextile. <pre><code>&gt;&gt;&gt; <span class="userinput">all_traits = arvados.api().traits().list(limit=1000).execute()</span></code></pre>
  43
  44 * @arvados.api()@ gets an object that provides access to the Arvados API server
  45 * @.traits()@ gets an object that provides access to the "traits" resource on the Arvados API server
  46 * @.list(limit=1000)@ constructs a query to list all elements of the "traits" resource, with a limit of 1000 entries returned
  47 * @.execute()@ executes the query and returns the result, which we assign to "all_traits"
  48
  49 notextile. <pre><code>&gt;&gt;&gt; <span class="userinput">cancer_traits = filter(lambda t: re.search('cancer', t['name']), all_traits['items'])</span></code></pre>
  50
  51 * @lambda t: re.search('cancer', t['name'])@ is an inline function that takes a parameter @t@ and uses a simple regular expression to test if @t['name']@ contains the substring 'cancer'
  52 * @all_traits['items']@ is the input sequence of traits
  53 * @filter@ tests each element @t@ and constructs a new sequence consisting only of the elements that pass the filter
  54 * @cancer_traits@ gets the result of @filter@
  55
  56 <notextile>
  57 <pre><code>&gt;&gt;&gt; <span class="userinput">for t in cancer_traits: print(t['uuid'], t['name'])</span>
  58 ...
  59 qr1hi-q1cn2-8q57g2diohwnzm0 Cervical cancer
  60 qr1hi-q1cn2-vqp4243janpjbyj Breast cancer
  61 qr1hi-q1cn2-v6usijujcpwqrn1 Non-melanoma skin cancer
  62 ...
  63 </code></pre>
  64 </notextile>
  65
  66 In this tutorial wil will use "Non-melanoma skin cancer" trait with uuid @qr1hi-q1cn2-v6usijujcpwqrn1@.
  67
  68 notextile. <pre><code>&gt;&gt;&gt; <span class="userinput">non_melanoma_cancer = 'qr1hi-q1cn2-v6usijujcpwqrn1'</code></pre>
  69
  70 h2. Finding humans with the selected trait
  71
  72 We query the "links" resource to find humans that report the selected trait.  Links are directional connections between Arvados data items, for example, from a human to their reported traits.
  73
  74 <notextile>
  75 <pre><code>&gt;&gt;&gt; <span class="userinput">trait_filter = [
  76     ['link_class', '=', 'human_trait'],
  77     ['tail_uuid', 'is_a', 'arvados#human'],
  78     ['head_uuid', '=', non_melanoma_cancer],
  79   ]
  80 </code></pre>
  81 </notextile>
  82
  83 * @['link_class', '=', 'human_trait']@ filters on links that connect phenotype traits to individuals in the database.
  84 * @['tail_uuid', 'is_a', 'arvados#human']@ filters that the "tail" must be a "human" database object.
  85 * @['head_uuid', '=', non_melanoma_cancer]@ filters that the "head" of the link must connect to the "trait" database object non_melanoma_cancer .
  86
  87 The query will return links that match all three conditions.
  88
  89 <notextile>
  90 <pre><code>&gt;&gt;&gt; <span class="userinput">trait_links = arvados.api().links().list(limit=1000, filters=trait_filter).execute()</span>
  91 </code></pre>
  92 </notextile>
  93
  94 * @arvados.api()@ gets an object that provides access to the Arvados API server
  95 * @.links()@ gets an object that provides access to the "links" resource on the Arvados API server
  96 * @.list(limit=1000, filters=trait_filter)@ constructs a query to elements of the "links" resource that match the criteria discussed above, with a limit of 1000 entries returned
  97 * @.execute()@ executes the query and returns the result, which we assign to "trait_links"
  98
  99 <notextile>
 100 <pre><code>&gt;&gt;&gt; <span class="userinput">human_uuids = map(lambda l: l['tail_uuid'], trait_links['items'])</span>
 101 &gt;&gt;&gt; <span class="userinput">human_uuids</span>
 102 [u'1h9kt-7a9it-c0uqa4kcdh29wdf', u'1h9kt-7a9it-x4tru6mn40hc6ah',
 103 u'1h9kt-7a9it-yqb8m5s9cpy88i8', u'1h9kt-7a9it-46sm75w200ngwny',
 104 u'1h9kt-7a9it-gx85a4tdkpzsg3w', u'1h9kt-7a9it-8cvlaa8909lgeo9',
 105 u'1h9kt-7a9it-as37qum2pq8vizb', u'1h9kt-7a9it-14fph66z2baqxb9',
 106 u'1h9kt-7a9it-e9zc7i4crmw3v69', u'1h9kt-7a9it-np7f35hlijlxdmt',
 107 u'1h9kt-7a9it-j9hqyjwbvo9cojn', u'1h9kt-7a9it-lqxdtm1gynmsv13',
 108 u'1h9kt-7a9it-zkhhxjfg2o22ywq', u'1h9kt-7a9it-nsjoxqd33lzldw9',
 109 u'1h9kt-7a9it-ytect4smzcgd4kg', u'1h9kt-7a9it-y6tl353b3jc4tos',
 110 u'1h9kt-7a9it-98f8qave4f8vbs5', u'1h9kt-7a9it-gd72sh15q0p4wq3',
 111 u'1h9kt-7a9it-zlx25dscak94q9h', u'1h9kt-7a9it-8gronw4rbgmim01',
 112 u'1h9kt-7a9it-wclfkjcb23tr5es', u'1h9kt-7a9it-rvp2qe7szfz4dy6',
 113 u'1h9kt-7a9it-50iffhmpzsktwjm', u'1h9kt-7a9it-ul412id5y31a5o8',
 114 u'1h9kt-7a9it-732kwkfzylmt4ik', u'1h9kt-7a9it-v9zqxegpblsbtai',
 115 u'1h9kt-7a9it-kmaraqduit1v5wd', u'1h9kt-7a9it-t1nwtlo1hru5vvq',
 116 u'1h9kt-7a9it-q3w6j9od4ibpoyl', u'1h9kt-7a9it-qz8vzkuuz97ezwv',
 117 u'1h9kt-7a9it-t1v8sjz6dm9jmjf', u'1h9kt-7a9it-qe8wrbyvuqs5jew']
 118 </code></pre>
 119 </notextile>
 120
 121 * @lambda l: l['tail_uuid']@ is an inline function that returns the 'tail_uuid' attribute of 'l'
 122 * @trait_links['items']@ is the input set from the query
 123 * @map@ converts each item in a sequence into a different item using the embedded function, in this case to produce a sequence of uuids which refer to humans that have the specified trait.
 124
 125 h2. Find Personal Genome Project identifiers from Arvados UUIDs
 126
 127 <notextile>
 128 <pre><code>&gt;&gt;&gt; <span class="userinput">human_filters = [
 129     ["link_class", "=", "identifier"],
 130     ["head_uuid", "in", human_uuids]
 131   ]</span>
 132 &gt;&gt;&gt; <span class="userinput">pgpid_links = arvados.api('v1').links().list(limit=1000, filters=human_filters).execute()</span>
 133 &gt;&gt;&gt; <span class="userinput">map(lambda l: l['name'], pgpid_links['items'])</span>
 134 [u'hu01024B', u'hu11603C', u'hu15402B', u'hu174334', u'hu1BD549', u'hu237A50',
 135  u'hu34A921', u'hu397733', u'hu414115', u'hu43860C', u'hu474789', u'hu553620',
 136  u'hu56B3B6', u'hu5917F3', u'hu599905', u'hu5E55F5', u'hu602487', u'hu633787',
 137  u'hu68F245', u'hu6C3F34', u'hu7260DD', u'hu7A2F1D', u'hu94040B', u'hu9E356F',
 138  u'huAB8707', u'huB1FD55', u'huB4883B', u'huD09050', u'huD09534', u'huD3A569',
 139  u'huDF04CC', u'huE2E371']
 140 </code></pre>
 141 </notextile>
 142
 143 These PGP IDs let us find public profiles, for example:
 144
 145 * "https://my.pgp-hms.org/profile/huE2E371":https://my.pgp-hms.org/profile/huE2E371
 146 * "https://my.pgp-hms.org/profile/huDF04CC":https://my.pgp-hms.org/profile/huDF04CC
 147 * ...
 148
 149 h2. Find genomic data from specific humans
 150
 151 Now we want to find collections in Keep that were provided by these humans.  We search the "links" resource for "provenance" links that point to entries in the list of humans with the non-melanoma skin cancer trait:
 152
 153 <notextile>
 154 <pre><code>&gt;&gt;&gt; <span class="userinput">provenance_links = arvados.api().links().list(limit=1000, filters=[
 155     ["link_class", "=", "provenance"],
 156     ["name", "=", "provided"],
 157     ["tail_uuid", "in", human_uuids]
 158   ]).execute()
 159 collection_uuids = map(lambda l: l['head_uuid'], provenance_links['items'])
 160
 161 # build map of human uuid -> PGP ID
 162 pgpid = {}
 163 for pgpid_link in pgpid_links['items']:
 164   pgpid[pgpid_link['head_uuid']] = pgpid_link['name']
 165
 166 # build map of collection uuid -> PGP ID
 167 for p_link in provenance_links['items']:
 168   pgpid[p_link['head_uuid']] = pgpid[p_link['tail_uuid']]
 169
 170 # get details (e.g., list of files) of each collection
 171 collections = arvados.api('v1').collections().list(filters=[
 172     ["uuid", "in", collection_uuids]
 173   ]).execute()
 174
 175 # print PGP public profile links with file locators
 176 for c in collections['items']:
 177   for f in c['files']:
 178     print "https://my.pgp-hms.org/profile/%s %s %s%s" % (pgpid[c['uuid']], c['uuid'], ('' if f[0] == '.' else f[0]+'/'), f[1])
 179 </span>
 180 https://my.pgp-hms.org/profile/hu43860C a58dca7609fa84c8c38a7e926a97b2fc var-GS00253-DNA_A01_200_37-ASM.tsv.bz2
 181 https://my.pgp-hms.org/profile/huB1FD55 ea30eb9e46eedf7f05ed6e348c2baf5d var-GS000010320-ASM.tsv.bz2
 182 https://my.pgp-hms.org/profile/huDF04CC 4ab0df8f22f595d1747a22c476c05873 var-GS000010427-ASM.tsv.bz2
 183 https://my.pgp-hms.org/profile/hu7A2F1D 756d0ada29b376140f64e7abfe6aa0e7 var-GS000014566-ASM.tsv.bz2
 184 https://my.pgp-hms.org/profile/hu553620 7ed4e425bb1c7cc18387cbd9388181df var-GS000015272-ASM.tsv.bz2
 185 https://my.pgp-hms.org/profile/huD09534 542112e210daff30dd3cfea4801a9f2f var-GS000016374-ASM.tsv.bz2
 186 https://my.pgp-hms.org/profile/hu599905 33a9f3842b01ea3fdf27cc582f5ea2af var-GS000016015-ASM.tsv.bz2
 187 https://my.pgp-hms.org/profile/hu43860C a58dca7609fa84c8c38a7e926a97b2fc+302 var-GS00253-DNA_A01_200_37-ASM.tsv.bz2
 188 https://my.pgp-hms.org/profile/huB1FD55 ea30eb9e46eedf7f05ed6e348c2baf5d+291 var-GS000010320-ASM.tsv.bz2
 189 https://my.pgp-hms.org/profile/huDF04CC 4ab0df8f22f595d1747a22c476c05873+242 var-GS000010427-ASM.tsv.bz2
 190 https://my.pgp-hms.org/profile/hu7A2F1D 756d0ada29b376140f64e7abfe6aa0e7+242 var-GS000014566-ASM.tsv.bz2
 191 https://my.pgp-hms.org/profile/hu553620 7ed4e425bb1c7cc18387cbd9388181df+242 var-GS000015272-ASM.tsv.bz2
 192 https://my.pgp-hms.org/profile/huD09534 542112e210daff30dd3cfea4801a9f2f+242 var-GS000016374-ASM.tsv.bz2
 193 https://my.pgp-hms.org/profile/hu599905 33a9f3842b01ea3fdf27cc582f5ea2af+242 var-GS000016015-ASM.tsv.bz2
 194 https://my.pgp-hms.org/profile/hu599905 d6e2e57cd60ba5979006d0b03e45e726+81 Witch_results.zip
 195 https://my.pgp-hms.org/profile/hu553620 ea4f2d325592a1272f989d141a917fdd+85 Devenwood_results.zip
 196 https://my.pgp-hms.org/profile/hu7A2F1D 4580f6620bb15b25b18373766e14e4a7+85 Innkeeper_results.zip
 197 https://my.pgp-hms.org/profile/huD09534 fee37be9440b912eb90f5e779f272416+82 Hallet_results.zip
 198 </code></pre>
 199 </notextile>
 200
 201 h3. Search for a variant
 202
 203 Now we will use crunch to issue a 'grep' job to look for variant rs1126809 in each of the "var-" files (these contain variant calls from WGS data).
 204
 205 <notextile>
 206 <pre><code>&gt;&gt;&gt; <span class="userinput">job = {}
 207 for c in collections['items']:
 208   if [] != filter(lambda f: re.search('^var-.*\.tsv\.bz2', f[1]), c['files']):
 209     job[c['uuid']] = arvados.api('v1').jobs().create(body={
 210       'script': 'grep',
 211       'script_parameters': {'input': c['uuid'], 'pattern': "rs1126809\\b"},
 212       'script_version': 'e7aeb42'
 213     }).execute()
 214     print "%s %s" % (pgpid[c['uuid']], job[c['uuid']]['uuid'])
 215 </span>
 216 hu43860C qr1hi-8i9sb-wbf3uthbhkcy8ji
 217 huB1FD55 qr1hi-8i9sb-scklkiy8dc27dab
 218 huDF04CC qr1hi-8i9sb-pg0w4rfrwfd9srg
 219 hu7A2F1D qr1hi-8i9sb-n7u0u0rj8b47168
 220 hu553620 qr1hi-8i9sb-k7gst7vyhg20pt1
 221 huD09534 qr1hi-8i9sb-4w65pm48123fte5
 222 hu599905 qr1hi-8i9sb-wmwa5b5r3eghnev
 223 hu43860C qr1hi-8i9sb-j1mngmakdh8iv9o
 224 huB1FD55 qr1hi-8i9sb-4j6ehiatcolaoxb
 225 huDF04CC qr1hi-8i9sb-n6lcmcr3lowqr5u
 226 hu7A2F1D qr1hi-8i9sb-0hwsdtojfcxjo40
 227 hu553620 qr1hi-8i9sb-cvvqzqea7jhwb0i
 228 huD09534 qr1hi-8i9sb-d0y0qtzuwzbrjj0
 229 hu599905 qr1hi-8i9sb-i9ec9g8d7rt70xg
 230 </code></pre>
 231 </notextile>
 232
 233
 234 Monitor job progress by refreshing the Jobs page in Workbench, or by using the API:
 235
 236 <notextile>
 237 <pre><code>&gt;&gt;&gt; <span class="userinput">map(lambda j: arvados.api('v1').jobs().get(uuid=j['uuid']).execute()['success'], job.values())
 238 [None, True, None, None, None, None, None, None, None, None, None, None, None, None]
 239 </code></pre>
 240 </notextile>
 241
 242 Unfinished jobs will appear as None, failed jobs as False, and completed jobs as True.
 243
 244 After the jobs have completed, check output file sizes.
 245
 246 <notextile>
 247 <pre><code>&gt;&gt;&gt; <span class="userinput">for collection_uuid in job:
 248   job_uuid = job[collection_uuid]['uuid']
 249   job_output = arvados.api('v1').jobs().get(uuid=job_uuid).execute()['output']
 250   output_files = arvados.api('v1').collections().get(uuid=job_output).execute()['files']
 251   # Test the output size.  If greater than zero, that means 'grep' found the variant
 252   if output_files[0][2] > 0:
 253     print("%s has variant rs1126809" % (pgpid[collection_uuid]))
 254   else:
 255     print("%s does not have variant rs1126809" % (pgpid[collection_uuid]))
 256 </span>
 257 hu553620 does not have variant rs1126809
 258 hu43860C does not have variant rs1126809
 259 hu599905 has variant rs1126809
 260 huD09534 has variant rs1126809
 261 hu553620 does not have variant rs1126809
 262 huB1FD55 does not have variant rs1126809
 263 huDF04CC has variant rs1126809
 264 hu7A2F1D has variant rs1126809
 265 hu7A2F1D has variant rs1126809
 266 hu599905 has variant rs1126809
 267 huDF04CC has variant rs1126809
 268 huB1FD55 does not have variant rs1126809
 269 huD09534 has variant rs1126809
 270 hu43860C does not have variant rs1126809
 271 </code></pre>
 272 </notextile>
 273
 274 Thus, of the 14 WGS results available for PGP participants reporting non-melanoma skin cancer, 8 include the rs1126809 variant.