GWAS表型数据处理

1 读取数据

sbw <- read.table("E:/bio/R/projects/gwas/df.txt",header = T,sep = "\t")
head(sbw)

2 描述性统计分析

均值、方差、最大值、最小值、变异系数、偏度、峰度、正态性检验p值计算

2.1 自行计算

存在问题，最后结果为一个嵌套的list，不好导出。
我的解决方法data.frame(analy)使其输出到屏幕，再复制粘贴到excel。
优点在于所求数据就是需要的数据，不需要过多的整理

sbw_mean <- read.table("E:/bio/R/projects/gwas/df.txt",header = T,sep = "\t")
sum(is.na(sbw_mean))
stastic <- function(x){
  min <-min(x)
  max <-max(x)
  m <- mean(x)
  n <- length(x)
  s <- sd(x)
  cv <- s/m
  skew <- sum((x-m)^3/s^3)/n
  kurt <- sum((x-m)^4/s^4)/n - 3
  sht <- shapiro.test(x)
  return(c(Num = n, Min = min, Max = max, Mean = m, SD = s,CV = cv, Kurtosis = kurt,Skewness = skew,shapiro.test = sht[2]))
  }
analy <- sapply(sbw_mean,stastic)
data.frame(analy)
write.table(unlist(data.frame(analy)),"test.txt")#输出的结果也还可以，但是shapiro.test这行不好处理

2.2 `pastecs`包（推荐）

简单方便，已于导出，数据稍加处理即可使用。
stat.desc(x,basic=T,desc=T,norm=F,p=0.95)
其中 x 代表数据框，basic默认值为TRUE，计算所有值、空值、缺失值的数量、最大值、最小值、值域、总和。desc默认值为TRUE，计算中位数、均值、均值的标准误差、均值置信度为95%的置信区间、方差、标准差、变异系数。若norm=TRUE，返回偏度和峰度（以及他们的统计显著程度）、Shapiro-Wilk正态检验结果。

install.packages("pastecs")
library(pastecs)
analy2 <- stat.desc(sbw_mean,basic = TRUE,desc = T,norm = T)

3.计算相关系数及可视化

3.1 相关系数计算

install.packages("Hmisc")
library(Hmisc)
cor <- rcorr(as.matrix(sbw_mean))
cor

3.2相关系数可视化

install.packages("PerformanceAnalytics")
library(PerformanceAnalytics)
chart.Correlation(sbw_mean, histogram=TRUE, pch="19")

4 BLUP分析&遗传力计算

4.1 BLUP计算

Error:package or filenames load failed for ‘lme4’ in inDL(x, as.logical(local), as.logical(now), ...):
程序包'Rcpp_precious_remove'不提供'Rcpp'这样的函数
遇到报错，把R包update即可，应该是由于调用的R包版本过旧导致的

install.packages("lme4")
library(lme4)
bl_sbw <- read.table("E:/bio/R/projects/gwas/df.txt",header = T,sep = "\t")
for(i in 1:4) bl_sbw[,i] <- as.factor(bl_sbw[,i])
blp =  lmer(SBW ~ var + (1|loca) +(1|loca:rep) + (1|year),data=bl_sbw)
summaryblp <- summary(blp)
anovablp <- anova(blp)
ranblp <- ranef(blp)
fixble <- as.data.frame(fixef(blp))

来自的公式

name_varcmp = lmer(phenotype_name~ (1|LINE) + (1|LOC) + (1|YEAR) + (1|REP %in% LOC:YEAR) + (1|LINE:LOC) + (1|LINE:YEAR),
                   control=lmerControl(check.nlev.gtr.1 = "ignore"))

注意：植物中，一般的BLUE值需要加上截距（Intercept）。因为BLUE值中，第一个水平会当做0，其它为相对值，可以手动进行相加，也可以使用lsmeans包中的lsmeans。

plot(sbw_mb$SBW_mean,sbw_mb$SBW_BLUP)
hist(sbw_mb$SBW_BLUP,main = "BLUP",xlab = "BLUP" )

4.2 遗传力计算

summary(blp_tt)

image.png

公式如下：