Collection of cassava landraces and associated farmers’ knowledge, genetic diversity and viral incidence assessment in western Kenya													
Ivan J Obare *1,2, Miriam K Charimbu2, Joseph Mafurah2, Christine K Mutoni3, Vincent W Woyengo4, Trushar Shah5 and Morag E Ferguson5 													
Supplementary file_S3: 'R' Scripts used in data analysis
													

#DartR


#install.packages("dartR")
library(dartR)
## I added this. Not sure if necessaray
 gl.install.vanilla.dartR()
setwd("C:/Users/ADMIN/Documents/R Workspace/Western Landraces")
getwd()
my_data<-gl.read.dart(filename ="Western landraces plus breeding lines.csv",ind.metafile ="Western landraces plus breeding lines _metadata.csv")


#To recalculate the metadata metrics
gl.recalc.metrics(my_data)

#To list the individuals
my_data@ind.names

#To list the markers / loci
my_data@loc.names

#To query number of individuals
nInd(my_data)

nPop(my_data)

levels(pop(my_data))
nlevels(pop(my_data))


######Filter based on different criteria - I added a new dataset after filtering
gl.report.monomorphs(my_data)
my_data1 <- gl.filter.monomorphs(my_data)
nLoc(my_data1)

#####To recalculate the metadata metrics
gl.recalc.metrics(my_data1)

nLoc(my_data1)


#Filter use different thresholds, works better when you specify a new file
my_data2 <- gl.filter.callrate(my_data1, method = "loc" , threshold = 0.99)
nLoc(my_data2)

####callrate is for SNP only, what to use for individuals??
my_data3 <- gl.filter.callrate(my_data2, method = "ind" , threshold = 0.95)
nInd(my_data3)



#####Calculate genetic distance
my_data3_distance_matrix<-gl.dist.ind(my_data3, method = "euclidean" )
my_data3_distance_matrix

distance_509_matrix<-as.matrix(my_data3_distance_matrix)
write.csv(distance_509_matrix,"509 distance matrix.csv")

#To generate the tree directly from the data, not distance matrix. You can also look at the distribution of distance measures. Use left and right arrows on plots tab.
#We tried this but can't add scale, so use GG tree
gl.tree.nj(my_data3, type ="phylogram")

#To generate the tree directly from the data, not distance matrix. You can also look at the distribution of distance measures. Use left and right arrows on plots tab.
gl.tree.nj(my_data3, type ="fan")
#This reads in the get duplicates script
source("GetDuplicate.R")
my_data3_distance_matrix
# to convert this DIS object to a matrix and divide everything by 100 so can use 0.95 or whatever cutoff. this is distance. we need to convert to similiraty by 1-
my_data3_matrix<-as.matrix(1-(my_data3_distance_matrix/100))
write.csv(my_data3_matrix,"ivan_datmatrix 509.csv")
#To identify row numbers un the distance matrix with known duplicates (TME14) the pipe means or. this is case sensitive

grep("TME-14|TME14|Abbey-Ife|TME14-Ug",rownames(my_data3_matrix))
##grep("NAROCASS-1|NAROCASS1 (TZ130)",rownames(my_data4_matrix))
##grep("Kibandameno forked kubwa|Kibandameno forked ndogo|Kibandameno straight",rownames(my_data3_matrix))

# to find smilarity
my_data3_matrix[c(93,130,142,154),c(93,130,142,154)]
#To select the most representative clone of the duplicates, we want to select the one that has the highest average similarity
#Enter each group duplicates, by copying and pasting from .csv, file then getting row numbers


