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

   1 ---
   2 layout: default
   3 navsection: userguide
   4 title: "Search PGP data by trait"
   5 navorder: 116
   6 ---
   7
   8 h1. Tutorial: Search PGP data by trait
   9
  10 Here you will use the Python SDK to find public WGS data for people who have a certain medical condition.
  11
  12 <!-- _Define WGS_ -->
  13
  14 <!-- _Explain the motivation in this example a little better.  If I'm
  15 reading this right, the workflow is
  16 traits -> people with those traits -> presense of a specific genetic
  17 variant in the people with the reported traits_ -->
  18
  19 <!-- _Rather than having the user do this through the Python command line,
  20 it might be easier to write a file that is going to do each step_ -->
  21
  22 *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*
  23
  24 If everything is set up correctly, you will be able to import the arvados SDK:
  25
  26 <pre>
  27 import arvados
  28 </pre>
  29
  30 ...and display your account information:
  31
  32 <pre>
  33 arvados.service.users().current().execute()
  34 </pre>
  35
  36 h3. More prerequisites
  37
  38 <pre>
  39 import re
  40 import json
  41 </pre>
  42
  43 h3. Find traits.
  44
  45 List traits containing the term "cancer":
  46
  47 <pre>
  48 for t in filter(lambda t: re.search('cancer', t['name']),
  49                 arvados.service.traits().list(limit=1000).execute()['items']):
  50   print t['uuid'], t['name']
  51
  52 </pre>
  53
  54 <!-- _Should break this down into steps instead of being clever and making
  55 it a python one-liner_ -->
  56
  57 %(darr)&darr;%
  58
  59 <pre>
  60 ...
  61 qr1hi-q1cn2-8q57g2diohwnzm0 Cervical cancer
  62 qr1hi-q1cn2-vqp4243janpjbyj Breast cancer
  63 qr1hi-q1cn2-v6usijujcpwqrn1 Non-melanoma skin cancer
  64 ...
  65 </pre>
  66
  67 We will use the "Non-melanoma skin cancer" trait with uuid @qr1hi-q1cn2-v6usijujcpwqrn1@.
  68
  69 <pre>
  70 trait_uuid = 'qr1hi-q1cn2-v6usijujcpwqrn1'
  71 </pre>
  72
  73 h3. Find humans.
  74
  75 List humans who report this condition:
  76
  77 <pre>
  78 trait_links = arvados.service.links().list(limit=1000,where=json.dumps({
  79     'link_class': 'human_trait',
  80     'tail_kind': 'arvados#human',
  81     'head_uuid': trait_uuid
  82   })).execute()['items']
  83 </pre>
  84
  85 <!-- _Same comment, break this out and describe each step_ -->
  86
  87 The "tail_uuid" attribute of each of these Links refers to a Human.
  88
  89 <pre>
  90 map(lambda l: l['tail_uuid'], trait_links)
  91 </pre>
  92
  93 %(darr)&darr;%
  94
  95 <pre>
  96 [u'1h9kt-7a9it-c0uqa4kcdh29wdf', u'1h9kt-7a9it-x4tru6mn40hc6ah',
  97 u'1h9kt-7a9it-yqb8m5s9cpy88i8', u'1h9kt-7a9it-46sm75w200ngwny',
  98 u'1h9kt-7a9it-gx85a4tdkpzsg3w', u'1h9kt-7a9it-8cvlaa8909lgeo9',
  99 u'1h9kt-7a9it-as37qum2pq8vizb', u'1h9kt-7a9it-14fph66z2baqxb9',
 100 u'1h9kt-7a9it-e9zc7i4crmw3v69', u'1h9kt-7a9it-np7f35hlijlxdmt',
 101 u'1h9kt-7a9it-j9hqyjwbvo9cojn', u'1h9kt-7a9it-lqxdtm1gynmsv13',
 102 u'1h9kt-7a9it-zkhhxjfg2o22ywq', u'1h9kt-7a9it-nsjoxqd33lzldw9',
 103 u'1h9kt-7a9it-ytect4smzcgd4kg', u'1h9kt-7a9it-y6tl353b3jc4tos',
 104 u'1h9kt-7a9it-98f8qave4f8vbs5', u'1h9kt-7a9it-gd72sh15q0p4wq3',
 105 u'1h9kt-7a9it-zlx25dscak94q9h', u'1h9kt-7a9it-8gronw4rbgmim01',
 106 u'1h9kt-7a9it-wclfkjcb23tr5es', u'1h9kt-7a9it-rvp2qe7szfz4dy6',
 107 u'1h9kt-7a9it-50iffhmpzsktwjm', u'1h9kt-7a9it-ul412id5y31a5o8',
 108 u'1h9kt-7a9it-732kwkfzylmt4ik', u'1h9kt-7a9it-v9zqxegpblsbtai',
 109 u'1h9kt-7a9it-kmaraqduit1v5wd', u'1h9kt-7a9it-t1nwtlo1hru5vvq',
 110 u'1h9kt-7a9it-q3w6j9od4ibpoyl', u'1h9kt-7a9it-qz8vzkuuz97ezwv',
 111 u'1h9kt-7a9it-t1v8sjz6dm9jmjf', u'1h9kt-7a9it-qe8wrbyvuqs5jew']
 112 </pre>
 113
 114 h3. Find PGP IDs.
 115
 116 For now we don't need to look up the Human objects themselves.
 117
 118 As an aside, we will look up "identifier" links to find PGP-assigned participant identifiers:
 119
 120 <pre>
 121 human_uuids = map(lambda l: l['tail_uuid'], trait_links)
 122 pgpid_links = arvados.service.links().list(limit=1000,where=json.dumps({
 123     "link_class": "identifier",
 124     "head_uuid": human_uuids
 125   })).execute()['items']
 126 map(lambda l: l['name'], pgpid_links)
 127 </pre>
 128
 129 %(darr)&darr;%
 130
 131 <pre>
 132 [u'hu01024B', u'hu11603C', u'hu15402B', u'hu174334', u'hu1BD549', u'hu237A50',
 133  u'hu34A921', u'hu397733', u'hu414115', u'hu43860C', u'hu474789', u'hu553620',
 134  u'hu56B3B6', u'hu5917F3', u'hu599905', u'hu5E55F5', u'hu602487', u'hu633787',
 135  u'hu68F245', u'hu6C3F34', u'hu7260DD', u'hu7A2F1D', u'hu94040B', u'hu9E356F',
 136  u'huAB8707', u'huB1FD55', u'huB4883B', u'huD09050', u'huD09534', u'huD3A569',
 137  u'huDF04CC', u'huE2E371']
 138 </pre>
 139
 140 These PGP IDs let us find public profiles:
 141
 142 * "https://my.personalgenomes.org/profile/huE2E371":https://my.personalgenomes.org/profile/huE2E371
 143 * "https://my.personalgenomes.org/profile/huDF04CC":https://my.personalgenomes.org/profile/huDF04CC
 144 * ...
 145
 146 h3. Find data.
 147
 148 Find Collections that were provided by these Humans.
 149
 150 <pre>
 151 provenance_links = arvados.service.links().list(where=json.dumps({
 152     "link_class": "provenance",
 153     "name": "provided",
 154     "tail_uuid": human_uuids
 155   })).execute()['items']
 156 collection_uuids = map(lambda l: l['head_uuid'], provenance_links)
 157
 158 # build map of human uuid -> PGP ID
 159 pgpid = {}
 160 for pgpid_link in pgpid_links:
 161   pgpid[pgpid_link['head_uuid']] = pgpid_link['name']
 162
 163 # build map of collection uuid -> PGP ID
 164 for p_link in provenance_links:
 165   pgpid[p_link['head_uuid']] = pgpid[p_link['tail_uuid']]
 166
 167 # get details (e.g., list of files) of each collection
 168 collections = arvados.service.collections().list(where=json.dumps({
 169     "uuid": collection_uuids
 170   })).execute()['items']
 171
 172 # print PGP public profile links with file locators
 173 for c in collections:
 174   for f in c['files']:
 175     print "https://my.personalgenomes.org/profile/%s %s %s%s" % (pgpid[c['uuid']], c['uuid'], ('' if f[0] == '.' else f[0]+'/'), f[1])
 176
 177 </pre>
 178
 179 %(darr)&darr;%
 180
 181 <pre>
 182 https://my.personalgenomes.org/profile/hu43860C a58dca7609fa84c8c38a7e926a97b2fc+302+K@qr1hi var-GS00253-DNA_A01_200_37-ASM.tsv.bz2
 183 https://my.personalgenomes.org/profile/huB1FD55 ea30eb9e46eedf7f05ed6e348c2baf5d+291+K@qr1hi var-GS000010320-ASM.tsv.bz2
 184 https://my.personalgenomes.org/profile/huDF04CC 4ab0df8f22f595d1747a22c476c05873+242+K@qr1hi var-GS000010427-ASM.tsv.bz2
 185 https://my.personalgenomes.org/profile/hu7A2F1D 756d0ada29b376140f64e7abfe6aa0e7+242+K@qr1hi var-GS000014566-ASM.tsv.bz2
 186 https://my.personalgenomes.org/profile/hu553620 7ed4e425bb1c7cc18387cbd9388181df+242+K@qr1hi var-GS000015272-ASM.tsv.bz2
 187 https://my.personalgenomes.org/profile/huD09534 542112e210daff30dd3cfea4801a9f2f+242+K@qr1hi var-GS000016374-ASM.tsv.bz2
 188 https://my.personalgenomes.org/profile/hu599905 33a9f3842b01ea3fdf27cc582f5ea2af+242+K@qr1hi var-GS000016015-ASM.tsv.bz2
 189 https://my.personalgenomes.org/profile/hu599905 d6e2e57cd60ba5979006d0b03e45e726+81+K@qr1hi Witch_results.zip
 190 https://my.personalgenomes.org/profile/hu553620 ea4f2d325592a1272f989d141a917fdd+85+K@qr1hi Devenwood_results.zip
 191 https://my.personalgenomes.org/profile/hu7A2F1D 4580f6620bb15b25b18373766e14e4a7+85+K@qr1hi Innkeeper_results.zip
 192 https://my.personalgenomes.org/profile/huD09534 fee37be9440b912eb90f5e779f272416+82+K@qr1hi Hallet_results.zip
 193 </pre>
 194
 195 h3. Search for a variant.
 196
 197 Look for variant rs1126809 in each of the "var" files (these contain variant calls from WGS data).
 198
 199 <pre>
 200 job = {}
 201 for c in collections:
 202   if [] != filter(lambda f: re.search('^var-.*\.tsv\.bz2', f[1]), c['files']):
 203     job[c['uuid']] = arvados.service.jobs().create(body={
 204       'script': 'grep',
 205       'script_parameters': {'input': c['uuid'], 'pattern': "rs1126809\\b"},
 206       'script_version': 'e7aeb42'
 207     }).execute()
 208     print "%s %s" % (pgpid[c['uuid']], job[c['uuid']]['uuid'])
 209
 210 </pre>
 211
 212 &darr;
 213
 214 <pre>
 215 hu43860C qr1hi-8i9sb-wyqq2eji4ehiwkq
 216 huB1FD55 qr1hi-8i9sb-ep68uf0jkj3je7q
 217 huDF04CC qr1hi-8i9sb-4ts4cvx6mbtcrsk
 218 hu7A2F1D qr1hi-8i9sb-5lkiu9sh7vdgven
 219 hu553620 qr1hi-8i9sb-nu4p6hjmziic022
 220 huD09534 qr1hi-8i9sb-bt9389e9g3ff0m1
 221 hu599905 qr1hi-8i9sb-ocg0i8r75luvke3
 222 </pre>
 223
 224 Monitor job progress by refreshing the Jobs page in Workbench, or by using the API:
 225
 226 <pre>
 227 map(lambda j: arvados.service.jobs().get(uuid=j['uuid']).execute()['success'], job.values())
 228 </pre>
 229
 230 %(darr)&darr;%
 231
 232 <pre>
 233 [True, True, True, True, True, True, True]
 234 </pre>
 235
 236 (Unfinished jobs will appear as None, failed jobs as False, and completed jobs as True.)
 237
 238 After the jobs have completed, check output file sizes.
 239
 240 <pre>
 241 for collection_uuid in job:
 242   job_uuid = job[collection_uuid]['uuid']
 243   job_output = arvados.service.jobs().get(uuid=job_uuid).execute()['output']
 244   output_files = arvados.service.collections().get(uuid=job_output).execute()['files']
 245   print "%s %3d %s" % (pgpid[collection_uuid], output_files[0][2], job_output)
 246
 247 </pre>
 248
 249 %(darr)&darr;%
 250
 251 <pre>
 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 </pre>
 260
 261 Thus, of the 7 WGS results available for PGP participants reporting non-melanoma skin cancer, 4 include the rs1126809 / TYR-R402Q variant.