R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

———————- 1. Transform .vcf to genlight

Import clean .vcf:

library(vcfR)

package <U+393C><U+3E31>vcfR<U+393C><U+3E32> was built under R version 3.5.3
   *****       ***   vcfR   ***       *****
   This is vcfR 1.8.0 
     browseVignettes('vcfR') # Documentation
     citation('vcfR') # Citation
   *****       *****      *****       *****

vcf_clean <- read.vcfR(file = "C:/Users/drob0006/Google Drive/Projects/1.MateChoice_BreedingManagement/Analyses2020/0_2020.02.11_SNPcalling_rawDArTs/3.sex-linked/p1r9_ld_50kb_1kb_removed_missing_0.61_no_sexlinked_HW.recode.vcf")

Scanning file to determine attributes.
File attributes:
  meta lines: 395
  header_line: 396
  variant count: 33168
  column count: 542

Meta line 395 read in.
All meta lines processed.
gt matrix initialized.
Character matrix gt created.
  Character matrix gt rows: 33168
  Character matrix gt cols: 542
  skip: 0
  nrows: 33168
  row_num: 0

Processed variant 1000
Processed variant 2000
Processed variant 3000
Processed variant 4000
Processed variant 5000
Processed variant 6000
Processed variant 7000
Processed variant 8000
Processed variant 9000
Processed variant 10000
Processed variant 11000
Processed variant 12000
Processed variant 13000
Processed variant 14000
Processed variant 15000
Processed variant 16000
Processed variant 17000
Processed variant 18000
Processed variant 19000
Processed variant 20000
Processed variant 21000
Processed variant 22000
Processed variant 23000
Processed variant 24000
Processed variant 25000
Processed variant 26000
Processed variant 27000
Processed variant 28000
Processed variant 29000
Processed variant 30000
Processed variant 31000
Processed variant 32000
Processed variant 33000
Processed variant: 33168
All variants processed

Transform .vcf to genlight:

gl <- vcfR2genlight(vcf_clean)

Loading required namespace: adegenet

Let’s inspect the gl object:

gl

 /// GENLIGHT OBJECT /////////

 // 533 genotypes,  33,168 binary SNPs, size: 9.2 Mb
 436558 (2.47 %) missing data

 // Basic content
   @gen: list of 533 SNPbin

 // Optional content
   @ind.names:  533 individual labels
   @loc.names:  33168 locus labels
   @chromosome: factor storing chromosomes of the SNPs
   @position: integer storing positions of the SNPs
   @other: a list containing: elements without names 

Change names (genotyping ID) for BandID

# Check that order is correct
sum(gl@ind.names == new$gl.ind.names)

[1] 533

———————– 2. Subset census per season

This census was modified manually to have each individual per season only once, and if they were present in more than one neighbourhood, I fused the neighbourhood label

census <- read.csv("C:/Users/drob0006/Google Drive/Projects/1.MateChoice_BreedingManagement/Analyses2020/11_2020.04.21_PCA/Census_all_seasons_neighbourhood.csv",
                  colClasses = c("factor",
                                 "factor",
                                 "NULL",
                                 "character",
                                 "factor",
                                 rep("NULL", 3)))
                   
censi <- list()
for (i in 1:length(levels(census$Season)))
  censi[[i]] <- droplevels(census[census$Season == levels(census$Season)[i], ])

———————- 3. Make recode table per season

seasons <- c("2011.2012",
             "2012.2013",
             "2013.2014",
             "2014.2015",
             "2015.2016",
             "2016.2017",
             "2017.2018")
for (i in 1:length(recodes))
  write.csv(recodes[[i]], file = paste("C:/Users/drob0006/Google Drive/Projects/1.MateChoice_BreedingManagement/Analyses2020/11_2020.04.21_PCA/recode_",
                                       seasons[i],
                                       ".csv",
                                       sep = ""), 
            row.names = FALSE,
            col.names = FALSE)

attempt to set 'col.names' ignoredattempt to set 'col.names' ignoredattempt to set 'col.names' ignoredattempt to set 'col.names' ignoredattempt to set 'col.names' ignoredattempt to set 'col.names' ignoredattempt to set 'col.names' ignored

