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¤t=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
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
> >http://i30.photobucket.com/albums/c334/Tarek325/totalbiomassN01IMP100...
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
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