前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >blupf90 VS Plink VS GCTA 基因型数据PCA分析

blupf90 VS Plink VS GCTA 基因型数据PCA分析

作者头像
邓飞
发布2019-06-13 20:21:35
1.3K0
发布2019-06-13 20:21:35
举报
文章被收录于专栏:育种数据分析之放飞自我

这篇是软件测试. 主要是用同一批数据, 测试不同软件和方法对结果的影响, 不同方法有:

  • BLUPF90构建G矩阵, 计算PCA
  • BLUPF90构建H矩阵, 计算PCA
  • PLINK构建G矩阵, 计算PCA
  • GCTA构建G矩阵, 计算PCA

结论: BLUPF90构建G, H, Plink构建G矩阵, 结果一致. GCTA画的PCA有出入, 怀疑是默认参数的问题.

模拟一套数据, 5个世代, 最后三代有基因型数据, 每个世代400个个体, SNP为50K.

1. blupf90构建G矩阵的PCA

blupf90如果想要进行GBLUP分析, 不写系谱信息即可, 示例par文件:

代码语言:javascript
复制
DATAFILE
dat_f90.txt
TRAITS
10         # This is column 10 (phenotype) from QMSim data file
FIELDS_PASSED TO OUTPUT
1          # This will copy the ID number to the renf90.dat data file
WEIGHT(S)  # WARNING: ALWAYS PUT AN EMPTY LINE AFTER THIS!!!!!

RESIDUAL_VARIANCE
0.9                # add starting values for residual variance
EFFECT
4 cross alpha    # Fit generation as a fixed effect, 'cross alpha' is a class in SAS
EFFECT
1 cross alpha    # Fit animal effect
RANDOM
animal           # Fit animal effect (A matrix) for the effect directly above it (column 1, animal)
#FILE
#ped_f90.txt  #  pedigree file (animal, sire, dam), 0's are missing always!!!
#FILE_POS
#1 2 3 0 0        # indicates that column 1 = Animal, column 2 = Sire, column 3 = Dam
SNP_FILE
yM.txt
(CO)VARIANCES    
0.1               # add starting values for additive animal effect
OPTION alpha_size 25               # Equal to the max number of characters within a column
OPTION max_string_readline 800  # maximum number of characters in one line of data file
OPTION max_field_readline 100   # maximum number of columns in the dataset
#OPTION saveHinvOrig
#OPTION saveHinv
#OPTION sol se
#OPTION use_yams
#OPTION missing -999
OPTION plotpca

运行preGSf90后, 会生成pc1vspc2文件, 里面包括PC1和PC2两列, 增加世代为pop, 然后使用R画图:

代码语言:javascript
复制
pca = read.table(("pc1vspc2"))
head(pca)
names(pca) = c("PC1","PC2")
pca$pop = rep(c("A","B","C"),each=400)
library(ggplot2)
p <- ggplot(pca, aes(x=PC1, y=PC2, colour=pop)) 
p <- p + geom_point(size=2)
p <- p + stat_ellipse(level = 0.95, size = 1)
p <- p + geom_hline(yintercept = 0) 
p <- p + geom_vline(xintercept = 0) 
p <- p + theme_bw()
p

结果:

2. blupf90构建H矩阵的PCA

需要定义系谱和基因型, 示例par文件:

代码语言:javascript
复制
DATAFILE
dat_f90.txt
TRAITS
10         # This is column 10 (phenotype) from QMSim data file
FIELDS_PASSED TO OUTPUT
1          # This will copy the ID number to the renf90.dat data file
WEIGHT(S)  # WARNING: ALWAYS PUT AN EMPTY LINE AFTER THIS!!!!!

RESIDUAL_VARIANCE
0.9                # add starting values for residual variance
EFFECT
4 cross alpha    # Fit generation as a fixed effect, 'cross alpha' is a class in SAS
EFFECT
1 cross alpha    # Fit animal effect
RANDOM
animal           # Fit animal effect (A matrix) for the effect directly above it (column 1, animal)
FILE
ped_f90.txt  #  pedigree file (animal, sire, dam), 0's are missing always!!!
FILE_POS
1 2 3 0 0        # indicates that column 1 = Animal, column 2 = Sire, column 3 = Dam
SNP_FILE
yM.txt
(CO)VARIANCES    
0.1               # add starting values for additive animal effect
OPTION alpha_size 25               # Equal to the max number of characters within a column
OPTION max_string_readline 800  # maximum number of characters in one line of data file
OPTION max_field_readline 100   # maximum number of columns in the dataset
#OPTION saveHinvOrig
#OPTION saveHinv
#OPTION sol se
#OPTION use_yams
#OPTION missing -999
OPTION plotpca

运行preGSf90后, 会生成pc1vspc2文件, 里面包括PC1和PC2两列, 增加世代为pop, 然后使用R画图:

代码语言:javascript
复制
pca = read.table(("pc1vspc2"))
head(pca)
names(pca) = c("PC1","PC2")
pca$pop = rep(c("A","B","C"),each=400)
library(ggplot2)
p <- ggplot(pca, aes(x=PC1, y=PC2, colour=pop)) 
p <- p + geom_point(size=2)
p <- p + stat_ellipse(level = 0.95, size = 1)
p <- p + geom_hline(yintercept = 0) 
p <- p + geom_vline(xintercept = 0) 
p <- p + theme_bw()
p
3. plink根据G矩阵做PCA

代码:

代码语言:javascript
复制
plink --file b --pca 3

结果生成:

代码语言:javascript
复制
plink.eigenval  plink.eigenvec  plink.log  plink.nosex

R语言作图:

代码语言:javascript
复制
library(ggplot2)
head(dd)
p <- ggplot(dd, aes(x=PC1, y=PC2, colour=pop)) 
p <- p + geom_point(size=2)
p <- p + stat_ellipse(level = 0.95, size = 1)
# p <- p + scale_color_manual(values = cols)
p <- p + geom_hline(yintercept = 0) 
p <- p + geom_vline(xintercept = 0) 
p <- p + theme_bw()
p
4. gcta64根据G矩阵做PCA

将ped文件转化为bed文件

代码语言:javascript
复制
plink --file b --make-bed --out c

生成grm文件

代码语言:javascript
复制
gcta64 --bfile c --autosome --make-grm --out grm

生成pca文件

代码语言:javascript
复制
gcta64 --grm grm --pca 3

根据PCA信息作图

结论

blupf90的G矩阵, H矩阵, plink的PCA结果一致. GCTA构建的PCA结果不太一致, 怀疑是参数默认的有问题, 回头查看一下.

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2019-04-28,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 育种数据分析之放飞自我 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1. blupf90构建G矩阵的PCA
  • 2. blupf90构建H矩阵的PCA
  • 3. plink根据G矩阵做PCA
  • 4. gcta64根据G矩阵做PCA
  • 结论
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档