|
關(guān)于分群的軟件，之前寫了structure 2.3.4 軟件使用指南，軟件雖然有windows版本，但是操作太麻煩了，也寫了Admixture使用說明文檔cookbook，但是只有Linux版本，使用起來有難度。難道不能使用R語言進(jìn)行structure繪圖么？結(jié)果來了：LEA！ 1. paperLEA: An R package for landscape and ecological association studies 使用說明文檔 不同格式的數(shù)據(jù)使用LEA 2. 軟件介紹This short tutorial explains how population structure analyses reproducing the results of the widely-used computer program structure can be performed using commands in the R language. The method works for any operating systems, and it does not require the installation
of structure or additional computer programs. The R program allows running population structure inference algorithms, choosing the number of clusters, and showing admixture coefficient bar-plots using a few commands. The methods used by R are fast and accurate, and they
are free of standard population genetic equilibrium hypotheses. In addition, these methods allow their users to play with a large panel of graphical functions for displaying pie-charts and interpolated admixture coefficients on geographic maps.
劃重點(diǎn): 可以在R語言中實(shí)現(xiàn)軟件Structure的功能 可以做類似admixture的圖 簡(jiǎn)單操作, 幾個(gè)命令實(shí)現(xiàn)相關(guān)功能 C語言開發(fā), 可以處理大數(shù)據(jù)
3. 軟件安裝install.packages(c("fields","RColorBrewer","mapplots"))
source("http:///biocLite.R")
biocLite("LEA")如果安裝不成功, 也可以通過CRAN把軟件包下載到本地, 進(jìn)行安裝: install.packages("LEA_1.4.0_tar.gz", repos = NULL, type ="source")載入兩個(gè)函數(shù), 進(jìn)行格式轉(zhuǎn)化以及可視化: source("http://membres-timc./Olivier.Francois/Conversion.R")
source("http://membres-timc./Olivier.Francois/POPSutilities.R")4. 測(cè)試數(shù)據(jù)plink格式的ped文件, 具體格式參考:plink格式的ped和map文件及轉(zhuǎn)化為012的方法 1 SAMPLE0 0 0 2 2 1 2 3 3 1 1 2 1
2 SAMPLE1 0 0 1 2 2 1 1 3 0 4 1 1
3 SAMPLE2 0 0 2 1 2 2 3 3 1 4 1 1前六列為:
家系ID
個(gè)體ID
父本
母本
性別
表型值
SNP1-1(SNP1的第一個(gè)位點(diǎn))
SNP1-2(SNP的第二個(gè)位點(diǎn)) 測(cè)試數(shù)據(jù)采用admixture的示例數(shù)據(jù), 使用plink將其轉(zhuǎn)化為ped文件 library(LEA)
# 結(jié)果會(huì)生成test.geno文件的數(shù)據(jù).
output = ped2lfmm("test.ped")# 使用LEA進(jìn)行structure進(jìn)行分析
library(LEA)
obj.snmf = snmf("test.geno", K = 3, alpha = 100, project = "new")
qmatrix = Q(obj.snmf, K = 3)
head(qmatrix)
barplot(t(qmatrix), col = rainbow(3), border = NA, space = 0,
xlab = "Individuals", ylab = "Admixture coefficients")
對(duì)比admixture的結(jié)果 # 對(duì)比admixture結(jié)果
qad = read.table("test.3.Q")
head(qad)
barplot(t(qad), col = rainbow(3), border = NA, space = 0,
xlab = "Individuals", ylab = "Admixture coefficients")
5. 使用snmf選擇最優(yōu)K值# 繪制折線圖, 選擇最優(yōu)K值.
plot(project, col = "blue", pch = 19, cex = 1.2)
可以看出, K=3時(shí), 最小, 因此選擇K=3.
|