grep("IJ20|IJ53|IJ101|IV126|IJ104|IJ112|IJ229|IJ09|IJ114|IJ195|IJ29|IJ43|IJ168|IJ91|IJ115|IJ144|IJ36|IJ79|IJ232|IJ13Fumbachai|IJ129|IJ15|IJ41|IJ158|IJ66|IJ231|IJ40|IJ89|IJ123|IJ75",rownames(my_data3_matrix))
Dupset1<-my_data3_matrix[c(6,8,10,11,19,36,59,92,106,131,137,158,160,173,190,199,207,211,216,223,224,234,235,239,244,245,246,249,257,258,268,271,325,371,404,427,428,429,437,439,443,453,475,478,484,489,490,497),c(6,8,10,11,19,36,59,92,106,131,137,158,160,173,190,199,207,211,216,223,224,234,235,239,244,245,246,249,257,258,268,271,325,371,404,427,428,429,437,439,443,453,475,478,484,489,490,497)]
summary(Dupset1)
grep("IV05|MM06/0082|IV181|IV109|IV187|MM96/9308",rownames(my_data3_matrix))
Dupset2<-my_data3_matrix[c(2,3,17,347,372,488),c(2,3,17,347,372,488)]
summary(Dupset2)
grep("IV30|IV33|IV31|IV54|IV56|IV57|IV80|IV34|IV25|IV53|IV48|IV28|IV24|MM98/1642|IV45|IV15",rownames(my_data3_matrix))
Dupset3<-my_data3_matrix[c(12,20,52,72,82,102,119,123,126,132,143,164,177,188,198,267,317,339,363,375,415,419,431,466,467,501),c(12,20,52,72,82,102,119,123,126,132,143,164,177,188,198,267,317,339,363,375,415,419,431,466,467,501)]
summary(Dupset3)
grep("IJ37|IJ57|IJ107|IJ58|IJ74|IJ152|IJ76|IJ131|IJ127|IJ35|IJ21|IJ153|IJ10|IJ98|IJ23|IJ46|IJ22|IJ24",rownames(my_data3_matrix))
Dupset4<-my_data3_matrix[c(18,23,33,56,68,89,91,109,116,120,136,171,189,192,201,208,212,224,232,234,236,244,252,256,259,268,270,321,325,327,328,331,332,336,340,343,346,350,351,352,355,356,357,359,362,369,370,371,374,376,379,384,386,387,388,390,393,395,396,397,404,406,408,414,455,479),c(18,23,33,56,68,89,91,109,116,120,136,171,189,192,201,208,212,224,232,234,236,244,252,256,259,268,270,321,325,327,328,331,332,336,340,343,346,350,351,352,355,356,357,359,362,369,370,371,374,376,379,384,386,387,388,390,393,395,396,397,404,406,408,414,455,479)]
summary(Dupset4)
grep("IJ42|IJ251|IJ06|IJ169|IJ39|IJ81",rownames(my_data3_matrix))
Dupset5<-my_data3_matrix[c( 5,24,47,161,329,464),c(5,24,47,161,329,464)]
summary(Dupset5)
grep("IJ02|IJ215|IJ03|IJ223|IJ252|IJ239|IJ214|IJ210|IJ14|IJ256|IJ82",rownames(my_data3_matrix))
Dupset6<-my_data3_matrix[c(7,25,37,172,195,207,209,219,231,233,243,254,265,276,350,351,356,365,387,409,455),c(7,25,37,172,195,207,209,219,231,233,243,254,265,276,350,351,356,365,387,409,455)]
summary(Dupset6)
grep("IV16|IJ63|IJ106|IJ102|IJ72|IJ62|IJ67|IJ125|IJ155|IJ103|IJ19|IJ136|IJ119|IJ124|IV10|IJ70|IJ161",rownames(my_data3_matrix))
Dupset7<-my_data3_matrix[c(30,32,88,145,147,157,179,180,191,212,251,253,256,258,259,262,281,293,296,299,302,308,315,320,322,344,367,368,380,391,398,401,402,424,426,438,440,448,452,461,472,476,484,486,488,494,505),c(30,32,88,145,147,157,179,180,191,212,251,253,256,258,259,262,281,293,296,299,302,308,315,320,322,344,367,368,380,391,398,401,402,424,426,438,440,448,452,461,472,476,484,486,488,494,505)]
summary(Dupset7)
grep("IV22|IJ211|IJ236|IJ207|IJ247|IJ235|IJ30|IV166|IV171|IV168|IJ212|IJ175|IV21",rownames(my_data3_matrix))
Dupset8<-my_data3_matrix[c(22,46,57,345,352,380,395,402,406,413,414,478,479,480),c(22,46,57,345,352,380,395,402,406,413,414,478,479,480)]
summary(Dupset8)
grep("IJ05|IJ31|IJ25|IJ248|IJ12|IJ242|IJ38|IJ27|IJ26|IJ213|IJ32|IJ194",rownames(my_data3_matrix))
Dupset9<-my_data3_matrix[c(34,35,45,58,61,69,77,80,192,203,204,215,216,227,239,251,262,273,328,329,342,365,376,378,386,399,409,472),c(34,35,45,58,61,69,77,80,192,203,204,215,216,227,239,251,262,273,328,329,342,365,376,378,386,399,409,472)]
summary(Dupset9)
grep("IJ113|IV39|IJ11|IV41|IJ111|IJ17|IJ28|IV194|IJ78|IV11|IJ162|IJ160",rownames(my_data3_matrix))
Dupset10<-my_data3_matrix[c(42,66,67,90,122,125,168,190,191,202,225,237,249,260,261,271,272,278,300,407,420,421,423,432,433,435,444,445,446,456,457,468,469,480,481,491,492,499,502,504),c(42,66,67,90,122,125,168,190,191,202,225,237,249,260,261,271,272,278,300,407,420,421,423,432,433,435,444,445,446,456,457,468,469,480,481,491,492,499,502,504)]
summary(Dupset10)
grep("IV68|IV36|IV141|IV130|IV134|IV81|IV185|IV43|IV143|IV201|IV09|IV184|IV196|IV70|IV195|IJ192|IV18|IV182|IV142|IV136|IV179|IV123|IV132|IV71|IV174|IV95|IV52|IV88|IV55|IV127",rownames(my_data3_matrix))
Dupset11<-my_data3_matrix[c(79,86,103,110,146,183,205,206,210,275,284,298,318,324,326,335,337,338,347,348,360,372,377,383,394,403,405,417,424,441,449,450,451,465,474,495,496),c(79,86,103,110,146,183,205,206,210,275,284,298,318,324,326,335,337,338,347,348,360,372,377,383,394,403,405,417,424,441,449,450,451,465,474,495,496)]
summary(Dupset11)
grep("IV203|IV150|IV04|IV176|IV204|IV85|IV106|IV154|IV12|IV93|MM96/7151	",rownames(my_data3_matrix))
Dupset12<-my_data3_matrix[c(84,87,247,279,330,334,400,425,436,437,440,447,449,460,466,467,473,485,495,506),c(84,87,247,279,330,334,400,425,436,437,440,447,449,460,466,467,473,485,495,506)]
summary(Dupset12)

