library(introgress) # load first to avoid masking select from dplyr
library(tidyverse)
library(forcats)
library(adegenet)
library(pegas)

Import data

G <- read_csv("Mytilus_crosses_genotypes.csv")

## Parsed with column specification:
## cols(
##   .default = col_character()
## )

## See spec(...) for full column specifications.

G$pop <- as_factor(G$pop)

Determine which markers are heterozygote for all F1 parents.

F1 <- G %>% filter(pop == "F1") %>% select(-c(1:3))

are_hz <- function(x){
    # are all genotypes in x heterozygotes?
    alleles <- as.data.frame(strsplit(x, "/"))
    return(all(alleles[1,] != alleles[2,]))
}
result_hz <- apply(F1, 2, are_hz)
hz_F1 <- names(result_hz)[result_hz == T] %>% .[-which(is.na(.))]

change_hz <- function(x){
    geno <- levels(as.factor(x))
    H <- geno[unlist(lapply(strsplit(geno, "/"), function(y){y[1] == y[2]})) == F]
    return(replace(x, x == H, "H"))
}

(a) 43 Markers, All.

Ga <- G %>% select(ind, pop, one_of(hz_F1)) # step 1

M <- apply(Ga %>% select(-ind, -pop), 2, change_hz)
M_F2 <- M[Ga$pop == "F2",]
dim(M_F2)

## [1] 132  43

M_BC <- M[Ga$pop %in% c("BCG","BCF1"),]
dim(M_BC)

## [1] 144  43

p12_F2 <- rowMeans(M_F2 == "H", na.rm = T)
p12_BC <- rowMeans(M_BC == "H", na.rm = T)

wilcox.test(p12_F2, mu = 0.5, conf.int = T)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  p12_F2
## V = 6238.5, p-value = 1.509e-06
## alternative hypothesis: true location is not equal to 0.5
## 95 percent confidence interval:
##  0.5465070 0.6046239
## sample estimates:
## (pseudo)median 
##      0.5699878

wilcox.test(p12_BC, mu = 0.5, conf.int = T)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  p12_BC
## V = 7214, p-value = 1.344e-05
## alternative hypothesis: true location is not equal to 0.5
## 95 percent confidence interval:
##  0.5397184 0.5930381
## sample estimates:
## (pseudo)median 
##      0.5697592

rm(M, M_F2, M_BC, p12_F2, p12_BC)

(b) 43 markers, No missing data.

Gb <- G %>% select(ind, pop, one_of(hz_F1)) # step 1
Gb <- Gb[rowSums(is.na(Gb)) == 0,] # step 4

M <- apply(Gb %>% select(-ind, -pop), 2, change_hz)
M_F2 <- M[Gb$pop == "F2",]
dim(M_F2)

## [1] 88 43

M_BC <- M[Gb$pop %in% c("BCG","BCF1"),]
dim(M_BC)

## [1] 94 43

p12_F2 <- rowMeans(M_F2 == "H", na.rm = T)
p12_BC <- rowMeans(M_BC == "H", na.rm = T)

wilcox.test(p12_F2, mu = 0.5, conf.int = T)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  p12_F2
## V = 2778, p-value = 0.0006411
## alternative hypothesis: true location is not equal to 0.5
## 95 percent confidence interval:
##  0.5233058 0.5814543
## sample estimates:
## (pseudo)median 
##      0.5581174

wilcox.test(p12_BC, mu  = 0.5, conf.int = T)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  p12_BC
## V = 2814.5, p-value = 0.0282
## alternative hypothesis: true location is not equal to 0.5
## 95 percent confidence interval:
##  0.5000284 0.5697639
## sample estimates:
## (pseudo)median 
##      0.5349032

rm(M, M_F2, M_BC, p12_F2, p12_BC)

(c) 33 markers, No missing data.

Gc <- G %>% select(ind, pop, one_of(hz_F1)) # step 1