———————- 4. Subset gl per season

gl_season <- list()
for (i in 1:length(recodes))
  gl_season[[i]] <- gl.recode.ind(gl, 
                                  ind.recode = paste("C:/Users/drob0006/Google Drive/Projects/1.MateChoice_BreedingManagement/Analyses2020/11_2020.04.21_PCA/recode_",
                                       seasons[i],
                                       ".csv",
                                       sep = ""))

Processing genlight object
  Relabelling individuals (=specimens) as per  C:/Users/drob0006/Google Drive/Projects/1.MateChoice_BreedingManagement/Analyses2020/11_2020.04.21_PCA/recode_2011.2012.csv 
Deleting individuals or samples flagged for deletion
Starting gl.filter.monomorphs: Deleting monomorphic loci
Completed gl.filter.monomorphs

Starting utils.recalc.avgpic: Recalculating OneRatioRef, OneRatioSnp, PICRef, PICSnp and AvgPIC
Completed utils.recalc.avgpic

Starting utils.recalc.callrate: Recalculating CallRate
Completed utils.recalc.callrate

Starting utils.recalc.freqhets: Recalculating frequency of heterozygotes
Completed utils.recalc.freqhets

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, reference allele
Completed utils.recalc.freqhomref

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, alternate allele
Completed utils.recalc.freqhomsnp

Processing genlight object
  Relabelling individuals (=specimens) as per  C:/Users/drob0006/Google Drive/Projects/1.MateChoice_BreedingManagement/Analyses2020/11_2020.04.21_PCA/recode_2012.2013.csv 
Deleting individuals or samples flagged for deletion
Starting gl.filter.monomorphs: Deleting monomorphic loci
Completed gl.filter.monomorphs

Starting utils.recalc.avgpic: Recalculating OneRatioRef, OneRatioSnp, PICRef, PICSnp and AvgPIC
Completed utils.recalc.avgpic

Starting utils.recalc.callrate: Recalculating CallRate
Completed utils.recalc.callrate

Starting utils.recalc.freqhets: Recalculating frequency of heterozygotes
Completed utils.recalc.freqhets

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, reference allele
Completed utils.recalc.freqhomref

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, alternate allele
Completed utils.recalc.freqhomsnp

Processing genlight object
  Relabelling individuals (=specimens) as per  C:/Users/drob0006/Google Drive/Projects/1.MateChoice_BreedingManagement/Analyses2020/11_2020.04.21_PCA/recode_2013.2014.csv 
Deleting individuals or samples flagged for deletion
Starting gl.filter.monomorphs: Deleting monomorphic loci
Completed gl.filter.monomorphs

Starting utils.recalc.avgpic: Recalculating OneRatioRef, OneRatioSnp, PICRef, PICSnp and AvgPIC
Completed utils.recalc.avgpic

Starting utils.recalc.callrate: Recalculating CallRate
Completed utils.recalc.callrate

Starting utils.recalc.freqhets: Recalculating frequency of heterozygotes
Completed utils.recalc.freqhets

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, reference allele
Completed utils.recalc.freqhomref

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, alternate allele
Completed utils.recalc.freqhomsnp

Processing genlight object
  Relabelling individuals (=specimens) as per  C:/Users/drob0006/Google Drive/Projects/1.MateChoice_BreedingManagement/Analyses2020/11_2020.04.21_PCA/recode_2014.2015.csv 
Deleting individuals or samples flagged for deletion
Starting gl.filter.monomorphs: Deleting monomorphic loci
Completed gl.filter.monomorphs

Starting utils.recalc.avgpic: Recalculating OneRatioRef, OneRatioSnp, PICRef, PICSnp and AvgPIC
Completed utils.recalc.avgpic

Starting utils.recalc.callrate: Recalculating CallRate
Completed utils.recalc.callrate

Starting utils.recalc.freqhets: Recalculating frequency of heterozygotes
Completed utils.recalc.freqhets

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, reference allele
Completed utils.recalc.freqhomref

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, alternate allele
Completed utils.recalc.freqhomsnp

