doc/user/tutorials/tutorial-trait-search.textile

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