Mapping between genome assemblies in R

From the archive: joining marker coordinates from two assemblies of the dog genome (canFam2 and CanFam3) without losing marker order.

NoteFrom the archive

First published on my old blog in March 2014. I have kept the idea and the original plyr solution, and added a version with today’s dplyr.

Genome arrays are designed against one particular version of a genome assembly. New assemblies are published, but the array stays the same, so you often need to map between the two. In my case, the task was to map between two assemblies of the canine genome, from canFam2 to CanFam3, for a SNP array. I had two .map files with genomic coordinates per marker, which look like this:

chrom name position
1 BICF2G630707759 14548
1 BICF2S2358127 68723
1 BICF2P1173580 80040

I wanted to join them so that the marker names come from the canFam2 file, while the chromosome and position come from the matching marker in the canFam3 file.

The pitfall of merge()

The first thing that came to mind was merge(). The problem is that a marker that has no match in canFam3 is moved to the end of the result, which breaks the marker order. What we want instead is NA in both the chromosome and the position when the marker cannot be found in canFam3.

The solution: a left join

The R package plyr has exactly that, in its join() function:

new.map <- join(map, canFam3, by = "marker", type = "left")

Today the same thing is usually written with dplyr. Here it is with a tiny made-up example, so you can see that the order is preserved and the missing marker gets NA:

library(dplyr)

canFam2 <- data.frame(
  marker = c("BICF2G630707759", "BICF2S2358127", "BICF2P1173580", "BICF2X0000001"),
  chr2 = 1, pos2 = c(14548, 68723, 80040, 90111)
)
canFam3 <- data.frame(
  marker = c("BICF2P1173580", "BICF2G630707759", "BICF2S2358127"),
  chr3 = 1, pos3 = c(81233, 15021, 69300)
)

left_join(canFam2, canFam3, by = "marker")

           marker chr2  pos2 chr3  pos3
1 BICF2G630707759    1 14548    1 15021
2   BICF2S2358127    1 68723    1 69300
3   BICF2P1173580    1 80040    1 81233
4   BICF2X0000001    1 90111   NA    NA

The unmatched marker BICF2X0000001 stays in its place and gets NA for the new coordinates, which is what we wanted. In base R, match() does the same job: canFam3$pos3[match(canFam2$marker, canFam3$marker)].

TipWhen marker names are not shared

This trick works because the same array markers appear in both assemblies. To convert coordinates that have no shared identifiers, use a chain-file based lift-over instead, for example liftOver() in the Bioconductor package rtracklayer.