Best way to extract scan data for a specific genomic range

[Saikat Bandyopadhyay]

Hi,

I have a large scan object consisting of several phenotypes. 

I am trying to extract the LOD values for some specific phenotypes (lodcolumns) within a specific genomic range (e.g. the confidence interval of the QTL peak for the phenotype with 1 Mb margin on both upstream and downstream).

I could come up with a way to get the markers within the genomic range of interest using the physical map, and then use the markers to extract the LOD data from the scan object for the phenotype of interest.

Is this the correct approach?

For some reason I fell my current process is a little bit cumbersome, are there any other potentially more direct ways of achieving the same?

My intention is to use the extracted data in Gviz plots.

Thank you so much for your help in advance.

Best,

Saikat

[Karl Broman]

Yes, I would grab the names of markers/pseudomarkers in the interval of interest and use those to grab the portion of the scan1 output.

The function find_marker has an argument interval that you can use to pull out the names of markers within a genomic interval.

Here's an example:

library(qtl2)
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))

# insert pseudomarkers every 1 cM
gmap <- insert_pseudomarkers(iron$gmap, step=1)

# find corresponding physical locations
pmap <- interp_map(gmap, iron$gmap, iron$pmap)

# genome scan
probs <- calc_genoprob(iron, gmap, error_prob=0.002)
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
Xcovar <- get_x_covar(iron)
out <- scan1(probs, pheno, addcovar=covar, Xcovar=Xcovar)

# lod interval for "spleen" phenotype on chr 9
lodint_c9 <- lod_int(out, pmap, lodcolumn="spleen", chr=9)

# find corresponding marker positions, adding 1 Mbp on each end
interval <- lodint_c9[c(1,3)] + c(-1, 1)
markers_in_interval <- find_marker(pmap, 9, interval=interval)

# pull out that part of lod scores
out[markers_in_interval,]

karl

[Karl Broman]

Oops in the line calling scan1(), it should have iron$pheno in place of pheno

out <- scan1(probs, iron$pheno, addcovar=covar, Xcovar=Xcovar)