1. Prepare data

maps <- readRDS("../2_LinMap/10_4_linkageMap.RData")
library(LinkageMapView)
library(janitor)
## 
## Attaching package: 'janitor'
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.2.0     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1
## Warning: package 'tidyr' was built under R version 4.0.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

Convert to LinkageMapView format

LMV_male <- maps %>%
  select(LG,male,SNPID) %>%
  drop_na() %>%
  rename(group = LG, position = male,
           locus = SNPID) %>%
  mutate(across(where(is.factor),
                  as.character)) %>%
  mutate(group = paste0(group, "M"))

LMV_female <- maps %>%
  select(LG,female,SNPID) %>%
  drop_na() %>%
  rename(group = LG, position = female,
           locus = SNPID) %>%
  mutate(across(where(is.factor),
                  as.character)) %>%
  mutate(group = paste0(group, "F"))

LMV_average <- maps %>%
  select(LG,sex.averaged,SNPID) %>%
  drop_na() %>%
  rename(group = LG, position = sex.averaged,
           locus = SNPID) %>%
  mutate(across(where(is.factor),
                  as.character)) %>%
  mutate(group = paste0(group, "A"))

HC <- rbind(LMV_male, LMV_female, LMV_average)

levels <- c("1F","1M","1A","2F","2M","2A","3F","3M","3A",
            "4F","4M","4A","5F","5M","5A","6F","6M","6A",
            "7F","7M","7A","8F","8M","8A","9F","9M","9A",
            "10F","10M","10A","11F","11M","11A",
            "12F","12M","12A","13F","13M","13A")

HC$group <- fct_relevel(HC$group,
                        levels = levels)
## Warning: Outer names are only allowed for unnamed scalar atomic inputs
HC <- HC[order(HC$group,HC$position),]

2. Plot density maps

outfile = file.path("../2_LinMap/12_1_maleMap_density.pdf")

lmv.linkage.plot(LMV_male,outfile,denmap=TRUE,
                 pdf.width = 15, pdf.height = 15,
                 cex.axis = 2,lgw = 0.30,cex.lgtitle = 2,
                 font.main = 2,cex.main = 3,
                 main="Male linkage maps")
## Required pdf.width = 11.67
## Required pdf.height = 4.56388888888891
## Using pdf.width = 15
## Using pdf.height = 15
outfile = file.path("../2_LinMap/12_1_femaleMap_density.pdf")

lmv.linkage.plot(LMV_female,outfile,denmap=TRUE,
                 pdf.width = 15, pdf.height = 15,
                 cex.axis = 2,lgw = 0.30,cex.lgtitle = 2,
                 font.main = 2,cex.main = 3,
                 main="Female linkage maps")
## Required pdf.width = 10.708
## Required pdf.height = 6.23402777777799
## Using pdf.width = 15
## Using pdf.height = 15
outfile = file.path("../2_LinMap/12_1_averageMap_density.pdf")

lmv.linkage.plot(LMV_average,outfile,denmap=TRUE,
                 pdf.width = 15, pdf.height = 15,
                 cex.axis = 2,lgw = 0.30,cex.lgtitle = 2,
                 font.main = 2,cex.main = 3,
                 main="Sex-averaged linkage maps")
## Required pdf.width = 10.9506666666667
## Required pdf.height = 13.5656250000006
## Using pdf.width = 15
## Using pdf.height = 15

Sex-averaged linkage map colour coded by density

3. Plot black and white maps

Now let’s plot the male, female and sex-averaged linkage maps showing only the position of markers as black horizontal lines.

outfile = file.path("../2_LinMap/12_2_maleMap.pdf")

lmv.linkage.plot(LMV_male,outfile, showonly = c("22896_41"),
                 ruler = TRUE,cex.axis = 2,
                 lgw = 0.30,
                 lgtitles = c("LG1","LG2","LG3","LG4",
                              "LG5","LG6","LG7","LG8",
                              "LG9","LG10","LG11","LG12",
                              "LG13"),
                 cex.lgtitle = 2,
                 pdf.width = 20, pdf.height = 18,
                 main="Male linkage maps",
                 font.main = 2,cex.main = 3)
## Required pdf.width = 14.8961666666667
## Required pdf.height = 1.38506666666667
## Using pdf.width = 20
## Using pdf.height = 18
outfile = file.path("../2_LinMap/12_2_femaleMap.pdf")

lmv.linkage.plot(LMV_female,outfile, showonly = c("22896_41"),
                 ruler = TRUE,cex.axis = 2,
                 lgw = 0.30,
                 lgtitles = c("LG1","LG2","LG3","LG4",
                              "LG5","LG6","LG7","LG8",
                              "LG9","LG10","LG11","LG12",
                              "LG13"),
                 cex.lgtitle = 2,
                 pdf.width = 20, pdf.height = 18,
                 main="Female linkage maps",
                 font.main = 2,cex.main = 3)
## Required pdf.width = 14.8961666666667
## Required pdf.height = 1.5526
## Using pdf.width = 20
## Using pdf.height = 18
outfile = file.path("../2_LinMap/12_2_averageMap.pdf")