grep("IJ51|IJ86|IJ48|IJ95|IJ92|IJ33|IJ230|IJ88|IJ128",rownames(my_data3_matrix))
Dupset13<-my_data3_matrix[c(70,97,117,138,139,152,162,204,359),c(70,97,117,138,139,152,162,204,359)]
summary(Dupset13)

grep("IV29|IJ99|IV158|IV163|IV44|IV189|IV32|IJ121|IV50|IV35|IV63|IV198|IV115|IJ118|IV190|IV175|IV19|IV100|IV209|IV65|IV17",rownames(my_data3_matrix))
Dupset14<-my_data3_matrix[c(9,44,98,99,148,169,174,215,220,238,240,272,293,318,319,323,324,333,334,337,338,341,345,349,358,361,363,364,373,381,382,385,392,401,403,407,412,492),c(9,44,98,99,148,169,174,215,220,238,240,272,293,318,319,323,324,333,334,337,338,341,345,349,358,361,363,364,373,381,382,385,392,401,403,407,412,492)]
summary(Dupset14)

grep("IJ52|IJ61|IJ59|IJ45|IJ68",rownames(my_data3_matrix))
Dupset15<-my_data3_matrix[c(101,105,108,144,167),c(101,105,108,144,167)]
summary(Dupset15)

grep("IV37|IV87|IJ34|IJ16",rownames(my_data3_matrix))
Dupset16<-my_data3_matrix[c(43,81,121,278,286,294,296,300,301,307,311,416,428,464),c(43,81,121,278,286,294,296,300,301,307,311,416,428,464)]
summary(Dupset16)