Processing genlight object
  Relabelling individuals (=specimens) as per  C:/Users/drob0006/Google Drive/Projects/1.MateChoice_BreedingManagement/Analyses2020/11_2020.04.21_PCA/recode_2015.2016.csv 
Deleting individuals or samples flagged for deletion
Starting gl.filter.monomorphs: Deleting monomorphic loci
Completed gl.filter.monomorphs

Starting utils.recalc.avgpic: Recalculating OneRatioRef, OneRatioSnp, PICRef, PICSnp and AvgPIC
Completed utils.recalc.avgpic

Starting utils.recalc.callrate: Recalculating CallRate
Completed utils.recalc.callrate

Starting utils.recalc.freqhets: Recalculating frequency of heterozygotes
Completed utils.recalc.freqhets

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, reference allele
Completed utils.recalc.freqhomref

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, alternate allele
Completed utils.recalc.freqhomsnp

Processing genlight object
  Relabelling individuals (=specimens) as per  C:/Users/drob0006/Google Drive/Projects/1.MateChoice_BreedingManagement/Analyses2020/11_2020.04.21_PCA/recode_2016.2017.csv 
Deleting individuals or samples flagged for deletion
Starting gl.filter.monomorphs: Deleting monomorphic loci
Completed gl.filter.monomorphs

Starting utils.recalc.avgpic: Recalculating OneRatioRef, OneRatioSnp, PICRef, PICSnp and AvgPIC
Completed utils.recalc.avgpic

Starting utils.recalc.callrate: Recalculating CallRate
Completed utils.recalc.callrate

Starting utils.recalc.freqhets: Recalculating frequency of heterozygotes
Completed utils.recalc.freqhets

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, reference allele
Completed utils.recalc.freqhomref

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, alternate allele
Completed utils.recalc.freqhomsnp

Processing genlight object
  Relabelling individuals (=specimens) as per  C:/Users/drob0006/Google Drive/Projects/1.MateChoice_BreedingManagement/Analyses2020/11_2020.04.21_PCA/recode_2017.2018.csv 
Deleting individuals or samples flagged for deletion
Starting gl.filter.monomorphs: Deleting monomorphic loci
Completed gl.filter.monomorphs

Starting utils.recalc.avgpic: Recalculating OneRatioRef, OneRatioSnp, PICRef, PICSnp and AvgPIC
Completed utils.recalc.avgpic

Starting utils.recalc.callrate: Recalculating CallRate
Completed utils.recalc.callrate

Starting utils.recalc.freqhets: Recalculating frequency of heterozygotes
Completed utils.recalc.freqhets

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, reference allele
Completed utils.recalc.freqhomref

Starting utils.recalc.freqhomref: Recalculating frequency of homozygotes, alternate allele
Completed utils.recalc.freqhomsnp

————- 5. Add individual metadata file with id, pop, sex

ind.metadata_season <- list()
for (i in 1:length(recodes)) {
  ind.metadata_season[[i]] <- merge(recodes[[i]],
                                    censi[[i]][, -1],
                                    by.x = "V1",
                                    by.y = "Band",
                                    sort = FALSE)
  ind.metadata_season[[i]] <- ind.metadata_season[[i]][, -2]
  names(ind.metadata_season[[i]]) <- c("id", "pop", "sex")
  # Check if they have same number of individuals
  print(nrow(ind.metadata_season[[i]]))
  print(length(gl_season[[i]]@ind.names))
  # Check if the order matches
  print(sum(ind.metadata_season[[i]]$id == gl_season[[i]]@ind.names))
}

[1] 36
[1] 36
[1] 36
[1] 27
[1] 27
[1] 27
[1] 49
[1] 49
[1] 49
[1] 64
[1] 64
[1] 64
[1] 79
[1] 79
[1] 79
[1] 94
[1] 94
[1] 94
[1] 118
[1] 118
[1] 118

Now we assign the ind.metadata as the @other$ind.metrics section of the genlight objects:

