# Script for pedigree reconstruction utilized for Piza Roca et al. 2019 "Fitness benefits of male dominance behaviours depend on the degree of inbreeding in a polyandrous lizard" # Method described in: Huisman, J. 2017 Pedigree reconstruction from SNP data: Parentage assignment, sibship clustering and beyond. Molecular Ecology Resources 17, 1009-1024. (doi:10.1111/1755-0998.12665) library(sequoia) Geno <- GenoConvert(InFile = "C:/Users/carme/Desktop/WORK/SNP data/generating sequoia inputs 180501/inputfile_sequoia_maf038_ld034.raw") #convert PLINK output into sequoia genome file Geno <- GenoConvert(InFile = "File path/349SNPs_pedigree_reconstruction.raw") #convert PLINK output into sequoia genome file str(Geno) LH <- read.csv("LH_sequoia.csv") #read in life history data # Run sequoia - duplicate check & parentage assignment only # (maximum number of sibship-clustering iterations = 0) ParOUT <- sequoia(GenoM = Geno, LifeHistData = LH[c("ID", "Sex", "BY")], MaxSibIter = 0) names(ParOUT) # Run sequoia - sibship clustering & grandparent assignment # use parents assigned above (in 'ParOUT$PedigreePar') SeqOUT <- sequoia(GenoM = Geno, SeqList = ParOUT, MaxSibIter = 20) Pedigree <- SeqOUT$Pedigree # Save the pedigree Pedigree$dam <- gsub("F00", NA, Pedigree$dam) #get rid of dummy mums Pedigree$sire <- gsub("M00", NA, Pedigree$sire) #get rid of dummy dads #LH <- subset(LH, !ID%in%c("Blondie", "Henry", "Jungle", "Mal", "Scruffy")) hatchlings <- LH[95:140,]$ID # calculate number of offspring dads <- setNames(as.data.frame(table(Pedigree$sire)), c("ID", "Paternities")) unsuccessful.dads <- data.frame(ID = unique(subset(LH, Sex == 2)$ID[which(!subset(LH, Sex == 2)$ID%in%dads$ID)]), Paternities = 0) dads <- rbind(dads, unsuccessful.dads) head(dads) tail(dads) # calculate the number of adult offspring offs <- setNames(as.data.frame(with(subset(Pedigree, !id%in%hatchlings), table(sire))), c("ID", "Adult.Offs")) dads<- merge(dads, offs, by = "ID", all.x = TRUE) dads$Adult.Offs <- ifelse(is.na(dads$Adult.Offs), 0, dads$Adult.Offs) head(dads) sum(dads$Paternities) #106 sum(dads$Adult.Offs) #65 #add PCA morph+dom morph <- read.csv("Morph.csv", stringsAsFactors = FALSE) morph <- morph[,c("ID", "TORSO.CMS", "JAW.WIDTH.MMS", "WEIGHT.GMS")] #only IDs and interesting morph morph <- na.omit(morph) #remove NAs, gets rid of all incomplete samples! morph <- morph[!duplicated(morph$ID),] #get only the newest morphology from each individual dom <- read.csv("Dominance.csv") dat <- merge(morph, dom) head(dat) pca <- princomp(dat[2:5], cor = TRUE) summary(pca) loadings(pca) # Dfreq is negatively loaded in PC2 md <- setNames(as.data.frame(cbind(dat$ID, pca$scores[,1], pca$scores[,2]*-1)), c("ID", "Dom.PC1", "Dom.PC2")) #merge PC1 and PC2 with IDs dat2 <- merge(dads, md, by = "ID") subset(dads,!dads$ID%in%dat2$ID) #Only 6 individuals removed (1 adult paternity) with no morph data with(subset(dat2, !is.na(Dom.PC1)), sum(Paternities)) #104 with(subset(dat2, !is.na(Dom.PC1)), sum(Adult.Offs)) #64 with(subset(dat2, !is.na(Dom.PC1)), length(unique(ID))) #111 #add heterozigozyty hetzy <- read.csv("inbreeding.csv") dat2 <- merge(dat2, hetzy[,c("IID", "Fhat3")], by.x = "ID", by.y = "IID") #add sexually maturating year, departure year and number of sightings dat2 <- merge(dat2, LH, by = "ID") dat2$Sex <- ifelse(dat2$Sex == 1, "Fe", "M") dat2$repr.time <- (dat2$Departure - dat2$BY)+1 # calculate sexually mature time dat2 <- subset(dat2, N_sight >= 15) # Restrict to minimum of 15 sightings head(dat2) write.csv(dat2, "Paternity_master.csv", row.names = FALSE)