育种数据分析之放飞自我分享 http://blog.sciencenet.cn/u/yijiaobai 关注:生物统计,数量遗传,混合线性模型,生物信息,R,Perl,Python,GWAS,GS相关方法,文章及代码

博文

可以做structure的R语言包:LEA

已有 1875 次阅读 2019-5-29 20:33 |个人分类:数量遗传学|系统分类:科研笔记| structure, LEA

关于分群的软件,之前写了structure 2.3.4 软件使用指南,软件虽然有windows版本,但是操作太麻烦了,也写了Admixture使用说明文档cookbook,但是只有Linux版本,使用起来有难度。难道不能使用R语言进行structure绘图么?结果来了:LEA!

1. paper

LEA: An R package for landscape and ecological association studies

使用说明文档

不同格式的数据使用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.

划重点:

  • 可以在R语言中实现软件Structure的功能

  • 可以做类似admixture的图

  • 简单操作, 几个命令实现相关功能

  • C语言开发, 可以处理大数据

3. 软件安装

install.packages(c("fields","RColorBrewer","mapplots"))
source("http://bioconductor.org/biocLite.R")
biocLite("LEA")

如果安装不成功, 也可以通过CRAN把软件包下载到本地, 进行安装:

install.packages("LEA_1.4.0_tar.gz", repos = NULL, type ="source")

载入两个函数, 进行格式转化以及可视化:

source("http://membres-timc.imag.fr/Olivier.Francois/Conversion.R")
source("http://membres-timc.imag.fr/Olivier.Francois/POPSutilities.R")

4. 测试数据

plink格式的ped文件, 具体格式参考:plink格式的ped和map文件及转化为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
个体ID
父本
母本
性别
表型值
SNP1-1(SNP1的第一个位点)
SNP1-2(SNP的第二个位点)

测试数据采用admixture的示例数据, 使用plink将其转化为ped文件

library(LEA)
# 结果会生成test.geno文件的数据.
output = ped2lfmm("test.ped")
# 使用LEA进行structure进行分析
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")

在这里插入图片描述

对比admixture的结果

# 对比admixture结果
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选择最优K值

# 绘制折线图, 选择最优K值.
plot(project, col = "blue", pch = 19, cex = 1.2)

在这里插入图片描述
可以看出, K=3时, 最小, 因此选择K=3.

r-breeding.png



http://blog.sciencenet.cn/blog-2577109-1181899.html

上一篇:我72岁开始学习python,花了一个星期--参加R语言会议有感
下一篇:通过基因组选择预测杂交水稻的表现(数据挖掘)

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2021-1-23 06:46

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部