1.5 LOD support intervals are too wide-what am I doing wrong?

661 views
Skip to first unread message

Tarek Elnaccash

unread,
Feb 2, 2010, 3:27:13 PM2/2/10
to R/qtl discussion
I've tried to determine the location of 1.5 LOD support intervals
around my QTL using two methods and they don't match. The first
method was to use lodint:

lodint(out.lcSG1d100.acovmat.imp,drop=1.5,lodcolumn=1,chr=2)

chr pos RDtotg01
AD.156C 2 0 0.777195
CH.65C 2 30 2.212948
EC.235L-Col/247C 2 71 1.866133

, which tells me the support interval around this QTL is from 0-70 cM
on chromosome two.

This didn't seem to match the plot which can be viewed here: <a
href="http://s30.photobucket.com/albums/c334/Tarek325/?
action=view&current=totalbiomassN01IMP100mat.jpg" target="_blank"><img
src="http://i30.photobucket.com/albums/c334/Tarek325/
totalbiomassN01IMP100mat.jpg" border="0" alt="Photobucket"></a>

(if the html link above doesn't work here is the direct link:
http://i30.photobucket.com/albums/c334/Tarek325/totalbiomassN01IMP100mat.jpg
)

A 1.5 LOD drop from the peak should be located at positions which have
a 0.713 LOD score or less. You can see from the picture that this
would be a smaller interval than the one above at 0-70 cM.

I then exported my scanone data using write.table:

write.table(out.lcSG1d100.acovmat.imp,
file="out.lcSG1d100.acovmat.imp.txt",sep=" ")

and the data in that output suggest that the support interval is from
22-42 cM instead of the 0-71 cM. I've pasted this output below.

Why don't these two agree and how can I make LODINT produce the
correct interval boundaries?

thanks,

Tarek


Data from write.table for chromosome 2 (marker, chromosome, position,
LOD) :

AD.156C 2 0 0.77719462
c2.loc1 2 1 0.729647281
c2.loc2 2 2 0.684397164
c2.loc3 2 3 0.636867981
c2.loc4 2 4 0.580117961
c2.loc5 2 5 0.529064689
c2.loc6 2 6 0.46718347
BF.325L 2 7 0.414588429
c2.loc8 2 8 0.4667815
c2.loc9 2 9 0.543176571
c2.loc10 2 10 0.616893748
c2.loc11 2 11 0.689375761
c2.loc12 2 12 0.739424512
c2.loc13 2 13 0.790611257
GH.580L 2 14 0.830255686
DF.225L 2 15 0.948339407
BF.226C/BH.58L 2 16 1.244207383
EC.495C-Col 2 17 1.166082496
BH.460L 2 17.1 1.166082496
c2.loc18 2 18 1.020975196
FD.81L 2 19 0.756464884
BF.105C 2 20 0.86430616
c2.loc21 2 21 0.803636857
CH.284C 2 22 0.622087224
c2.loc23 2 23 0.735259225
c2.loc24 2 24 0.85913975
c2.loc25 2 25 1.010600481
c2.loc26 2 26 1.151798211
FD.222L-Col 2 27 1.280301824
c2.loc28 2 28 1.427848727
CD.246L 2 29 1.578176272
CH.65C 2 30 2.212947598
c2.loc31 2 31 2.025395202
c2.loc32 2 32 1.825944826
c2.loc33 2 33 1.581208696
CH.1500C 2 34 1.382551665
BF.221L 2 35 1.323024071
c2.loc36 2 36 1.307786249
c2.loc37 2 37 1.239925696
c2.loc38 2 38 1.05926906
c2.loc39 2 39 0.907871742
FD.85C 2 40 0.825079372
c2.loc41 2 41 0.665622389
c2.loc42 2 42 0.393612847
GB.150L-Col 2 43 0.27517926
c2.loc44 2 44 0.467761669
FD.150C 2 45 0.689451673
c2.loc46 2 46 0.757343671
GD.460L-Col 2 47 0.81239517
c2.loc48 2 48 1.04764787
CH.145L-Col/150C 2 49 1.188805034
c2.loc50 2 50 1.129491903
c2.loc51 2 51 1.059133869
c2.loc52 2 52 1.003829835
BH.195L-Col 2 53 0.946511605
GD.298C 2 54 0.935983334
c2.loc55 2 55 0.908297746
GH.247L 2 56 0.891839949
c2.loc57 2 57 1.122085443
c2.loc58 2 58 1.25220947
BH.120L-Col 2 59 1.372749712
c2.loc60 2 60 1.397826111
c2.loc61 2 61 1.427431009
c2.loc62 2 62 1.427357865
DF.140C 2 63 1.445231925
c2.loc64 2 64 1.567013336
c2.loc65 2 65 1.707448056
c2.loc66 2 66 1.855858523
c2.loc67 2 67 1.900587261
c2.loc68 2 68 2.083107121
c2.loc69 2 69 2.050827611
c2.loc70 2 70 2.010399981
EC.235L-Col/247C 2 71 1.866132803

