allopolyVignette.R
In polysat: Tools for Polyploid Microsatellite Analysis

### R code from vignette source 'allopolyVignette.Rnw'

###################################################
### code chunk number 1: allopolyVignette.Rnw:129-134
###################################################
library(polysat)
data(AllopolyTutorialData)
summary(AllopolyTutorialData)
# make a copy of the dataset to modify
mydata <- AllopolyTutorialData


###################################################
### code chunk number 2: allopolyVignette.Rnw:151-161
###################################################
# Calculate the length of each genotype vector (= the number of alleles) and
# construct a TRUE/FALSE matrix of whether that number is greater than four.
tooManyAlleles <- apply(Genotypes(mydata), c(1,2), 
                        function(x) length(x[[1]])) > 4
# Find position(s) in the matrix that are TRUE.
which(tooManyAlleles, arr.ind=TRUE) # 43rd sample, second locus
# Look at the identified genotype, then replace it with missing data.
Genotype(mydata, 43, 2)
Genotype(mydata, 43, 2) <- Missing(mydata)
Genotype(mydata, 43, 2)


###################################################
### code chunk number 3: allopolyVignette.Rnw:168-170 (eval = FALSE)
###################################################
## mydist <- meandistance.matrix(mydata, distmetric=Lynch.distance,
##                               progress=FALSE)


###################################################
### code chunk number 4: allopolyVignette.Rnw:173-174
###################################################
load("vignettebuild/AllopolyTutorialDist.RData")


###################################################
### code chunk number 5: allopolyVignette.Rnw:177-180
###################################################
require(ape)
mynj <- nj(mydist)
plot(mynj, type="unrooted")


###################################################
### code chunk number 6: allopolyVignette.Rnw:190-191
###################################################
mydata <- deleteSamples(mydata, c("301","302","303"))


###################################################
### code chunk number 7: allopolyVignette.Rnw:197-201
###################################################
par(mfrow=c(2,1))
mypca <- cmdscale(mydist[Samples(mydata), Samples(mydata)])
plot(mypca[,1], mypca[,2])
hist(mypca[,1], breaks=30)


###################################################
### code chunk number 8: allopolyVignette.Rnw:210-214
###################################################
pop1ind <- Samples(mydata)[mypca[,1] <= 0]
pop2ind <- Samples(mydata)[mypca[,1] > 0]
PopInfo(mydata)[pop1ind] <- 1
PopInfo(mydata)[pop2ind] <- 2


###################################################
### code chunk number 9: allopolyVignette.Rnw:274-276 (eval = FALSE)
###################################################
## myassign <- processDatasetAllo(mydata, n.subgen = 2, SGploidy = 2, 
##                                usePops = TRUE, R = 500)


###################################################
### code chunk number 10: allopolyVignette.Rnw:279-280
###################################################
load("vignettebuild/AllopolyTutorialAssign.RData")


###################################################
### code chunk number 11: allopolyVignette.Rnw:298-300
###################################################
myassign$AlCorrArray[["Loc6", 1]]$significant.pos
myassign$AlCorrArray[["Loc6", 2]]$significant.pos


###################################################
### code chunk number 12: allopolyVignette.Rnw:310-311
###################################################
mydata <- deleteLoci(mydata, loci="Loc6")


###################################################
### code chunk number 13: allopolyVignette.Rnw:327-329
###################################################
par(mfrow = c(1,1))
plotSSAllo(myassign$AlCorrArray)


###################################################
### code chunk number 14: allopolyVignette.Rnw:356-358
###################################################
heatmap(myassign$AlCorrArray[["Loc5", "Pop1"]]$heatmap.dist,
        main = "Loc5, Pop1")


###################################################
### code chunk number 15: allopolyVignette.Rnw:363-371
###################################################
# A plot to show how the colors correspond to p-values in the
# heat map; you can repeat this for the other heat maps in this
# tutorial if you wish.
plot(x=seq(min(myassign$AlCorrArray[["Loc5", "Pop1"]]$heatmap.dist),
           max(myassign$AlCorrArray[["Loc5", "Pop1"]]$heatmap.dist), 
           length.out=12),
     y=rep(1,12), xlab="P-values", ylab="", bg=heat.colors(12), 
     pch=22, cex=3)


###################################################
### code chunk number 16: allopolyVignette.Rnw:373-375
###################################################
heatmap(myassign$AlCorrArray[["Loc6", "Pop1"]]$heatmap.dist,
        main = "Loc6, Pop1")


###################################################
### code chunk number 17: allopolyVignette.Rnw:380-382
###################################################
heatmap(myassign$AlCorrArray[["Loc7", "Pop1"]]$heatmap.dist,
        main = "Loc7, Pop1")


###################################################
### code chunk number 18: allopolyVignette.Rnw:384-386
###################################################
heatmap(myassign$AlCorrArray[["Loc7", "Pop2"]]$heatmap.dist,
        main = "Loc7, Pop2")