lmv.linkage.plot(LMV_average,outfile, showonly = c("22896_41"),
                 ruler = TRUE,cex.axis = 2,
                 lgw = 0.30,
                 lgtitles = c("LG1","LG2","LG3","LG4",
                              "LG5","LG6","LG7","LG8",
                              "LG9","LG10","LG11","LG12",
                              "LG13"),
                 cex.lgtitle = 2,
                 pdf.width = 20, pdf.height = 18,
                 main="Sex-averaged linkage maps",
                 font.main = 2,cex.main = 3)
## Required pdf.width = 14.8961666666667
## Required pdf.height = 1.5526
## Using pdf.width = 20
## Using pdf.height = 18

Sex-averaged linkage map

4. Plot male and female maps next to each other for size comparison

First, let’s create a dataframe with male and female maps, and then order them.

temp <- rbind(LMV_male, LMV_female) %>% arrange(group)

temp1 <- temp %>% filter(group != "10M" & group != "11M" &
                           group != "12M" & group != "13M" &
                           group != "10F" & group != "11F" &
                           group != "12F" & group != "13F")
temp2 <- temp %>% filter(group == "10M" | group == "11M" |
                           group == "12M" | group == "13M" |
                           group == "10F" | group == "11F" |
                           group == "12F" | group == "13F")

temp <- rbind(temp1, temp2)

And then plot

outfile = file.path("../2_LinMap/12_3_maleVSfemaleMap.pdf")

lmv.linkage.plot(temp,outfile, showonly = c("22896_41"),
                 ruler = TRUE,cex.axis = 3,
                 lgw = 0.30,
                 lgtitles = c("LG1f","LG1m","LG2f","LG2m",
                              "LG3f","LG3m","LG4f","LG4m",
                              "LG5f","LG5m","LG6f","LG6m",
                              "LG7f","LG7m","LG8f","LG8m",
                              "LG9f","LG9m","LG10f","LG10m",
                              "LG11f","LG11m","LG12f","LG12m",
                              "LG13f","LG13m"),
                 cex.lgtitle = 3,
                 pdf.width = 50, pdf.height = 20,
                 main="Female and male linkage maps",
                 font.main = 2,cex.main = 4)
## Required pdf.width = 28.6923333333333
## Required pdf.height = 1.7526
## Using pdf.width = 50
## Using pdf.height = 20

Female vs male linkage map

And now let’s plot the density map for LG1 for both the female and male next to each other to show the heterochiasmy.

temp <- temp %>% filter(group == "1F" | group == "1M")
outfile = file.path("../2_LinMap/12_4_maleVSfemaleLG1_density.pdf")

lmv.linkage.plot(temp, outfile, denmap=TRUE,
                 pdf.width = 3, pdf.height = 7,
                 cex.axis = 1, lgw = 0.30, cex.lgtitle = 1,
                 font.main = 2, cex.main = 1,
                 lgtitles = c("Female","Male"),
                 main="Female and male LG1")
## Required pdf.width = 1.80366666666667
## Required pdf.height = 6.03402777777799
## Using pdf.width = 3
## Using pdf.height = 7

Female versus male LG1 ## 5. Plot position in male versus female maps

rm(list = ls())

Read in map data

maps <- readRDS("../2_LinMap/10_4_linkageMap.RData")

Relevel in order of LG name

maps <- maps %>%
  mutate(LG = fct_relevel(LG, "1", "2", "3", "4",
                          "5", "6", "7", "8", "9",
                          "10", "11", "12", "13"))

And plot

p1 <- ggplot(maps,
             aes(x=female, y=male,
                 group = LG, colour = LG)) +
  geom_point(size=0.5) + theme_classic() +
  theme(legend.position = "none") +
  xlab("Female map position (cM)") +
  ylab("Male map position (cM)") +
  scale_y_continuous(limits = c(0, 170),
                     breaks = seq(0, 150, by = 50)) +
  scale_x_continuous(limits = c(0,265),
                     breaks = seq(0, 250, by = 50))

p1
## Warning: Removed 7015 rows containing missing values (geom_point).

Save the resulting plot

ggsave(filename = "../2_LinMap/12_5_MaleVsFemale_All.jpg",
       plot = p1, width = 6, height = 4)
## Warning: Removed 7015 rows containing missing values (geom_point).

That figure is a bit confusing, showing many LGs at once. Let’s plot each LG individually.

maps.NA <- maps %>% drop_na()

p2 <- ggplot(maps.NA,
             aes(x=female, y=male,
                         group = LG, colour = LG)) +
  geom_point() + facet_wrap(~LG, ncol = 5) +
  xlab("Female map position (cM)") +
  ylab("Male map position (cM)") +
  theme_classic() +
  theme(legend.position = "none")

p2

Save the resulting figure

ggsave(filename = "../2_LinMap/12_6_MaleVsFemale.jpg",
       plot = p2, width = 8, height = 6)
rm(p1,p2)

6. Plot marker density

Now, let’s plot marker density of male versus female maps to show the high level of heterochiasmy.