grep("IJ73|IV159|IJ71|IJ01|IV124|IV49|IV137|IV139|IV01",rownames(my_data3_matrix))
Dupset17<-my_data3_matrix[c(1,13,113,124,170,375,463,498,506),c(1,13,113,124,170,375,463,498,506)]
summary(Dupset17)

grep("IV77|IV151|IV153|IV180|IV192|Abbey-Ife|IV205|IV199|TME14-Ug|IJ171|IV188|IV164|IV169|IV165|IV206|IV76|IV152|IV03|IV197|IV193|IV167|IV144|IV202|IV177|TME-14|TME14",rownames(my_data3_matrix))
Dupset18<-my_data3_matrix[c(73,93,130,142,154,185,197,322,326,335,341,344,354,358,364,366,368,373,385,389,391,419,431,432,477,501),c(73,93,130,142,154,185,197,322,326,335,341,344,354,358,364,366,368,373,385,389,391,419,431,432,477,501)]
summary(Dupset18)

grep("IV107|IV191|IV138|IV145|IV210|IV26|IV135|IV86|IV104|IV46|IV98",rownames(my_data3_matrix))
Dupset19<-my_data3_matrix[c(48,135,280,288,315,361,413,452,487,500,508),c(48,135,280,288,315,361,413,452,487,500,508)]
summary(Dupset19)

grep("IJ94|IV02|IV51|IJ122|IV23",rownames(my_data3_matrix))
Dupset20<-my_data3_matrix[c(49,71,140,181,227),c(49,71,140,181,227)]
summary(Dupset20)

grep("IJ87|IJ85",rownames(my_data3_matrix))
Dupset21<-my_data3_matrix[c(115,150),c(115,150)]
summary(Dupset21)

grep("IJ49|MM96/4684|IJ50",rownames(my_data3_matrix))
Dupset22<-my_data3_matrix[c(75,151,163),c(75,151,163)]
summary(Dupset22)

grep("IV58|IV59|IV133|IV66|IV75|IV117|IV69|IV79|IV101|IV72|IV128|IV83|IV183|IV116|IV122|IV140|IV13|IV121|MH95/0183|IV99|IV186|IV90|IV97|IV89",rownames(my_data3_matrix))
Dupset23<-my_data3_matrix[c(28,31,186,200,218,248,263,264,269,277,287,290,299,306,314,360,394,436,446,450,451,460,462,463,473,474,487,496,498,504,507,508,509),c(28,31,186,200,218,248,263,264,269,277,287,290,299,306,314,360,394,436,446,450,451,460,462,463,473,474,487,496,498,504,507,508,509)]
summary(Dupset23)

grep("IJ120|IJ108|IJ166|IV91",rownames(my_data3_matrix))
Dupset24<-my_data3_matrix[c(189,203,309,311),c(189,203,309,311)]
summary(Dupset24)

grep("IJ148|IJ156|IJ132|IJ116|IV84|IJ204|IJ150|IJ165|IJ163",rownames(my_data3_matrix))
Dupset25<-my_data3_matrix[c(187,199,202,209,245,274,286,307,429),c(187,199,202,209,245,274,286,307,429)]
summary(Dupset25)

grep("IV61|IV73|MM96/2480	",rownames(my_data3_matrix))
Dupset26<-my_data3_matrix[c(214,255),c(214,255)]
summary(Dupset26)

grep("IV62|MM96/4271(NASE 14)|IV74|IV60|IV94",rownames(my_data3_matrix))
Dupset27<-my_data3_matrix[c(213,226,266,292),c(213,226,266,292)]
summary(Dupset27)

grep("IJ100|IJ219|IJ216|IJ218|IJ97|IJ217|IJ105|IJ180|IJ177|IJ142|IJ154|IJ183|IJ221|IJ250|IJ255|IJ179",rownames(my_data3_matrix))
Dupset28<-my_data3_matrix[c(196,232,236,246,276,321,331,343,355,378,379,399,458,469,502,503),c(196,232,236,246,276,321,331,343,355,378,379,399,458,469,502,503)]
summary(Dupset28)

grep("IJ139|IJ137|IJ135|IJ134",rownames(my_data3_matrix))
Dupset29<-my_data3_matrix[c(194,229,241,242),c(194,229,241,242)]
summary(Dupset29)