Karl Broman

unread,
Feb 2, 2010, 5:07:48 PM2/2/10
to rqtl...@googlegroups.com
In lodint(), we take all positions with LOD within 1.5 of the maximum
for that chromosome, and merge into a single contiguous interval, and
so for your results, that second peak leads to the long interval.

If you want it to ignore the 2nd peak, you'd need to replace LOD
scores with NAs at some point along the chromosome.

Karl

On Feb 2, 2010, at 3:27 PM, Tarek Elnaccash
<tarek.e...@gmail.com> wrote:

> --
> You received this message because you are subscribed to the Google
> Groups "R/qtl discussion" group.
> To post to this group, send email to rqtl...@googlegroups.com.
> To unsubscribe from this group, send email to rqtl-disc+...@googlegroups.com
> .
> For more options, visit this group at http://groups.google.com/group/rqtl-disc?hl=en
> .
>

Julin Maloof

unread,
Feb 3, 2010, 9:58:32 AM2/3/10
to R/qtl discussion
This brings up a philosophical issue that I have always been curious
about. Many of the summary functions in rqtl are geared towards the
assumption of 1 QTL per chromosome (or at least only report the
largest single QTL per chromosome). Why is that? Do you feel that
realistically this is the most that can be confidently resolved?

Julin

> >http://i30.photobucket.com/albums/c334/Tarek325/totalbiomassN01IMP100...

Karl Broman

unread,
Feb 5, 2010, 1:15:41 AM2/5/10
to rqtl...@googlegroups.com, Julin Maloof

Dear Julin,

Interval mapping (as in scanone) involves the fitting of a series of single-QTL model. While the presence of multiple well-separated LOD peaks is suggestive of the presence of multiple QTL, the interval mapping results, on their own, are not sufficient to discriminate between the case of one vs two or more QTL. For that, it is best to fit multiple-QTL models, either with scantwo or makeqtl/fitqtl/refineqtl. And that's the reason that I've not put any effort into writing code to identify the bumps in scanone results.

As an example, consider the hyper data (Sugiyama et al. 2001). There are two peaks on chr 1 (and a big peak on chr 4). If one is willing to assume the presence of two QTL on chr 1, the best way to get interval estimates of their locations is to use makeqtl to create a QTL object, refineqtl to get improved estimates of the QTL locations, in the context of the multiple-QTL model, and then lodint.

data(hyper)
hyper <- sim.geno(hyper, step=2, n.draws=128)
qtl <- makeqtl(hyper, chr=c(1,1,4), pos=c(49, 80, 30))
rqtl <- refineqtl(hyper, qtl=qtl)
summary(fitqtl(hyper, qtl=rqtl))

lodint(rqtl, qtl.index=1) # 21.3-73.3
lodint(rqtl, qtl.index=2) # 64.5-91.3
lodint(rqtl, qtl.index=3) # 28.4-31.7

The performance of these LOD intervals has not been well studied, and they don't fully capture the uncertainty in the QTL locations. For example, the interval for the proximal QTL on chr 1 is derived assuming that the locations of the distal chr 1 QTL and the chr 4 QTL are known exactly.

regards,
karl

brassica

unread,
Feb 5, 2010, 10:35:21 AM2/5/10
to R/qtl discussion
Dear Karl,
On the two QTL on one chromesome and large intervals, I normally
consider two things and in my case it works
1. If two QTL peak intervals larger than 50cM with relative large LOD,
then I consider two QTL on that chromesome; if within 50cM, there are
two QTL peak then the interval mapping results might be not tricky.
2. About the LOD1.0/1.5/2.0, whatever you are using, if the LOD below
1, you cann't trust.
Is my considering right?
Best regards
Ping

Karl Broman

unread,
Feb 5, 2010, 10:41:15 AM2/5/10
to rqtl...@googlegroups.com, brassica
Dear Ping,

That sounds reasonable, but it's hard to form a hard-and-fast rule for the interpretation of multiple peaks in a LOD curve, and I still insist that one should look at the fit of the two-QTL model versus the one-QTL model.

And indeed, for LOD support intervals, you need the actual LOD peak to be large; otherwise the interval would just be the whole chromosome. And the interval doesn't have much meaning if there is not evidence for a QTL.

karl

Reply all
Reply to author
Forward
0 new messages