1.R包的下載(R version 3.6)
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("LEA")
如果R版本是老版本的刻获，請查看網(wǎng)站Bioconductor release的要求進行R包下載。
注：具體怎么在R中進行操作在以下網(wǎng)站中有詳細說明炕婶。
http://membres-timc.imag.fr/Olivier.Francois/tutoRstructure.pdf
2.加載其他函數(shù)
source("http://membres-timc.imag.fr/Olivier.Francois/Conversion.R")
source("http://membres-timc.imag.fr/Olivier.Francois/POPSutilities.R")
3導入輸入文件
#struct2geno函數(shù)的講解
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:意味著每個個體使用兩行數(shù)據(jù)對標記進行編碼柠掂。
FORMAT=1:意味著每個個體使用一行數(shù)據(jù)對標記進行編碼。
extra.row/.col:表示額外行/列數(shù)依沮。
"./genotype.geno"：輸出文件的路徑涯贞。
4.使用R包進行群體結構分析
4.1example1：Allelic markers(等位基因標記)
分析在歐洲二次contract后通過空間顯式結合模擬所產生的群體遺傳數(shù)據(jù)。在最初的一個分離的階段后危喉，一個物種開始從遙遠的南方殖民到歐洲宋渔，一些在伊利比亞群島，一些在土耳其辜限。第二次contract發(fā)生在中世紀的歐洲皇拣，靠近德國。數(shù)據(jù)由10個二倍體個體的60個族群樣本組成列粪，用100個等位基因標記對10個二倍體個體進行基因分型审磁。
4.1.1數(shù)據(jù)的下載和轉換
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函數(shù)設置完成谈飒。
library(LEA)
obj.snmf = snmf("secondary_contact.geno", K = 3, alpha = 100, project = "new")
qmatrix = Q(obj.snmf, K = 3)
4.1.3 600個個體祖先系數(shù)柱狀圖 {Q-matrix(用barplot呈現(xiàn))}
barplot(t(qmatrix), col = c("orange","violet","lightgreen"), border = NA, space = 0,
xlab = "Individuals", ylab = "Admixture coefficients")
4.1.4 通過在歐洲地理地圖上疊加餅圖法進行admixture群體估計
為達到這個目標岂座，我們需要讀取樣本的地理坐標和創(chuàng)建樣本標志符。
共有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 選擇集群數(shù)
在LEA里選擇集群數(shù)是基于交叉熵準則，這一準則也被admixture程序所使用手素。交叉熵準則基于一部分被掩蓋基因型的預測和交叉驗證的方法鸳址。更小的交叉熵準則的值通常意味著能夠更好的運行。我們執(zhí)行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數(shù)據(jù).這些數(shù)據(jù)是從Horton等人所研究的大量數(shù)據(jù)中抽出來的一部分。調查了一號染色體上1096個近交生態(tài)型(52001個SNPs)的群體結構崩哩。
4.2.1SNP數(shù)據(jù)的下載
數(shù)據(jù)保存在磁盤中巡球，并且還為每次加入加載單獨的坐標。
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函數(shù)評估K值
注：此時我們認為這些物種是這些數(shù)據(jù)的單倍體邓嘹。
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時候祖先系數(shù)矩陣
qmatrix = Q(obj.at, K =5)
因為只有個別數(shù)據(jù)酣栈，沒有群體樣本。因此我們不能像第一個例子一樣使用圖表映射汹押。相反我們應該計算admixture系數(shù)的空間估計矿筝。我們在一個地理圖上表示空間預測。為達到這一目的棚贾，需要一個代表歐洲的光柵網(wǎng)格窖维。光柵可以從GIS應用程序下載或者網(wǎng)上榆综。
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函數(shù)組中的maps函數(shù)來執(zhí)行祖先系數(shù)的空間插值铸史。輸出一個漂亮顏色的admixture系數(shù)圖奖年。使用max選項，只有對祖先的最大貢獻在地圖的每個地理點上表示沛贪。
maps(matrix = qmatrix, coord.at, grid, constraints, method = "max",
main = "Ancestry coefficients", xlab = "Longitude", ylab = "Latitude", cex = .5)
map(add = T, interior = F)