grep("IJ109|IJ07|IJ138|IJ140|IJ141",rownames(my_data3_matrix))
Dupset30<-my_data3_matrix[c(53,201,230,254,265),c(53,201,230,254,265)]
summary(Dupset30)

grep("IV208|IV148|IV06|IV172|IV155|IV173|IV156|GG52_Seruruseke_MM96/5280",rownames(my_data3_matrix))
Dupset31<-my_data3_matrix[c(29,297,317,381,392,411,415,442),c(29,297,317,381,392,411,415,442)]
summary(Dupset31)

grep("IJ225|IJ202|IJ203|IJ197|IV170|IJ167|IJ201|IJ200|IJ254|IJ243|IJ198|IJ176|IJ226|IJ206|IJ199",rownames(my_data3_matrix))
Dupset32<-my_data3_matrix[c(333,342,346,369,397,416,426,427,438,439,475,486,489,491,497),c(333,342,346,369,397,416,426,427,438,439,475,486,489,491,497)]
summary(Dupset32)

grep("IV14|IV178|IJ227|IJ18",rownames(my_data3_matrix))
Dupset33<-my_data3_matrix[c(55,78,370,382,417,418,422,430,434,441,442,454,458,459,465,470,471,477,482,483,493,500,503,509),c(55,78,370,382,417,418,422,430,434,441,442,454,458,459,465,470,471,477,482,483,493,500,503,509)]
summary(Dupset33)

grep("IJ170|IJ224|IJ237|IJ209",rownames(my_data3_matrix))
Dupset34<-my_data3_matrix[c(357,396,420,443),c(357,396,420,443)]
summary(Dupset34)

grep("IV207|MIGYERA|IV120|IJ126|IV82",rownames(my_data3_matrix))
Dupset35<-my_data3_matrix[c(74,222,273,410,447),c(74,222,273,410,447)]
summary(Dupset35)

grep("IV119|IJ190|IJ181|IJ193|IJ178|IJ186|IJ184",rownames(my_data3_matrix))
Dupset36<-my_data3_matrix[c(422,435,448,457,470,493,494),c(422,435,448,457,470,493,494)]
summary(Dupset36)

grep("IV113|IV111|IV147",rownames(my_data3_matrix))
Dupset37<-my_data3_matrix[c(421,430,445),c(421,430,445)]
summary(Dupset37)

grep("IJ173|IJ205|IJ228",rownames(my_data3_matrix))
Dupset38<-my_data3_matrix[c(393,453,456),c(393,453,456)]
summary(Dupset38)

grep("IV108|MM98/3567|IJ187|IV129",rownames(my_data3_matrix))
Dupset39<-my_data3_matrix[c(4,459,476,485),c(4,459,476,485)]
summary(Dupset39)

grep("IV114|IV112",rownames(my_data3_matrix))
Dupset40<-my_data3_matrix[c(433,481),c(433,481)]
summary(Dupset40)

grep("IJ189|IV146",rownames(my_data3_matrix))
Dupset41<-my_data3_matrix[c(418,483),c(418,483)]
summary(Dupset41)

grep("MM06/0083|MM06/0143",rownames(my_data3_matrix))
Dupset42<-my_data3_matrix[c(15,63),c(15,63)]
summary(Dupset42)

grep("IV47|IV40",rownames(my_data3_matrix))
Dupset43<-my_data3_matrix[c(159,178),c(159,178)]
summary(Dupset43)

grep("IJ08|IJ110",rownames(my_data3_matrix))
Dupset44<-my_data3_matrix[c(76,225),c(76,225)]
summary(Dupset44)

grep("IJ83|IJ130",rownames(my_data3_matrix))
Dupset45<-my_data3_matrix[c(182,228),c(182,228)]
summary(Dupset45)

grep("IV78|IJ146",rownames(my_data3_matrix))
Dupset46<-my_data3_matrix[c(221,231),c(221,231)]
summary(Dupset46)

grep("MM06/0138|IV64",rownames(my_data3_matrix))
Dupset47<-my_data3_matrix[c(50,250),c(50,250)]
summary(Dupset47)

