Bootstrapping in rqtl2

[Eden McQueen]

Hi Karl,

I am trying to run bootstraps by resampling my individuals with replacement in a back cross. I need the full lod by position output for each bootstrap.

In rqtl, I modified the existing scanoneboot function (subsetting the cross with resampled individuals and then running scanone) by just asking it to print the data that I needed at the end.

In rqtl2, I am having trouble doing something similar, because now my id column seems to be tied to a row names attribute, and therefore the individuals cannot repeat. I came up with a very inefficient workaround by generating new data files and arbitrarily renaming the id column prior to reading the cross but, given that we have over 14,000 markers, generating all of those new files takes a lot of time and therefore the lack of efficiency is pretty detrimental for us. 

I wanted to know whether there is a more efficient way to perform bootstrapping in rqtl2. Is there a workaround that would be better than what I am doing?

Thanks so much,

Eden 

[Karl Broman]

Yes, a bootstrap will be a bit of a pain, due to the system for ensuring that individuals are properly aligned.

You'd have to sample the set of individuals, do the subsetting, and then stick in a set of unique IDs for all individuals.

Here's a stab at it: scan1boot.R (on gist).

karl

[Eden McQueen]

Wow! Many thanks for the quick reply and for the script! 

I gave it a try with my data today, and also attempted it with the iron data. Right now running scan1boot returns a function (in the example, the object "z" is a function). However, the data manipulation you did is definitely more efficient than what I was trying to do, so this is going to be a much better solution for us once I can get it to work for me. If you spot something in the script that might be causing the issue, let me know!

Perhaps I should start a new thread, but I also wondered whether there are plans to implement mqm in rqtl2? 

Thanks again,

Eden

[Karl Broman]

The output should be a list, each component being the output of scan1() from one bootstrap replicate.

You could grab the positions with maximum LOD on chromosome 9 for the "spleen" phenotype like this:

    sapply(z, function(a) max(a, iron$gmap, lod="spleen", chr=9)$pos)

I don't have immediate plans to implement MQM. I do want to introduce multiple-QTL mapping methods, but I don't like the 1d scan with windowing in MQM and CIM.

karl

[Eden McQueen]

Hi Karl, 

Thank you again. Apologies for the delay in my response. This has been extremely helpful.

This is a little off topic for my original question, but since you mentioned your frustration with the current iteration of MQM in rqtl, may I ask a quick question about that? 

For our data the Multiple-QTL Mapping approach would really be better, as we believe we have quite a few QTL per chromosome in some cases. One thing that is odd, however, is that the MQM model that is generated by mqmscan does not seem to match up well with the lod profile in the plot output. Specifically, in some cases the model will indicate that there are two QTL next to each other, but the lod profile will look like one peak with a lod max somewhere between the two. Do you know why that may be happening? I have not been able to find any other examples like this.

All the best,

Eden  

[Karl Broman]

It would be best to start a separate line of discussion, and back at R/qtl discussion group if about mqm. And it’s hard to say anything without seeing the details. 

karl

On Jun 23, 2020, at 8:59 AM, Eden McQueen <eden.m...@gmail.com> wrote:



--
You received this message because you are subscribed to the Google Groups "R/qtl2 discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to rqtl2-disc+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/rqtl2-disc/6716cc91-d01f-4834-a477-7855d83575c3o%40googlegroups.com.