# step 2
gt <- loci2genind(as.loci(Gc %>% select(-ind, -pop)))
gt@pop <- Gc$pop
indNames(gt) <- Gc$ind
freq <- makefreq(gt, quiet = T)
f_edu <- colMeans(freq[gt@pop == "Edu",], na.rm = T)
f_gallo <- colMeans(freq[gt@pop == "Gallo",], na.rm = T)
f_diff <- abs(f_edu - f_gallo)
high_freq <- sub("\\..*$", "", names(f_diff[f_diff > 0.9])) %>% .[duplicated(.)]
Gc <- Gc %>% select(ind, pop, one_of(high_freq))

Gc <- Gc[rowSums(is.na(Gc)) == 0,] # step 4

M <- apply(Gc %>% select(-ind, -pop), 2, change_hz)
M_F2 <- M[Gc$pop == "F2",]
dim(M_F2)

## [1] 91 33

M_BC <- M[Gc$pop %in% c("BCG","BCF1"),]
dim(M_BC)

## [1] 105  33

p12_F2 <- rowMeans(M_F2 == "H", na.rm = T)
p12_BC <- rowMeans(M_BC == "H", na.rm = T)

wilcox.test(p12_F2, mu = 0.5, conf.int = T)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  p12_F2
## V = 2834.5, p-value = 0.003304
## alternative hypothesis: true location is not equal to 0.5
## 95 percent confidence interval:
##  0.5151528 0.5757294
## sample estimates:
## (pseudo)median 
##      0.5454429

wilcox.test(p12_BC, mu = 0.5, conf.int = T)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  p12_BC
## V = 3378, p-value = 0.05689
## alternative hypothesis: true location is not equal to 0.5
## 95 percent confidence interval:
##  0.4999780 0.5606414
## sample estimates:
## (pseudo)median 
##       0.530371

rm(M, M_F2, M_BC, p12_F2, p12_BC)

(d) 33 markers, Subsampled.

# heterozygote in F1
Gd <- G %>% select(ind, pop, one_of(hz_F1))

# high freq
gt <- loci2genind(as.loci(Gd %>% select(-ind, -pop)))
gt@pop <- Gd$pop
indNames(gt) <- Gd$ind
freq <- makefreq(gt, quiet = T)
f_edu <- colMeans(freq[gt@pop == "Edu",], na.rm = T)
f_gallo <- colMeans(freq[gt@pop == "Gallo",], na.rm = T)
f_diff <- abs(f_edu - f_gallo)
high_freq <- sub("\\..*$", "", names(f_diff[f_diff > 0.9])) %>% .[duplicated(.)]
Gd <- Gd %>% select(ind, pop, one_of(high_freq))

# missing data
Gd <- Gd[rowSums(is.na(Gd)) == 0,]

# get edulis major allele
get.edulis.major <- function(z)
{
    a <- apply(z,MARGIN=2,FUN=function(f) {t <- table(f); return(names(t)[which.max(t)])})
    return(a)
}
em <- get.edulis.major(Gd %>% filter(pop == "Edu") %>% select(-c(1,2)))

M <- apply(Gd %>% select(-ind, -pop), 2, change_hz)
M_BC <- M[Gd$pop %in% c("BCG","BCF1"),]
# remove F2-like BCs
to.keep <- rep(0,dim(M_BC)[1])
for(i in 1:length(em)) {
    #print(i)
    if(any(M_BC[,i]==em[i],na.rm=T)) {
        to.keep[which(M_BC[,i]==em[i])] <- to.keep[which(M_BC[,i]==em[i])] + 1
    }
}
to.keep <- which(to.keep <= 0)
M_BC <- M_BC[to.keep,]
dim(M_BC)

## [1] 56 33

# downsample F2
M_F2 <- M[Gd$pop == "F2",]
M_F2 <- M_F2[1:nrow(M_BC),]
dim(M_F2)

## [1] 56 33

p12_F2 <- rowMeans(M_F2 == "H", na.rm = T)
p12_BC <- rowMeans(M_BC == "H", na.rm = T)

