1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
|
<!-- header fragment for html documentation -->
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2 Final//EN">
<HTML>
<HEAD>
<META NAME="description" CONTENT="Estimation of population parameters using genetic data using a maximum likelihood approach with Metropolis-Hastings Monte Carlo Markov chain importance sampling">
<META NAME="keywords" CONTENT="MCMC, Markov chain, Monte Carlo, Metropolis-Hastings, population, parameters, migration rate, population size, recombination rate, maximum likelihood">
<TITLE>LAMARC Documentation: Mapping</title>
</HEAD>
<BODY BGCOLOR="#FFFFFF">
<!-- coalescent, coalescence, Metropolis-Hastings, Markov chain Monte Carlo
simulation, migration rate, effective population size, recombination rate,
maximum likelihood, trait mapping -->
(<A HREF="parallel.html">Previous</A> | <A HREF="index.html">Contents</A> |
<A HREF="compiling.html">Next</A>)
<H2>Mapping Traits Using LAMARC</H2>
<P>As of LAMARC version 2.1, we now have the capacity to perform fine-scale
mapping of traits and diseases that have been localized to a single stretch
of DNA (in LAMARC terms, to one genomic region). The way it works is that
as LAMARC searches among trees, it calculates the liklihood of the observed
trait data being created at each site in the sequence. Assuming that each
site has an equal chance of being the true site, it then converts these
liklihoods into probabilities, an averages the resulting values over all the
trees LAMARC collects.
</P>
<P>We have one submitted paper and another in the works (as of early 2007)
that take this approach. Until those appear, more detailed theory can be
provided upon request.
<P>To use LAMARC as a mapping tool, you'll have to collect the data, convert
it into a LAMARC input file, set up the parameters for your search, and run
the analysis. We'll go through each step in succession.
During our discussion, we'll make use of example data file and converter
command for a simple trait mapping case. Included are
<a href="trait_mapping/traitData.mig">sample trait data</a>
(in migrate format),
and a <a href="trait_mapping/traitCmd.html">sample converter command</a> file.
(actual xml is <a href="trait_mapping/traitCmd.xml">here</a>.)
<UL>
<LI><A HREF="mapping.html#collect">Collecting Data</A></LI>
<LI> <A HREF="mapping.html#convert">Converting Data</a></LI>
<LI> <A HREF="mapping.html#move">Moving to LAMARC</a></LI>
<LI> <A HREF="mapping.html#setup">Setting up parameters</a></LI>
<LI> <A HREF="mapping.html#analyze">Analyzing your results</a></LI>
</UL>
<h3><A NAME="collect">Collecting Data</h3>
<P>One of the assumptions of LAMARC is that the samples are collected at
random from the population. This affects the types of trees one would
expect to see that relate your samples to one another. However, researchers
that study a particular disease often have one set of samples from affected
individuals, and another set of samples from unaffected individuals. This
introduces an ascertainment bias that makes the true tree less likely to be
produced than it should. We have not yet investigated how strongly this
effects the accuracy of mapping results, but it is important to keep in
mind.
<P>The ideal trait for a LAMARC analysis would be one that is ubiquitously
diverse in the studied population (like eye color in humans). A random
sample from the population as a whole could then be classified into
different phenotypic groups. More information on the ascertainment question
will be provided as our own experiments proceed.
<H3><A NAME="convert">Converting Data</H3>
<P>Once you have collected data (of any type--marker SNPs, full sequence
information, or microsatellites), you will need to get the following
information into lamarc:
<ul>
<li> trait and allele names,
<li> the specific genomic region in which the trait is to be searched, and
<li> the correspondence between data samples, individuals, and the
phenotypes they express.
</ul>
<p>At this writing, most of this information can only be entered
into the lamarc file by writing a
<a href="converter_cmd.html">converter command file</a>.
You should refer to that document for the complete syntax of that
file. Below is a brief overview of the mapping-specific components
of that file.
<H4><a name="trait-info-defs">Defining Trait and Allele Names</a></H4>
Below is the section of our example
<a href="converter_cmd.html">converter command file</a>
which defines our trait and its alleles.
<pre>
<lamarc-converter-cmd>
...
<traits>
<trait-info>
<name>funny-nose</name>
<allele>normal</allele>
<allele>affected</allele>
</trait-info>
<phenotype>
...
</phenotype>
</traits>
...
</lamarc-converter-cmd>
</pre>
<P>This would define the 'funny-nose' trait, with its two alleles, normal
and affected.</P>
<H4><a name="phenotype-defs">Defining Phenotypic Groups</a></H4>
<p>Phenotypic groups are collections of individuals who all
display the same phenotype. This requires that one know or have a
reasonable model for the mechanism and penetrance of the trait in question.
In a simple dominant/recessive case, individuals displaying the dominant
phenotype could be any of AA, Aa, or aA, while individuals displaying the
recessive phenotype all must be aa.</P>
<p>If your data set includes many individuals with the same phenotype, the
easiest way to specify them is to define <pheontype> tags in the
converter command file. Below are the definitions for our simple example.
<pre>
<lamarc-converter-cmd>
...
<traits>
...
<trait-info>
...
</trait-info>
<phenotype>
<name>straight</name>
<genotype-resolutions>
<trait>funny-nose</trait>
<haplotypes>
<penetrance> 1.0 </penetrance>
<alleles> normal normal </alleles>
</haplotypes>
</genotype-resolutions>
</phenotype>
<phenotype>
<name>bent</name>
<genotype-resolutions>
<trait>funny-nose</trait>
<haplotypes>
<penetrance> 1.0 </penetrance>
<alleles> affected normal </alleles>
</haplotypes>
</genotype-resolutions>
<haplotypes>
<penetrance> 1.0 </penetrance>
<alleles> normal affected </alleles>
</haplotypes>
</phenotype>
<phenotype>
<name>broken</name>
<genotype-resolutions>
<trait>funny-nose</trait>
<haplotypes>
<penetrance> 1.0 </penetrance>
<alleles> affected affected </alleles>
</haplotypes>
</genotype-resolutions>
</phenotype>
</traits>
...
</lamarc-converter-cmd>
</pre>
<P>Note that the 'penetrance' in all instances is the chance that an
individual with those alleles displaying the observed phenotype. If
heterozygotes were 80% likely to display the same phenotype as all
normal/normal individuals, and 20% likely to display the same phenotype as
all broken/broken individuals, there would be only 2 phenotypes to define
and they would be as follows:
<pre>
<phenotype>
<name>straight</name>
<genotype-resolutions>
<trait>funny-nose</trait>
<haplotypes>
<penetrance> 1.0 </penetrance>
<alleles> normal normal </alleles>
</haplotypes>
<haplotypes>
<penetrance> 0.8 </penetrance>
<alleles> normal affected </alleles>
</haplotypes>
<haplotypes>
<penetrance> 0.8 </penetrance>
<alleles> affected normal </alleles>
</haplotypes>
</genotype-resolutions>
</phenotype>
<phenotype>
<name>broken</name>
<genotype-resolutions>
<trait>funny-nose</trait>
<haplotypes>
<penetrance> 1.0 </penetrance>
<alleles> affected affected </alleles>
</haplotypes>
<haplotypes>
<penetrance> 0.2 </penetrance>
<alleles> normal affected </alleles>
</haplotypes>
<haplotypes>
<penetrance> 0.2 </penetrance>
<alleles> affected normal </alleles>
</haplotypes>
</genotype-resolutions>
</phenotype>
</pre>
<P>If additional information is available about an individual that provides a
clue to their genotype, that information can also be included as a part of
the phenotype. For example, an individual with a dominant phenotype who has
a recessive parent or child must be a heterozygote, and should be classified
separately from other individuals displaying a dominant phenotype.</P>
<h4><a name="pheno-to-ind">Assigning phenotypes to individuals</a></h4>
In order to identify which data samples correspond to an individual
and what the phenotype of that individual is, use the <individual>
tag.
<pre>
<lamarc-converter-cmd>
...
<individuals>
<individual>
<name>ind_a</name>
<sample><name>s0</name></sample>
<sample><name>s1</name></sample>
<has-phenotype>broken</has-phenotype>
</individual>
...
</individuals>
...
</lamarc-converter-cmd>
</pre>
<p>This indicates that the samples "s0" and "s1" belong
to the same individual and that that individual displays the
"broken" phenotype. (If your individuals are named in your data
input file but do not have sample names, as in migrate-format microsatellite
data, you do not need to provide extra sample names here.)</P>
<p>
It is also possible to explicitly specify genotype resolutions for
an individual instead of using the <tt><has-phenotype></tt>
tag.
See the section
<a href="converter_cmd.html#phase">
Specifying Relationships Between Individuals and Data Samples</a> in the
<a href="converter_cmd.html">Converter Command File</a> documentation.
</p>
<h4>Identifying a trait with a genomic region</h4>
<p>Finally, the trait mapper needs to know in which genomic region to search
for your trait. This is an additional tag inside the "region" tag.
<pre>
<lamarc-converter-cmd>
...
<regions>
<region>
<name>region1</name>
...
<trait-location>
<trait-name>funny-nose</trait-name>
</trait-location>
...
</region>
...
</regions>
</lamarc-converter-cmd>
</pre>
<P>Data from other regions may be included in your analysis, but will not
affect the results from your mapping analysis one way or another. If you
have two unlinked regions, either of which might contain the trait allele,
LAMARC will be able to tell you the most likely place within either region
where the trait is likely to map, but will not be able to provide any
information about whether one region fits the data better than the
other.</P>
<h3><A NAME="move">Moving to Lamarc</h3>
<P>At this point, you are done creating a converter input file, and you can
run the converter and output a LAMARC input file (<a
HREF="trait_mapping/lamarc-trait-input.html">lamarc-trait-input.xml</a> [actual xml is <a HREF="trait_mapping/lamarc-trait-input.xml">here]</a>).
Run the converter in batch mode using this command file with:</P>
<pre>
lam_conv -b -c traitCmd.xml
</pre>
<P>You should also be able to further modify your run within the converter
by leaving off the '-b' to run in interactive mode. Once you're done,
select 'Write Lamarc File' from the 'File' menu.</P>
<P>Then, start up lamarc, using your newly-minted file as the input.</P>
<h3><A NAME="setup">Setting up parameters</h3>
<P>The settings for mapping are under the Analysis menu ('A', from the main
menu), then 'A' again for the 'Trait Data analysis' menu. The next menu
shows a list of all the traits you are mapping; the settings for each can be
changed independently. Select the trait for which you want to change the
settings.</P>
<P>There are two things you can do from this menu. The first is to restrict
the range of possible sites for your trait. For some analyses, you might
have sequence data in the same general region as the area where you've
mapped your trait, but you know that the trait alleles are not there. In
this case, you can Add ('A') or Remove ('R') sites from the genomic region
from consideration as being a possible site for the trait you're mapping.
The numbering scheme used here is the region-wide numbering, as defined by
the 'map-position' numbers both in the converter and in the lamarc input
file. Also note that by default, there is no 'site 0' in this scheme: the
position to the left of site 1 is defined to be site -1. (This can be
changed in the lamarc input file XML tag <tt><A
HREF="xmlinput.html#options">convert-output-to-eliminate-zero</a></tt>.</P>
<P>The second option here is to select the type of mapping analysis you wish to
perform, a 'floating' ('F') or a 'jumping' ('J') analysis.</P>
<P>The 'floating' analysis (the default) will not affect the search through
tree-space at all. Instead, as trees are collected, each is analyzed, and
the likelihood that the observed trait data could have been produced at each site is
calculated and stored. This can be a somewhat time-consuming process, but
the results are more robust. The final result is a complete analysis, for
every tree collected in the final chain, of the relative likelihoods of each
site in the genomic region under analysis.</P>
<P>In the 'jumping' analysis, the trait alleles are actually placed at a
particular site. A new arranger is turned on that moves the trait
alleles from site to site based on the relative likelihoods of the data
being created at each site. This arranger is a 'Gibbs' arranger, in that it
calculates the likelihood at each site, and then chooses one according to
their relative probabilities. In this way, it always accepts its new
choice, without regard to where it used to be. The results of this analysis
are an average of where the trait alleles were placed during the run.</P>
<P>Based on our preliminary analyses, we recommend the 'floating' analysis for
most studies. This analysis is more thorough, and while the analysis takes
somewhat longer while it is running, it only needs to run during the final
chain, so on balance, the time difference is relatively small. In addition,
by never incorporating the trait data into tree acceptance or rejection, the
trait data itself does not influence the tree and 'lock' it into place by
shoving the trait into a fairly reasonable position and then modifying the tree to fit
the data. With known data, this is exactly what you want, but in this case,
we want to know where on the sequence the allele resides, without unduly
biasing the results.</P>
<P>Still, the 'jumping' analysis has its uses. Since its arranger is on at all
times, it provides updates on the current state of the mapper at the end of
each chain. And while normally we believe that not incorporating trait
data into tree acceptance and rejection is good for mapping, there may be
some cases (perhaps with sparse markers) where using the data from the trait
may lead to finding trees with overall better likelihood at the correct
position, with a corresponding lower likelihood at incorrect positions. In
this case, the range of sites where the trait might be located would narrow
(a good thing, when mapping).</P>
<P>When performing a Jumping analysis, two new rearrangers are created, and the
relative amount of time spent on each can be set from the Rearrangers menu
('S', then 'R' from the main menu.) The basic arranger here is the Trait
Location rearranger ('L'), which accomplishes the 'jumping' part of the
analysis by moving the trait alleles to various positions in the sequence.
The other arranger you might find here is the Trait haplotypes rearranger
('M', for no good reason), which will rearrange the unknown haplotypes (if
any) of your trait data. This rearranger will take the various
genotype resolutions you defined for your phenotypes and individuals, and
swap among them in proportion to their relative likelihoods. If there are a
lot of individuals in your data with unknown haplotypes, this can be another
reason to use the 'jumping' analysis, since this analysis searches among
possible haplotype resolutions, while the 'floating' analysis sums over all
possible combinations of haplotype resolutions for each analyzed tree, which
can be very slow.</P>
<h3><A NAME="analyze">Analyzing your results</h3>
<P>Running LAMARC in interactive mode (i.e. non-batch mode) will result in
output to the screen that looks something like this:</P>
<pre>
15:44:03 Most likely site(s) for funny-nose: 919:923. Relative data likelihood = 0.0025891
The top 5% of all sites in this region: 919:936, 946:947
The top 50% of all sites in this region: 678:709, 776:803, 811:953
The top 95% of all sites in this region: 1:212, 214:215, 217:220, 641:953
You have a total of 531 sites in your 95% range.
</pre>
<P>These <A HREF="trait_mapping/outfile.txt">results</a> were for a mock
data set for the 'funny-nose' trait, in a region 1000 base pairs long. The
top scorers were sites 837:841 although each site in this range had only a
0.22% chance of actually being the correct one. The last line informs us
that 566 sites are included at 95% confidence, so about half the sequence
was excluded. The true site for this trait happened to be site 905, which
was included in the top 50% of all sites. This was the result of a default
run, which is on the short side; it may be the case that we'd get better
results from spending more time on a longer run of LAMARC.</P>
<P>Note that these results are disjointed--the most likely sites are not
guaranteed to be next to one another. This can be due to both random chance
and shared heritage. Random chance can produce two different historical
pathways which are both about equally good explanations of the trait.
Shared genetic heritage can sometimes mean that distant segments of DNA
might share a crucial bit of genetic history, while the intervening segment
experienced a widely divergent history.</P>
<P>More detailed results are available in the mapping output file for each
trait. For our example, the default name for this file is
<a HREF="trait_mapping/mapfile_funny-nose.txt">mapfile_funny-nose.txt</a>.
There, you will find a list like this:</P>
<pre>
Site Data likelihood
1 0.00090735
2 0.00090871
3 0.00090877
4 0.00091575
5 0.00091603
6 0.00091603
7 0.00091603
8 0.00091756
[...]
</pre>
<P>This column of numbers should add up to 1.0, and indicates the percent
chance that your trait maps to each site (so, for the above, a 0.069191%
chance of being at site 1, a 0.069390% chance for site 2, and so on.) The
list can be imported into a graphing utility or spreadsheet program to
better visualize your results.</p>
<P>(<A HREF="parallel.html">Previous</A> | <A HREF="index.html">Contents</A> |
<A HREF="compiling.html">Next</A>)</P>
<!--
//$Id: mapping.html,v 1.23 2016/04/19 21:01:32 lpsmith Exp $
-->
</BODY>
</HTML>
|