grep("IV200|IJ208|IJ234",rownames(my_data3_matrix))
Dupset48<-my_data3_matrix[c(353,384,490),c(353,384,490)]
summary(Dupset48)

grep("98/3875|IV149",rownames(my_data3_matrix))
Dupset49<-my_data3_matrix[c(312,454),c(312,454)]
summary(Dupset49)

grep("IJ172|IJ188",rownames(my_data3_matrix))
Dupset50<-my_data3_matrix[c( 444,471),c(444,471)]
summary(Dupset50)

grep("MM06/0139|IV110",rownames(my_data3_matrix))
Dupset51<-my_data3_matrix[c(26,499),c(26,499)]
summary(Dupset51)

#To identify duplicates in the dataset based on threshold of 0.9
install.packages("caret")
library(caret)

duplicateslist <-GetDuplicate(my_data3_matrix,0.80,F)
duplicateslist

#To write this list out to a file. You will need to sort in Excel. You can do a left to right sort in Excel. Data, Sort, OPtions, Left
write.csv(duplicateslist,file="Duplicates_509.csv",na="")


#To look at a table of categories
table(pop(my_data3))

#CONSTRUCT TREE IN GGTREE
#read in the label data to pass it to the tree
cassava_labels<-read.csv("./Western landraces plus breeding lines _metadata.csv", header = TRUE)
#This shows you its a list
typeof(cassava_labels)
cassava_labels

#Now let's construct the tree in phyclust
install.packages("phyclust")
library(phyclust)
#convert from matrix back to distance object type
IBS.num.dist<-as.dist(my_data3_distance_matrix) # QC0
#constructing tree using NJ / wards clustering method
NJ.IBS<-nj(IBS.num.dist) # QC0
HC.IBS<-hclust(IBS.num.dist, method = "average") # QC0
HC.IBS

#Now the plotting with ggtree
BiocManager::install("ggtree")
install.packages(ggtree)
library(ggtree)
#converting the cluster file (hclust) into a phylogenetic object that ggtree can recognize
phyfile<-as.phylo(NJ.IBS)
#phyfile<-as.phylo.hclust(HC.IBS)
phyfile$tip.label

#To get cladogram (the tree is p)
p <-ggtree(phyfile, branch.length = 'none', size = 0.1) + layout_dendrogram()
p

p$data$label

#To see what your dataset looks like
p$data


#To export the plots as PNG
pdf("cladogram_ivan11.pdf", width = 11, height = 8)
p + geom_tiplab(linesize = 0.05,size=0.2, aes(label = label, angle = 270)) + geom_treescale(width = 1, linesize = 2) + geom_vline(aes(xintercept = -21), colour = "Red")

#,col = c("blue","brown","magenta"))
dev.off()
BiocManager::install(c("SNPRelate", "qvalue"))
BiocManager::install("ggtree")

install.packages("dartR")
library(dartR)
## I added this. Not sure if necessaray
gl.install.vanilla.dartR()
getwd()
setwd("C:/Users/ADMIN/Documents/R Workspace/Western Landraces")



#Read in file from working directory (do not include Metadata from DArT for now). The metadata should be in a csv format with id and pop columns. This is required even if no pop. 
my_data <- gl.read.dart(filename = "DAPC Western plus reference.csv", ind.metafile = "DAPC Western plus reference_metadata.csv")

#To recalculate the metadata metrics
gl.recalc.metrics(my_data)

#To list the individuals
my_data@ind.names

#To list the markers / loci
my_data@loc.names

#To query number of individuals
nInd(my_data)

nPop(my_data)

levels(pop(my_data))
nlevels(pop(my_data))


######Filter based on different criteria - I added a new dataset after filtering
gl.report.monomorphs(my_data)
my_data1 <- gl.filter.monomorphs(my_data)
nLoc(my_data1)


#####To recalculate the metadata metrics
gl.recalc.metrics(my_data1)

nLoc(my_data1)


#Filter use different thresholds, works better when you specify a new file
my_data2 <- gl.filter.callrate(my_data1, method = "loc" , threshold = 0.99)
nLoc(my_data2)