wilcox.test(p12_F2, mu = 0.5, conf.int = T)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  p12_F2
## V = 1177, p-value = 0.001981
## alternative hypothesis: true location is not equal to 0.5
## 95 percent confidence interval:
##  0.5151995 0.5909200
## sample estimates:
## (pseudo)median 
##      0.5605687

wilcox.test(p12_BC, mu = 0.5, conf.int = T)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  p12_BC
## V = 866, p-value = 0.5815
## alternative hypothesis: true location is not equal to 0.5
## 95 percent confidence interval:
##  0.4545617 0.5606001
## sample estimates:
## (pseudo)median 
##      0.5151206

rm(M, M_F2, M_BC, p12_F2, p12_BC, BC)

## Warning in rm(M, M_F2, M_BC, p12_F2, p12_BC, BC): objet 'BC' introuvable

Figure S6

G_ordered <- G
G_ordered$pop <- fct_relevel(G_ordered$pop, "Edu","BCF1","BCG","F2","F1","Gallo")

admix <- t(G_ordered %>% select(-c(1:3)))
P1 <- t(G_ordered %>% filter(pop == "Edu") %>% select(-c(1:3)))
P2 <- t(G_ordered %>% filter(pop == "Gallo") %>% select(-c(1:3)))
loci <- cbind(locus = colnames(G_ordered)[-c(1:3)], type = rep("c", nrow(admix)))

dat <- prepare.data(admix.gen = admix, loci.data = loci, parental1 = P1,
                    parental2 = P2, pop.id = FALSE, ind.id = FALSE)

indice.h <- est.h(introgress.data=dat, loci.data=loci[,c('locus','type')])

marker_order <- c(loci[loci[,1] %in% high_freq, 1],
                  loci[loci[,1] %in% hz_F1 & !(loci[,1] %in% high_freq), 1])
marker_order <- c(marker_order, loci[!loci[,1] %in% marker_order, 1])
marker_colors <- c(rep("red", 33), rep("orange", 10), rep("black", 55))
pop_labels <- c(expression(italic("M. edulis")),
                expression(BC[12]),
                expression(BC[21]), 
                "F2", 
                "F1",
                expression(italic("M. gallo.")))

source("~/These/Scripts/mk.imageMOD.R")
mk.imageMOD(introgress.data = dat, loci.data = loci,
            hi.index=indice.h, ind.touse=NULL, loci.touse=NULL,
            ylab.image="Populations", main.image="",
            xlab.h="population 2 ancestry",
            pdf=TRUE, out.file = "output/FigS3.pdf",
            pop = G_ordered$pop, pop.labels = pop_labels,
            marker.order = marker_order, marker.colors = marker_colors)

devtools::session_info()

## Session info -------------------------------------------------------------

##  setting  value                       
##  version  R version 3.5.0 (2018-04-23)
##  system   x86_64, darwin15.6.0        
##  ui       X11                         
##  language (EN)                        
##  collate  fr_FR.UTF-8                 
##  tz       Europe/London               
##  date     2018-06-28

## Packages -----------------------------------------------------------------