for (i in 1:length(gl_season)) {
  # Assign as individual metafile
  gl_season[[i]]@other$ind.metrics <- ind.metadata_season[[i]]
  # And as pop
  gl_season[[i]]@pop <- ind.metadata_season[[i]]$pop
  # Check if everything is correct
  print(sum(gl_season[[i]]@ind.names == gl_season[[i]]@other$ind.metrics$id))
}

[1] 36
[1] 27
[1] 49
[1] 64
[1] 79
[1] 94
[1] 118

pop(gl_season[[1]])

 [1] Blue       Blue       Blue       Blue       Blue       Blue       Blue       Green      Blue      
[10] Blue       Green      Green      Green      Blue/Green Blue       Blue/Green Red        Green     
[19] Green      Blue/Green Green      Blue       Blue       Blue       Blue/Green Green      Green     
[28] Blue       Green      Blue       Blue       Blue       Blue       Blue       Blue       Green     
Levels: Blue Blue/Green Green Red

gl_season[[1]]

 /// GENLIGHT OBJECT /////////

 // 36 genotypes,  23,570 binary SNPs, size: 5.4 Mb
 15287 (1.8 %) missing data

 // Basic content
   @gen: list of 36 SNPbin

 // Optional content
   @ind.names:  36 individual labels
   @loc.names:  23570 locus labels
   @chromosome: factor storing chromosomes of the SNPs
   @position: integer storing positions of the SNPs
   @pop: population of each individual (group size range: 1-20)
   @other: a list containing: loc.metrics  ind.metrics 

———————————– 6. PCA

Modifying dartR function for generating PCA biplot

#install.packages("directlabels")
#library(directlabels)
biplot.pop <- function(glPca, data, scale = FALSE, ellipse = FALSE, p = 0.95, 
                       labels = "pop", hadjust = 1.5, vadjust = 1, xaxis = 1, 
                       yaxis = 2) {
  # Sanity checks
  if (class(glPca) != "glPca" | class(data) != "genlight") {
    cat("Fatal Error: glPca and genlight objects required for glPca and data parameters respectively!\n")
    stop()
  }
  if (labels != "pop") {
    cat("Fatal Error: This function is only for plotting populations!\n")
    stop()
  }
  library(ggplot2)
  #library(directlabels)
  # Parameters
  m <- cbind(glPca$scores[, xaxis], glPca$scores[, yaxis])
  df <- data.frame(m)
  s <- sum(glPca$eig)
  e <- round(glPca$eig * 100/s, 1)
  xlab <- paste("PCoA  axis  ", xaxis, "  (", e[xaxis], "%)", sep = "")
  ylab <- paste("PCoA  axis  ", yaxis, "  (", e[yaxis], "%)", sep = "")
  
  # Plotting
  if (labels == "pop") {
    cat("Plotting populations\n")
    ind <- indNames(data)
    Neighbourhood <- factor(pop(data))
    df <- cbind(df, ind, factor(pop(data)))
    colnames(df) <- c("PCoAx", "PCoAy", "ind", "pop")
    
    p <- ggplot(df, 
                aes(x = df$PCoAx, y = df$PCoAy, group = Neighbourhood, colour = Neighbourhood)) + 
      geom_point(size = 3, aes(colour = Neighbourhood)) +
      scale_color_manual(values = c("Blue"           = "deepskyblue",
                                    "Blue/Green"     = "turquoise",
                                    "Blue/Green/Red" = "blueviolet",
                                    "Blue/Red"       = "magenta",
                                    "Blue/Woori1"    = "mediumblue",
                                    "Blue/Woori3"    = "lawngreen",
                                    "Green"          = "forestgreen",
                                    "Green/Red"      = "goldenrod",
                                    "JESWA"          = "gray",
                                    "Red"            = "red",
                                    "Woori1"         = "black",
                                    "Woori1/Woori2"  = "sienna",
                                    "Woori2"         = "darkorange",
                                    "Woori3"         = "yellow",
                                    "WY&CCC"         = "pink")) + 
      #geom_dl(aes(label = pop), method = "smart.grid") + 
      theme(axis.title   = element_text(face = "bold", size = "16", 
                                          color = "black"), 
            axis.text.x  = element_text(face = "bold", angle = 0, vjust = 0.5, 
                                          size = 59), 
            axis.text.y  = element_text(face = "bold", angle = 0, vjust = 0.5, 
                                          size = 14), 
            legend.title = element_text(colour = "black", size = 18, face = "bold"), 
            legend.text  = element_text(colour = "black", size = 18, face = "bold")) + 
      labs(x = xlab, y = ylab) + 
      geom_hline(yintercept = 0) + 
      geom_vline(xintercept = 0) + 
      theme(legend.position = "none") + 
      theme_bw()
    if (scale == TRUE) {
      p <- p + coord_fixed(ratio = e[yaxis]/e[xaxis])
    }
    if (ellipse == TRUE) {
      p <- p + stat_ellipse(aes(colour = "black"), type = "norm", 
        level = 0.95)
    }
  }
  p
  return(p)
}

