1.R包的下载(R version 3.6)
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("LEA")
如果R版本是老版本的，请查看网站Bioconductor release的要求进行R包下载。
注：具体怎么在R中进行操作在以下网站中有详细说明。
http://membres-timc.imag.fr/Olivier.Francois/tutoRstructure.pdf
2.加载其他函数
source("http://membres-timc.imag.fr/Olivier.Francois/Conversion.R")
source("http://membres-timc.imag.fr/Olivier.Francois/POPSutilities.R")
3导入输入文件
#struct2geno函数的讲解
struct2geno(file = input.file, TESS = FALSE, diploid = TRUE, FORMAT = 2,
extra.row = 0, extra.col = 0, output = "./genotype.geno")
注：input.file(导入文件):structure或tess 2.3格式到LEA所使用的.geno和.lfmm格式。
TESS=FALSE:当额外的列包含所有不包含基因型信息的列时候为FALSE.
TESS=TRUE:当额外的列包含所有不包含基因型信息和地理信息的列时候为FALSE.
diploid=TRUE:代表为二倍体。
FORMAT=2:意味着每个个体使用两行数据对标记进行编码。
FORMAT=1:意味着每个个体使用一行数据对标记进行编码。
extra.row/.col:表示额外行/列数。
"./genotype.geno"：输出文件的路径。
4.使用R包进行群体结构分析
4.1example1：Allelic markers(等位基因标记)
分析在欧洲二次contract后通过空间显式结合模拟所产生的群体遗传数据。在最初的一个分离的阶段后，一个物种开始从遥远的南方殖民到欧洲，一些在伊利比亚群岛，一些在土耳其。第二次contract发生在中世纪的欧洲，靠近德国。数据由10个二倍体个体的60个族群样本组成，用100个等位基因标记对10个二倍体个体进行基因分型。
4.1.1数据的下载和转换
input.file = "http://membres-timc.imag.fr/Olivier.Francois/secondary_contact.str"
struct2geno(file = input.file, TESS = TRUE, diploid = TRUE, FORMAT = 2,
extra.row = 0, extra.col = 0, output = "secondary_contact.geno")
4.1.2k值设置
k值的设置可以通过LEA包中的snmf函数设置完成。
library(LEA)
obj.snmf = snmf("secondary_contact.geno", K = 3, alpha = 100, project = "new")
qmatrix = Q(obj.snmf, K = 3)
4.1.3 600个个体祖先系数柱状图 {Q-matrix(用barplot呈现)}
barplot(t(qmatrix), col = c("orange","violet","lightgreen"), border = NA, space = 0,
xlab = "Individuals", ylab = "Admixture coefficients")
4.1.4 通过在欧洲地理地图上叠加饼图法进行admixture群体估计
为达到这个目标，我们需要读取样本的地理坐标和创建样本标志符。
共有60个群体样本，10个在每个群体中。
coord = read.table("coordinates.coord")
pop = rep(1:60, each = 10)
K = 3
Npop = length(unique(pop))
qpop = matrix(NA, ncol = K, nrow = Npop)
coord.pop = matrix(NA, ncol = 2, nrow = Npop)
for (i in unique(pop)){
qpop[i,] = apply(qmatrix[pop == i,], 2, mean)
coord.pop[i,] = apply(coord[pop == i,], 2, mean)}
library(mapplots)
plot(coord, xlab = "Longitude", ylab = "Latitude", type = "n")
map(add = T, col = "grey90", fill = TRUE)
for (i in 1:Npop){
add.pie(z = qpop[i,], x = coord.pop[i,1], y = coord.pop[i,2], labels = "",
col = c("orange","violet","lightgreen"))}
pop = scan("mypop.txt")
4.1.5 选择集群数
在LEA里选择集群数是基于交叉熵准则，这一准则也被admixture程序所使用。交叉熵准则基于一部分被掩盖基因型的预测和交叉验证的方法。更小的交叉熵准则的值通常意味着能够更好的运行。我们执行8个K值，然后选择交叉熵的曲线到达平台期时候的K值。
obj.snmf = snmf("/Users/bcl/secondary_contact.geno", K = 1:8, ploidy = 2, entropy = T,
alpha = 100, project = "new")
plot(obj.snmf, col = "blue4", cex = 1.4, pch = 19)
4.2example2:SNP data
接下来我们考虑来自于欧洲的拟南芥的SNP数据.这些数据是从Horton等人所研究的大量数据中抽出来的一部分。调查了一号染色体上1096个近交生态型(52001个SNPs)的群体结构。
4.2.1SNP数据的下载
数据保存在磁盘中，并且还为每次加入加载单独的坐标。
url = "http://membres-timc.imag.fr/Olivier.Francois/Arabidopsis/A_thaliana_chr1.geno"
download.file(url = url, destfile = "./A_thaliana_chr1.geno")
url = "http://membres-timc.imag.fr/Olivier.Francois/Arabidopsis/at_coord.coord"
download.file(url = url, destfile = "./at_coord.coord")
4.2.2使用snmf函数评估K值
注：此时我们认为这些物种是这些数据的单倍体。
obj.at = snmf("/Users/bcl/A_thaliana_chr1.geno", K = 1:10, ploidy = 1, entropy = T,
CPU = 1, project = "new")
plot(obj.at, col = "blue4", cex = 1.4, pch = 19)
4.2.3可视化K=5时候祖先系数矩阵
qmatrix = Q(obj.at, K =5)
因为只有个别数据，没有群体样本。因此我们不能像第一个例子一样使用图表映射。相反我们应该计算admixture系数的空间估计。我们在一个地理图上表示空间预测。为达到这一目的，需要一个代表欧洲的光栅网格。光栅可以从GIS应用程序下载或者网上。
asc.raster="http://membres-timc.imag.fr/Olivier.Francois/RasterMaps/Europe.asc"
grid=createGridFromAsciiRaster(asc.raster)
constraints=getConstraintsFromAsciiRaster(asc.raster, cell_value_min=0)
coord.at = read.table("at_coord.coord")
接下来，我们使用POPSutilities.R函数组中的maps函数来执行祖先系数的空间插值。输出一个漂亮颜色的admixture系数图。使用max选项，只有对祖先的最大贡献在地图的每个地理点上表示。
maps(matrix = qmatrix, coord.at, grid, constraints, method = "max",
main = "Ancestry coefficients", xlab = "Longitude", ylab = "Latitude", cex = .5)
map(add = T, interior = F)