##  package      * version date       source        
##  ade4         * 1.7-11  2018-04-05 CRAN (R 3.5.0)
##  adegenet     * 2.1.1   2018-02-02 CRAN (R 3.5.0)
##  ape          * 5.1     2018-04-04 CRAN (R 3.5.0)
##  assertthat     0.2.0   2017-04-11 CRAN (R 3.5.0)
##  backports      1.1.2   2017-12-13 CRAN (R 3.5.0)
##  base         * 3.5.0   2018-04-24 local         
##  bindr          0.1.1   2018-03-13 CRAN (R 3.5.0)
##  bindrcpp     * 0.2.2   2018-03-29 CRAN (R 3.5.0)
##  boot           1.3-20  2017-08-06 CRAN (R 3.5.0)
##  broom          0.4.4   2018-03-29 CRAN (R 3.5.0)
##  cellranger     1.1.0   2016-07-27 CRAN (R 3.5.0)
##  cli            1.0.0   2017-11-05 CRAN (R 3.5.0)
##  cluster        2.0.7-1 2018-04-13 CRAN (R 3.5.0)
##  coda           0.19-1  2016-12-08 CRAN (R 3.5.0)
##  colorspace     1.3-2   2016-12-14 CRAN (R 3.5.0)
##  combinat     * 0.0-8   2012-10-29 CRAN (R 3.5.0)
##  compiler       3.5.0   2018-04-24 local         
##  crayon         1.3.4   2017-09-16 CRAN (R 3.5.0)
##  datasets     * 3.5.0   2018-04-24 local         
##  deldir         0.1-15  2018-04-01 CRAN (R 3.5.0)
##  devtools       1.13.5  2018-02-18 CRAN (R 3.5.0)
##  digest         0.6.15  2018-01-28 CRAN (R 3.5.0)
##  dplyr        * 0.7.4   2017-09-28 CRAN (R 3.5.0)
##  evaluate       0.10.1  2017-06-24 CRAN (R 3.5.0)
##  expm           0.999-2 2017-03-29 CRAN (R 3.5.0)
##  forcats      * 0.3.0   2018-02-19 CRAN (R 3.5.0)
##  foreign        0.8-70  2017-11-28 CRAN (R 3.5.0)
##  gdata        * 2.18.0  2017-06-06 CRAN (R 3.5.0)
##  genetics     * 1.3.8.1 2013-09-03 CRAN (R 3.5.0)
##  ggplot2      * 2.2.1   2016-12-30 CRAN (R 3.5.0)
##  glue           1.2.0   2017-10-29 CRAN (R 3.5.0)
##  gmodels        2.16.2  2015-07-22 CRAN (R 3.5.0)
##  graphics     * 3.5.0   2018-04-24 local         
##  grDevices    * 3.5.0   2018-04-24 local         
##  grid           3.5.0   2018-04-24 local         
##  gtable         0.2.0   2016-02-26 CRAN (R 3.5.0)
##  gtools       * 3.5.0   2015-05-29 CRAN (R 3.5.0)
##  haven          1.1.1   2018-01-18 CRAN (R 3.5.0)
##  hms            0.4.2   2018-03-10 CRAN (R 3.5.0)
##  htmltools      0.3.6   2017-04-28 CRAN (R 3.5.0)
##  httpuv         1.4.2   2018-05-03 CRAN (R 3.5.0)
##  httr           1.3.1   2017-08-20 CRAN (R 3.5.0)
##  igraph         1.2.1   2018-03-10 CRAN (R 3.5.0)
##  introgress   * 1.2.3   2012-10-29 CRAN (R 3.5.0)
##  jsonlite       1.5     2017-06-01 CRAN (R 3.5.0)
##  knitr          1.20    2018-02-20 CRAN (R 3.5.0)
##  later          0.7.2   2018-05-01 CRAN (R 3.5.0)
##  lattice        0.20-35 2017-03-25 CRAN (R 3.5.0)
##  lazyeval       0.2.1   2017-10-29 CRAN (R 3.5.0)
##  LearnBayes     2.15.1  2018-03-18 CRAN (R 3.5.0)
##  lubridate      1.7.4   2018-04-11 CRAN (R 3.5.0)
##  magrittr       1.5     2014-11-22 CRAN (R 3.5.0)
##  MASS         * 7.3-49  2018-02-23 CRAN (R 3.5.0)
##  Matrix         1.2-14  2018-04-13 CRAN (R 3.5.0)
##  memoise        1.1.0   2017-04-21 CRAN (R 3.5.0)
##  methods      * 3.5.0   2018-04-24 local         
##  mgcv           1.8-23  2018-01-21 CRAN (R 3.5.0)
##  mime           0.5     2016-07-07 CRAN (R 3.5.0)
##  mnormt         1.5-5   2016-10-15 CRAN (R 3.5.0)
##  modelr         0.1.1   2017-07-24 CRAN (R 3.5.0)
##  munsell        0.4.3   2016-02-13 CRAN (R 3.5.0)
##  mvtnorm      * 1.0-8   2018-05-31 CRAN (R 3.5.0)
##  nlme           3.1-137 2018-04-07 CRAN (R 3.5.0)
##  nnet         * 7.3-12  2016-02-02 CRAN (R 3.5.0)
##  parallel       3.5.0   2018-04-24 local         
##  pegas        * 0.10    2017-05-03 CRAN (R 3.5.0)
##  permute        0.9-4   2016-09-09 CRAN (R 3.5.0)
##  pillar         1.2.2   2018-04-26 CRAN (R 3.5.0)
##  pkgconfig      2.0.1   2017-03-21 CRAN (R 3.5.0)
##  plyr           1.8.4   2016-06-08 CRAN (R 3.5.0)
##  promises       1.0.1   2018-04-13 CRAN (R 3.5.0)
##  psych          1.8.4   2018-05-06 CRAN (R 3.5.0)
##  purrr        * 0.2.4   2017-10-18 CRAN (R 3.5.0)
##  R6             2.2.2   2017-06-17 CRAN (R 3.5.0)
##  RColorBrewer * 1.1-2   2014-12-07 CRAN (R 3.5.0)
##  Rcpp           0.12.16 2018-03-13 CRAN (R 3.5.0)
##  readr        * 1.1.1   2017-05-16 CRAN (R 3.5.0)
##  readxl         1.1.0   2018-04-20 CRAN (R 3.5.0)
##  reshape2       1.4.3   2017-12-11 CRAN (R 3.5.0)
##  rlang          0.2.0   2018-02-20 CRAN (R 3.5.0)
##  rmarkdown      1.9     2018-03-01 CRAN (R 3.5.0)
##  rprojroot      1.3-2   2018-01-03 CRAN (R 3.5.0)
##  rstudioapi     0.7     2017-09-07 CRAN (R 3.5.0)
##  rvest          0.3.2   2016-06-17 CRAN (R 3.5.0)
##  scales         0.5.0   2017-08-24 CRAN (R 3.5.0)
##  seqinr         3.4-5   2017-08-01 CRAN (R 3.5.0)
##  shiny          1.0.5   2017-08-23 CRAN (R 3.5.0)
##  sp             1.2-7   2018-01-19 CRAN (R 3.5.0)
##  spData         0.2.8.3 2018-03-25 CRAN (R 3.5.0)
##  spdep          0.7-7   2018-04-03 CRAN (R 3.5.0)
##  splines        3.5.0   2018-04-24 local         
##  stats        * 3.5.0   2018-04-24 local         
##  stringi        1.2.2   2018-05-02 CRAN (R 3.5.0)
##  stringr      * 1.3.0   2018-02-19 CRAN (R 3.5.0)
##  tibble       * 1.4.2   2018-01-22 CRAN (R 3.5.0)
##  tidyr        * 0.8.0   2018-01-29 CRAN (R 3.5.0)
##  tidyverse    * 1.2.1   2017-11-14 CRAN (R 3.5.0)
##  tools          3.5.0   2018-04-24 local         
##  utils        * 3.5.0   2018-04-24 local         
##  vegan          2.5-1   2018-04-14 CRAN (R 3.5.0)
##  withr          2.1.2   2018-03-15 CRAN (R 3.5.0)
##  xml2           1.2.0   2018-01-24 CRAN (R 3.5.0)
##  xtable         1.8-2   2016-02-05 CRAN (R 3.5.0)
##  yaml           2.1.19  2018-05-01 CRAN (R 3.5.0)

Crossing experiment analysis

Alexis Simon

28/06/2018

Import data

(a) 43 Markers, All.

(b) 43 markers, No missing data.

(c) 33 markers, No missing data.

(d) 33 markers, Subsampled.

Figure S6