###################################################
### code chunk number 19: allopolyVignette.Rnw:410-412
###################################################
plotParamHeatmap(myassign$propHomoplasious, popname = "Pop1", 
                 main = "Proportion homoplasious loci:")


###################################################
### code chunk number 20: allopolyVignette.Rnw:423-425
###################################################
plotParamHeatmap(myassign$propHomoplasious, popname = "Pop2", 
                 main = "Proportion homoplasious loci:")


###################################################
### code chunk number 21: allopolyVignette.Rnw:437-439
###################################################
plotParamHeatmap(myassign$propHomoplMerged, popname = "Merged across populations", 
                 main = "Proportion homoplasious loci:")


###################################################
### code chunk number 22: allopolyVignette.Rnw:461-463
###################################################
plotParamHeatmap(myassign$missRate, popname = "All Individuals", 
                 main = "Missing data after recoding:")


###################################################
### code chunk number 23: allopolyVignette.Rnw:481-484
###################################################
myBestAssign <- myassign$bestAssign
myBestAssign
myBestAssign <- myBestAssign[1:5]


###################################################
### code chunk number 24: allopolyVignette.Rnw:496-497
###################################################
apply(myassign$missRate, 1, min)


###################################################
### code chunk number 25: allopolyVignette.Rnw:507-511
###################################################
myassign$AlCorrArray[["Loc1", "Pop1"]]$Kmeans.groups
myassign$AlCorrArray[["Loc1", "Pop2"]]$Kmeans.groups
myassign$AlCorrArray[["Loc1", "Pop1"]]$UPGMA.groups
myassign$AlCorrArray[["Loc1", "Pop2"]]$UPGMA.groups


###################################################
### code chunk number 26: allopolyVignette.Rnw:520-532
###################################################
myassign$TAGarray[["Loc1", "Pop1", 1]]$assignments
myassign$TAGarray[["Loc1", "Pop2", 1]]$assignments
myassign$mergedAssignments[["Loc1", 1]]$assignments
myassign$TAGarray[["Loc1", "Pop1", 2]]$assignments
myassign$TAGarray[["Loc1", "Pop2", 2]]$assignments
myassign$mergedAssignments[["Loc1", 2]]$assignments
myassign$TAGarray[["Loc1", "Pop1", 3]]$assignments
myassign$TAGarray[["Loc1", "Pop2", 3]]$assignments
myassign$mergedAssignments[["Loc1", 3]]$assignments
myassign$TAGarray[["Loc1", "Pop1", 4]]$assignments
myassign$TAGarray[["Loc1", "Pop2", 4]]$assignments
myassign$mergedAssignments[["Loc1", 4]]$assignments


###################################################
### code chunk number 27: allopolyVignette.Rnw:551-554
###################################################
corrLoc1AllInd <- alleleCorrelations(mydata, locus = "Loc1", n.subgen = 2)
corrLoc1AllInd$Kmeans.groups
corrLoc1AllInd$UPGMA.groups


###################################################
### code chunk number 28: allopolyVignette.Rnw:565-573
###################################################
TaLoc1.param2 <- testAlGroups(mydata, corrLoc1AllInd, SGploidy = 2,
                              null.weight = 0.5, tolerance = 0.05, 
                              swap = FALSE)
TaLoc1.param5 <- testAlGroups(mydata, corrLoc1AllInd, SGploidy = 2,
                              null.weight = 0.5, tolerance = 0.1, 
                              swap = FALSE)
TaLoc1.param2$assignments
TaLoc1.param5$assignments


###################################################
### code chunk number 29: allopolyVignette.Rnw:581-582
###################################################
myBestAssign[[1]] <- TaLoc1.param5


###################################################
### code chunk number 30: allopolyVignette.Rnw:590-592
###################################################
recodedData <- recodeAllopoly(mydata, myBestAssign)
summary(recodedData)


###################################################
### code chunk number 31: allopolyVignette.Rnw:600-604
###################################################
for(L in Loci(recodedData)){
  proportionmissing <- mean(isMissing(recodedData, loci=L))
  cat(paste(L,":",proportionmissing,"missing"),sep="\n")
}


###################################################
### code chunk number 32: allopolyVignette.Rnw:612-613
###################################################
table(Ploidies(recodedData))


###################################################
### code chunk number 33: allopolyVignette.Rnw:618-619
###################################################
Ploidies(recodedData)[,"Loc7"] <- 4


###################################################
### code chunk number 34: allopolyVignette.Rnw:627-630
###################################################
myfreq <- simpleFreq(recodedData)
myGst <- calcPopDiff(myfreq, metric = "Gst")
myGst


###################################################
### code chunk number 35: allopolyVignette.Rnw:635-636
###################################################
write.GeneMapper(recodedData, file = "tutorialRecodedData.txt")


###################################################
### code chunk number 36: allopolyVignette.Rnw:648-656
###################################################
catResults <- list()
length(catResults) <- length(Loci(mydata))
names(catResults) <- Loci(mydata)

for(L in Loci(mydata)){
  cat(L, sep="\n")
  catResults[[L]] <- catalanAlleles(mydata, locus=L, verbose=TRUE)
}