2011.2012

#pc1 <- gl.pcoa(gl_season[[1]], nfactors=5)
barplot(pc1$eig/sum(pc1$eig)*100, 
        xlab = "PC",
        ylab = "Percentage of represented variation")

      
biplot.pop(pc1, gl_season[[1]], labels = "pop", xaxis = 1, yaxis = 2)

Plotting populations

2012.2013

#pc2 <- gl.pcoa(gl_season[[2]], nfactors=5)
barplot(pc2$eig/sum(pc2$eig)*100, 
        xlab = "PC",
        ylab = "Percentage of represented variation")

      
biplot.pop(pc2, gl_season[[2]], labels = "pop", xaxis = 1, yaxis = 2)

Plotting populations

2013.2014

#pc3 <- gl.pcoa(gl_season[[3]], nfactors = 5)
barplot(pc3$eig/sum(pc3$eig)*100, 
        xlab = "PC",
        ylab = "Percentage of represented variation")

      
biplot.pop(pc3, gl_season[[3]], labels = "pop", xaxis = 1, yaxis = 2)

Plotting populations

2014.2015

#pc4 <- gl.pcoa(gl_season[[4]], nfactors = 5)
barplot(pc4$eig/sum(pc4$eig)*100, 
        xlab = "PC",
        ylab = "Percentage of represented variation")

      
biplot.pop(pc4, gl_season[[4]], labels = "pop", xaxis = 1, yaxis = 2)

Plotting populations

2015.2016

#pc5 <- gl.pcoa(gl_season[[5]], nfactors = 5)
barplot(pc5$eig/sum(pc5$eig)*100, 
        xlab = "PC",
        ylab = "Percentage of represented variation")

      
biplot.pop(pc5, gl_season[[5]], labels = "pop", xaxis = 1, yaxis = 2)

Plotting populations

2016.2017

pc6 <- gl.pcoa(gl_season[[6]], nfactors = 5)

Performing a PCoA, individuals as entities, SNP loci as attributes
Ordination yielded 26 informative dimensions from 93 original dimensions
  PCoA Axis 1 explains 6.9 % of the total variance
  PCoA Axis 1 and 2 combined explain 13.2 % of the total variance
  PCoA Axis 1-3 combined explain 18 % of the total variance

barplot(pc6$eig/sum(pc6$eig)*100, 
        xlab = "PC",
        ylab = "Percentage of represented variation")

      
biplot.pop(pc6, gl_season[[6]], labels = "pop", xaxis = 1, yaxis = 2)

Plotting populations

2017.2018

pc7 <- gl.pcoa(gl_season[[7]], nfactors = 5)

Performing a PCoA, individuals as entities, SNP loci as attributes
Ordination yielded 31 informative dimensions from 117 original dimensions
  PCoA Axis 1 explains 6.5 % of the total variance
  PCoA Axis 1 and 2 combined explain 12.1 % of the total variance
  PCoA Axis 1-3 combined explain 16.7 % of the total variance

barplot(pc7$eig/sum(pc7$eig)*100, 
        xlab = "PC",
        ylab = "Percentage of represented variation")

      
biplot.pop(pc7, gl_season[[7]], labels = "pop", xaxis = 1, yaxis = 2)

Plotting populations