####callrate is for SNP only, what to use for individuals??
my_data3 <- gl.filter.callrate(my_data2, method = "ind" , threshold = 0.95)
nInd(my_data3)

my_data3

my_data4_distance_matrix <- gl.dist.ind(my_data3, method = 'euclidean')

library(devtools)
library(dartR)
my_data4 <-  gl.dist.ind(my_data3)
my_data5 <- gl.pcoa(my_data4,nfactors=2)
output <- gl.pcoa.plot(my_data5, my_data3)

install.packages("ggplot2")
library(ggplot2)
library(pca3d)
library(dartR)

pca <- prcomp(my_data4, scale.=TRUE)
gr <- factor(my_data4)
summary(gr)

pca2d(pca, group = gr , biplot=TRUE, biplot.vars=3)

pop(my_data4)


pcoa_obj<-gl.pcoa(my_data3)
gl.pcoa.plot(pcoa_obj,my_data3)
pcoa_obj



#to do the DAPC we need to convert genlight to genind
#You will have to enter the number of PCs to retain in Q3 (consol). First look at variance explained. the lower the number you put in, the closer together your populations will be. The higher the value, the more discrimination you'll get. I used 50 and 8. Also enter number of eigenvalues.I put 8.
gm_genind<-gl2gi(my_data3)
dapc_gm_genind<-dapc(gm_genind)



scatter(dapc_gm_genind, legend = TRUE, clabel = 0.5, posi.leg = "topright", scree.da = FALSE)
#To export the plot as a high res PNG file; firstly name of file and size and res; then plot and colour details and the produce the plot
png("dapc_western.png", width = 5, height = 5, units = "in", res = 600)
scatter(dapc_gm_genind, legend = TRUE, clabel = 0.5, posi.leg = "topright",posi.da = "topleft",ratio.da = 0.2, cleg = 0.75) 
#,col = c("blue","brown","magenta"))
dev.off()

#Have a look at hierfstat
library(hierfstat)

#Now calculate the Stats using hierfstat program
#Converts genind objects from adegenet into a hierfstat data frame
Mydata2<-genind2hierfstat(gm_genind)
Mydata2$pop

#WC84 can calculate a distance among populations (ie Fst); we're doing two methods (Weir and Cockerham table of Fst stats)
wcdist<-genet.dist(gm_genind, method = "WC84")
nei_dist<-genet.dist(gm_genind, method = "Nei87")
wcdist
nei_dist
#Converting it to matrix and Export WC matrix into .csv file
wcdistmatrix<-as.matrix(wcdist)
nei_distmatrix<-as.matrix(nei_dist)

# just put upper triangle blank
wcdistmatrix[upper.tri(wcdistmatrix,diag = TRUE)]<-""
write.csv(wcdistmatrix,"Fst_WC_567ind_8pops.csv")


neidistmatrix<-as.matrix(nei_dist)
neidistmatrix[upper.tri(neidistmatrix,diag = TRUE)]<-""
write.csv(neidistmatrix, "Nei_567ind_8pops.csv")

#Basic stats overall and by locus, not by population
Mydata2_basic_stats<-hierfstat::basic.stats(data = Mydata2)
Mydata2_basic_stats
Mydata2_basic_stats$overall

#Install popgenreport
library(PopGenReport)
library(tinytex)

Mydata2_popgen<-popgenreport(gm_genind, mk.pdf = FALSE)
Mydata2_popgen
pop.freq(Mydata2)



gm_genind
popNames(gm_genind)

summary(Mydata2$pop)
dim(Mydata2[Mydata2$pop=="1",])
Mydata2[Mydata2$pop=="1",1]
#To get population level stats (try changing to Pop1 etc to run)
basicstat_Pop1 <- basic.stats(Mydata2[Mydata2$pop=="1",], diploid = TRUE)
basicstat_Pop2 <- basic.stats(Mydata2[Mydata2$pop=="2",], diploid = TRUE)
basicstat_Pop3 <- basic.stats(Mydata2[Mydata2$pop=="3",], diploid = TRUE)
# for overall
basicstat_overall <- basic.stats(Mydata2, diploid = TRUE)

basicstat_overall$overall