```{r setup, include=FALSE}
knitr::opts_chunk$set(results="hide", warning=FALSE, message=FALSE)
options(width=110, scipen=100)
```

This analysis serves to do three things:

- Shift the
  [Liu et al (2014)](https://doi.org/10.1534/genetics.114.161653)
  mouse genetic map so that 0 cM corresponds to 0 Mbp.

- Calculate a grid of markers across the genome.

- Get interpolated cM positions for the markers on the GigaMUGA array.


I'll use the [R/qtl2scan](https://github.com/rqtl/qtl2scan) package to
do interpolations,
[R/qtl2convert](https://github.com/rqtl/qtl2convert) to do some data
manipulation, and [R/broman](https:github.com/kbroman/broman) for a
few minor things. I'll also use a few other packages below
([rvest](https://github.com/hadley/rvest),
[AnnotationHub](http://bioconductor.org/packages/release/bioc/html/AnnotationHub.html),
and [devtools](https://github.com/hadley/devtools)).


```{r load_libraries}
library(qtl2geno)
library(qtl2scan)
library(qtl2convert)
library(broman)
```

### Download and prepare the maps

I first need to download the
[Liu et al (2014)](https://doi.org/10.1534/genetics.114.161653) maps
and do a bit of preparations.

I'll first ensure that I've got directories to contain the files I
download plus those I create.

```{r make_directories}
in_dir <- "files"
out_dir <- "results"
for(dir in c(in_dir, out_dir)) {
    if(!dir.exists(dir)) dir.create(dir)
}
```

Now I'll download and unzip the
[Liu et al (2014)](https://doi.org/10.1534/genetics.114.161653) map.

```{r download_liu_map}
url <- "http://cgd.jax.org/mousemapconverter/G2F1.anchor.maps.zip"
file <- file.path(in_dir, basename(url))
if(!file.exists(file)) {
    download.file(url, file)
}
unzip(file, exdir=in_dir)
```

I'll read in the sex-averaged, female, and male genetic maps.
I'm also going to give each row a name that is like `1_3036178` for
chromosome 1 at 3,036,178 bp.

```{r read_maps}
files <- c("avg.map.csv", "female.map.csv", "male.map.csv")
liu_map <- lapply(files, function(f) read.csv(file.path(in_dir, f)))
liu_map <- lapply(liu_map, function(m) { cbind(m, marker=paste(m$Chr, m$Pos, sep="_")) })
names(liu_map) <- lapply(strsplit(files, "\\."), "[", 1)
liu_map <- lapply(liu_map, function(m) {
    colnames(m)[1:2] <- tolower(colnames(m)[1:2])
    m$chr[m$chr==20] <- "X"
    m$chr <- factor(m$chr, c(1:19,"X"))
    m })
```

Strangely, they have different numbers of positions: the sex-averaged,
female, and male maps have `r add_commas(nrow(liu_map$avg))`
`r add_commas(nrow(liu_map$female))`, and `r add_commas(nrow(liu_map$male))`
positions, respectively.

And it turns out that the basepair positions are from build 37, so we
first need to convert them to build 38. (I will separately use
[liftOver](https://genome.ucsc.edu/cgi-bin/hgLiftOver) to convert the
positions to build 38.)

```{r save_liu_build37_to_file}
write.table(cbind(chr=paste0("chr", liu_map$avg$chr), pos=paste(liu_map$avg$pos,liu_map$avg$pos,sep="-")),
            file.path(in_dir, "liu_build37_bp.txt"), sep=":",
            row.names=FALSE, col.names=FALSE, quote=FALSE)
```

Now reading the liftOver results back in:

```{r read_liu_build38}
liu_build38 <- scan(file.path(out_dir, "liu_build38.txt"), what=character())
liu_build38 <- as.data.frame(matrix(unlist(strsplit(liu_build38, "[:\\-]")), ncol=3, byrow=TRUE),
                             stringsAsFactors=FALSE)[,1:2]
colnames(liu_build38) <- c("chr", "bp_build38")
liu_build38$chr <- factor(substr(liu_build38$chr, 4, nchar(liu_build38$chr)), levels=c(1:19,"X"))
liu_build38$bp_build38 <- as.numeric(liu_build38$bp_build38)
stopifnot(all(liu_build38$chr == liu_map$avg$chr))
liu_map$avg$build38 <- liu_build38$bp_build38
```

As it turns out, there's one pair of positions that are inverted in
build 38 vs build 37 (on chr 4), the middle two markers here:

```{r build38_inversion, results="markup"}
liu_map$avg[liu_map$avg$chr==4,][447:450,]
```

These markers are present in the male map but not the female map.
We'll omit the marker at build 37 position
`r add_commas(liu_map$avg$pos[liu_map$avg$chr==4][449])`.

```{r liu_drop_one_position}
drop <- 42229812
for(m in c("avg", "female", "male")) {
    map <- liu_map[[m]]
    map <- map[!(map$chr==4 & map$pos==drop),]
    liu_map[[m]] <- map
}
```

I next add the build 38 positions to the female and male maps.

```{r liu_add_build38}
for(m in c("female", "male")) {
    liu_map[[m]]$build38 <- liu_map$avg$build38[match(liu_map[[m]]$marker, liu_map$avg$marker)]
}
```

Now to align the maps and make a single data frame. I'm going to split
them into lists and interpolate positions in the female and male maps
so that they are all the same length. Note that the male map doesn't
include the X chromosome

```{r align_maps}
liu_pmap_list <- lapply(liu_map, map_df_to_list, chr_column="chr", pos_column="build38")
liu_gmap_list <- lapply(liu_map, map_df_to_list, chr_column="chr", pos_column="cM")
liu_gmap_list$female <- interp_map(liu_pmap_list$avg, liu_pmap_list$female, liu_gmap_list$female)
liu_gmap_list$male <- interp_map(liu_pmap_list$avg[1:19], liu_pmap_list$male, liu_gmap_list$male)
```

Now I'll combine everything into one data frame.

```{r combine_maps}
liu_map <- cbind(liu_map$avg[,c("marker", "chr", "build38", "pos", "cM")],
                 cM_female = map_list_to_df(liu_gmap_list$female, pos_column="cM_female")[,"cM_female"],
                 cM_male = c(map_list_to_df(liu_gmap_list$male, pos_column="cM_male")[,"cM_male"],
                             rep(NA, sum(liu_map$avg$chr=="X"))))
colnames(liu_map)[3:4] <- c("bp_build38", "bp_build37")
```


### Get the chromosome lengths

Now we want to anchor the maps at 0 cM. I also want to anchor them at
the telomeres. So I first need to get the chromosome lengths.
I'll start by grabbing the build 38 from the web,
at <https://www.ncbi.nlm.nih.gov/grc/mouse/data">.

I use [rvest](https://github.com/hadley/rvest) to scrape the
chromosome lengths from the web page.

```{r build38_chr_lengths}
library(rvest)
url <- "https://www.ncbi.nlm.nih.gov/grc/mouse/data"
file <- file.path(in_dir, "grc_mouse_genome_assembly.html")
if(!file.exists(file)) download.file(url, file)
z <- read_html(file)
tab <- html_nodes(z, css="table.ui-ncbigrid")
tab <- html_table(tab)[[1]]
mm10_L <- as.numeric(sapply(strsplit(tab[,2], ","), paste, collapse=""))[1:20]
names(mm10_L) <- c(1:19,"X")
```

Let's double-check these using
[AnnotationHub](http://bioconductor.org/packages/release/bioc/html/AnnotationHub.html).
I'm grabbing the cytobands which seem to be the smallest thing that
will give me full lengths.

```{r annotation_hub_lengths}
file <- file.path(in_dir, "ah51380.rds")
if(file.exists(file)) {
    zz <- readRDS(file)
} else {
    library(AnnotationHub)
    ah <- AnnotationHub()
    z <- query(ah, "mm10")
    zz <- as.data.frame(z[["AH53180"]]) # cytoband track
    saveRDS(zz, file)
}
mm10_Lb <- tapply(zz$end, zz$seqnames, max)
mm10_Lb <- setNames(mm10_Lb, substr(names(mm10_Lb), 4, nchar(mm10_Lb)))[1:20]
```

These **`r ifelse(all(mm10_L == mm10_Lb), "are", "are not")`**
the same as those from the web.

And are all longer than Liu map? `r all(mm10_L > tapply(liu_map$bp_build38, liu_map$chr, max))`.
`r ifelse(all(mm10_L > tapply(liu_map$bp_build38, liu_map$chr, max)), "Yay!", "")`


### Shift the maps

Okay, so now we've got the Liu et al. maps into build 38 and with the
sex-averaged, female, and male maps aligned. Next we want to add the 0
position plus the telomere and then determine the correponding cM
positions.

I'm going to do the split, interpolate, recombine, again.

First split:

```{r split_liu_map}
pmap <- map_df_to_list(liu_map, pos_column="bp_build38")
gmap_ave <- map_df_to_list(liu_map, pos_column="cM")
gmap_fem <- map_df_to_list(liu_map, pos_column="cM_female")
gmap_mal <- map_df_to_list(liu_map, pos_column="cM_male")
```

Now anchor the chromosomes at 0 and the telomere:

```{r anchor_liu_map}
pmap_anchored <- pmap
for(chr in seq_along(pmap)) {
    pmap_anchored[[chr]] <- c(0, pmap[[chr]], mm10_L[chr])
    names(pmap_anchored[[chr]])[c(1, length(pmap_anchored[[chr]]))] <-
        paste0(c("zero", "telo"), names(pmap[chr]))
}
```

Now I do the interpolation to get the cM positions for the zero and
telomere points.

```{r interpolate_liu_map}
gmap_ave <- interp_map(pmap_anchored, pmap, gmap_ave)
gmap_fem <- interp_map(pmap_anchored, pmap, gmap_fem)
gmap_mal <- interp_map(pmap_anchored, pmap, gmap_mal)
```

Now we shift the genetic maps so they start at 0.

```{r shift_liu_gmaps}
gmap_ave <- lapply(gmap_ave, function(a) a-min(a))
gmap_fem <- lapply(gmap_fem, function(a) a-min(a))
gmap_mal <- lapply(gmap_mal, function(a) a-min(a))
```

Finally, we bring them back together into a data frame.

```{r shifted_maps_to_df}
liu_map <- cbind(map_list_to_df(pmap_anchored, pos_column="bp_build38"),
                 cM=map_list_to_df(gmap_ave, pos_column="cM")[,"cM"],
                 cM_female=map_list_to_df(gmap_ave, pos_column="cM_female")[,"cM_female"],
                 cM_male=map_list_to_df(gmap_ave, pos_column="cM_male")[,"cM_male"])
```

We also need to change the position names to use the build38
positions.

```{r shifted_map_fix_markers}
wh <- grep("_", liu_map$marker)
rownames(liu_map)[wh] <- liu_map$marker[wh] <- paste(liu_map$chr[wh], liu_map$bp_build38[wh], sep="_")
```

Let's grab the physical and genetic maps as lists, for later use.

```{r liu_map_to_lists}
liu_pmap <- map_df_to_list(liu_map, pos_column="bp_build38")
liu_gmap <- map_df_to_list(liu_map, pos_column="cM")
```


### Interpolate GigaMUGA cM

Let's now get the interpolated cM positions for the GigaMUGA markers.
We first download the SNP map file and load it into R.

```{r download_gigamuga_map}
url <- "ftp://ftp.jax.org/MUGA/GM_snps.Rdata"
file <- file.path(in_dir, basename(url))
if(!file.exists(file)) {
    download.file(url, file)
}
load(file)
```

Now we split the physical positions into a list, first subsetting to
the autosomes and the X chromosome. Also, we convert the Mbp positions
to basepairs.

```{r subset_GM_snps}
GM_pmap <- map_df_to_list(GM_snps[GM_snps$chr %in% c(1:19,"X"),], pos_column="pos")
GM_pmap <- lapply(GM_pmap, function(a) round(a*1e6))
```

And now let's get interplated cM positions using the Liu et al. map.

```{r interpolate_GM}
GM_gmap <- interp_map(GM_pmap, liu_pmap, liu_gmap)
```

Finally, we paste the cM positions back into the `GM_snps` object.

```{r paste_in_GM_gmap}
GM_snps$cM[GM_snps$chr %in% c(1:19,"X")] <- map_list_to_df(GM_gmap, pos_column="cM")[,"cM"]
```

I'll also make the `chr` column a factor, with levels `1`, `2`, ...,
`19`, `X`, `Y`, and `M`.

```{r GM_snps_chr_make_factor}
GM_snps$chr <- factor(GM_snps$chr, levels=c(1:19,"X", "Y", "M"))
```

### Interpolate MegaMUGA cM

Let's do the same thing with the MegaMUGA markers.
We first download the SNP map file and load it into R.

```{r download_megamuga_map}
url <- "ftp://ftp.jax.org/MUGA/MM_snps.Rdata"
file <- file.path(in_dir, basename(url))
if(!file.exists(file)) {
    download.file(url, file)
}
load(file)
```

Note that the `MM_snps` object has markers on chr 1-19, X, Y, P, and
22. No missing values in the `chr` column, and no notations for the
mitochondria. But there are missing values in the `pos` column.

We split the physical positions into a list, first subsetting to
the autosomes and the X chromosome (and omitting snps missing
physical positions). Also, we convert the Mbp positions
to basepairs.

```{r subset_MM_snps}
MM_keep <- !is.na(MM_snps$pos) & MM_snps$chr %in% c(1:19,"X")
MM_pmap <- map_df_to_list(MM_snps[MM_keep,], pos_column="pos")
MM_pmap <- lapply(MM_pmap, function(a) round(a*1e6))
```

And now let's get interplated cM positions using the Liu et al. map.

```{r interpolate_MM}
MM_gmap <- interp_map(MM_pmap, liu_pmap, liu_gmap)
```

Finally, we paste the cM positions back into the `MM_snps` object.

```{r paste_in_MM_gmap}
MM_snps$cM[MM_keep] <- map_list_to_df(MM_gmap, pos_column="cM")[,"cM"]
```



### Genome grid in cM

Let's now calculate a grid along the genetic map, and find the corresponding
bp positions. On each chromosome, I'll have the grid start at **3
Mbp** and go to the telomere. (Traditionally, the mouse chromosomes
sequence all start at 3 Mbp, leaving space for the centromere.)

So first we need to add a position at 3 Mb for each chromosome and
interpolate to get the corresponding cM location.

```{r liu_add_3Mb}
liu_pmap_3Mbp <- lapply(liu_pmap, function(a) sort(c(a, start=3e6)))
liu_gmap_3Mbp <- interp_map(liu_pmap_3Mbp, liu_pmap, liu_gmap)
```

Now we want to form a grid starting at the second position on each
chromosome and ending at the telomere. KB Choi's 64k grid is a nice
round number, but the inter-marker distances are slightly different
between chromosomes. I'd be inclined to go with a 0.02 cM grid, which
would give
`r add_commas(sum(sapply(liu_gmap_3Mbp, function(a) length(seq(a[2], a[length(a)], by=0.02)))))`
positions, or a 0.01 cM grid, which would give
`r add_commas(sum(sapply(liu_gmap_3Mbp, function(a) length(seq(a[2], a[length(a)], by=0.01)))))`
positions.

Let's go with a 0.02 cM grid for now.

```{r form_the_grid}
grid_gmap <- lapply(liu_gmap_3Mbp, function(a) seq(a[2], a[length(a)], by=0.02))
```

And next get interpolated physical positions, rounding to basepairs.

```{r get_grid_bp}
grid_pmap <- interp_map(grid_gmap, liu_gmap_3Mbp, liu_pmap_3Mbp)
grid_pmap <- lapply(grid_pmap, round)
```

Now, we add names.

```{r grid_add_names}
for(i in seq_along(grid_pmap)) {
    names(grid_pmap[[i]]) <- names(grid_gmap[[i]]) <-
        paste(names(grid_pmap)[i], grid_pmap[[i]], sep="_")
}
```

And we put them back together as a data frame.

```{r grid_to_df}
grid <- cbind(map_list_to_df(grid_pmap, pos_column="bp"),
              cM=map_list_to_df(grid_gmap, pos_column="cM")[,"cM"])
rownames(grid) <- grid$marker
grid$pos <- grid$bp/1e6
grid <- grid[,c("marker", "chr", "pos", "cM", "bp")]
```

The resulting `grid` object has `r add_commas(nrow(grid))` rows.
The columns are slightly different than the `snps` object in the
`DO378_islet_v2.RData` file, in that the `snps` object has a column
called `bp` that is actually Mbp and a column `pos` that is basepairs.
The new `grid` object has the same content, but the third column with
Mbp positions is called `pos` and the final column with basepairs is
called `bp`.

### RNAseq gene positions on grid

We won't download the islet RNA-seq data, because it's not yet
officially public, but we grabbed the file `DO378_islet_v2.RData` from
the Dan Gatti's directory on the JAX ftp site.

```{r load_rnaseq, eval=FALSE}
file <- file.path(in_dir, "DO378_islet_v2.RData")
load(file)
```

```{r load_rnaseq_really, echo=FALSE}
file <- file.path(in_dir, "DO378_islet_v2.RData")
sub_file <- file.path(in_dir, "DO378_islet_v2_sub.RData")
if(file.exists(sub_file)) {
    load(sub_file)
} else {
    load(file)
    save(snps, annot.mrna, file=sub_file)
}
```

The object `snps` is the previous 64k grid with positions evenly spaced in cM.
(Again note that the spacing is slightly different between chromosomes.)

The object we want is `annot.mrna`. The `chr` column has 1-19, `X`, `Y`,
and `MT`. There are columns `start`, `end`, and `middle_point`, all in
basepairs, with `middle_point` being `round((start+end)/2)`.
The column `nearest_snp` has integers in the range 1-64k,
with missing values for the genes on the Y and MT; this is the index
of the nearest point on the 64k SNP grid. We want to replace those
values with the indices corresponding to our new grid.


Here's a function to do the job for a single position:

```{r find_nearest_grid_pos}
find_nearest <-
    function(chr, pos, map, index)
{
    if(!(chr %in% names(map))) return(NA)
    index[[chr]][ which.min( abs(pos - map[[chr]]) ) ]
}
```

Now we pull out the grid positions as a list and calculate a
corresponding list of grid indices.

```{r grid_pmap_and_index}
grid_pmap <- map_df_to_list(grid, pos_column="bp")
grid_index <- map_df_to_list(cbind(grid, index=1:nrow(grid)), pos_column="index")
```

Finally, we calculate the new grid indices and plug them into the object.

```{r calc_grid_indices}
index <- rep(NA, nrow(annot.mrna))
for(i in 1:nrow(annot.mrna)) {
    index[i] <- find_nearest(annot.mrna$chr[i], annot.mrna$middle_point[i],
                             grid_pmap, grid_index)
}
annot.mrna_newgrid <- annot.mrna
annot.mrna_newgrid$nearest_snp <- index
```


### One concern about the new grid

One potential issue to keep in mind: the Liu et al. map has some
long regions with very low recombination, and so this new 0.02 cM grid
has some pretty big physical gaps, particularly on chromosome 14.
Across the whole genome, there are
`r sum(sapply(grid_pmap, function(a) sum(diff(a) > 1e6)))`
gaps > 1 Mbp, and
`r sum(sapply(grid_pmap, function(a) sum(diff(a) > 2e6)))`
gaps > 2 Mbp, and then chromosome 14 has two gaps > 5 Mbp.

```{r calc_snps_pmap, echo=FALSE}
snps_pmap <- map_df_to_list(snps, pos_column="pos")
```

In contrast, KB Choi's original grid had
`r sum(sapply(snps_pmap, function(a) sum(diff(a) > 2e6)))`
gaps > 2 Mbp and just
`r sum(sapply(snps_pmap, function(a) sum(diff(a) > 1e6)))`
gaps > 1 Mbp.

One result of these physical gaps in the grid is that the many of the
genes in `annot.mrna` are rather far away from the closest grid position.

```{r calculate_stuff, echo=FALSE}
new_gridpos <- grid$bp[annot.mrna_newgrid$nearest_snp]
old_gridpos <- snps$pos[annot.mrna$nearest_snp]
```

With the original grid, of the `r add_commas(sum(!is.na(annot.mrna$middle_point)))`
genes on the autosomes or X chromosome,
`r sum(abs(annot.mrna$middle_point - old_gridpos) > 1e6, na.rm=TRUE)`
are more than 1 Mbp away from a grid position, and
just `r sum(abs(annot.mrna$middle_point - old_gridpos) > 0.5e6, na.rm=TRUE)`
are more than 500 kbp away from a grid position.

With the new grid,
`r sum(abs(annot.mrna$middle_point - new_gridpos) > 1e6, na.rm=TRUE)`
are more than 1 Mbp away from a grid position, and
`r sum(abs(annot.mrna$middle_point - new_gridpos) > 0.5e6, na.rm=TRUE)`
are more than 500 kbp away from a grid position. There are
`r sum(abs(annot.mrna$middle_point - new_gridpos) > 2e6, na.rm=TRUE)`
genes more than 2 Mbp from a grid position, all on chr 14.
The maximum distance is
`r round(max(abs(annot.mrna$middle_point - new_gridpos) , na.rm=TRUE)/1e6, 1)`
Mbp.


### Add more pseudomarkers

To deal with the physical gaps between grid points, I'm inclined to
add additional pseudomarkers so that no two points are more than 0.5
Mbp apart.

```{r add_more_pseudomarkers}
grid_pmap <- lapply(grid_pmap, function(a) a/1e6) # pmap in Mbp rather than bp

grid_pmap_plus <- insert_pseudomarkers(grid_pmap, step=0.5, stepwidth="max")
for(i in seq_along(grid_pmap_plus)) {
    names(grid_pmap_plus[[i]]) <- paste(names(grid_pmap)[i], round(grid_pmap_plus[[i]]*1e6), sep="_")
}
grid_gmap_plus <- interp_map(grid_pmap_plus, grid_pmap, grid_gmap)

grid_plus <- cbind(map_list_to_df(grid_pmap_plus, pos_column="pos"),
              cM=map_list_to_df(grid_gmap_plus, pos_column="cM")[,"cM"])
rownames(grid_plus) <- grid_plus$marker
grid_plus$bp <- round(grid_plus$pos*1e6)
grid_plus <- grid_plus[,c("marker", "chr", "pos", "cM", "bp")]
```

This increases the number of grid points from `r nrow(grid)` to `r nrow(grid_plus)`
(that is, an additional `r nrow(grid_plus)-nrow(grid)` points).


### Save stuff to files

Now let's save everything to files. I'm going to use `.rds` files rather
than `.RData` files, as they each will contain one object.

```{r save_maps}
saveRDS(liu_map, file.path(out_dir, "liu_map.rds"))
saveRDS(GM_snps, file.path(out_dir, "GM_snps_v2.rds"))
saveRDS(MM_snps, file.path(out_dir, "MM_snps_v2.rds"))
saveRDS(grid, file.path(out_dir, "grid_0.02cM.rds"))
saveRDS(grid_plus, file.path(out_dir, "grid_0.02cM_plus.rds"))
saveRDS(annot.mrna_newgrid, file.path(out_dir, "annot_mrna_0.02cM_grid.rds"))
```



### Session Info

```{r session_info, results="markup"}
devtools::session_info()
```
