校正表型
getwd()
setwd("/path/")
## Reading raw phenotypic set data: including animal ID, potential fixed effects
library(openxlsx)
phenotype <- read.table("Target_fat_deposition_traits.txt", header = T, row.names = 1)
dim(phenotype)
colnames(phenotype)[1:10]
phenotype
指定固定效应 注意确定固定效应不共线
fit_phenotype <- phenotype
covars <- c("Farm","season")
# 制造一个数据框用于存放校正后的表型
adjust_phenotype_df<-matrix(0,nrow(fit_phenotype),ncol(fit_phenotype))
rownames(adjust_phenotype_df)<-rownames(fit_phenotype)
colnames(adjust_phenotype_df)<-colnames(fit_phenotype)
adjust_phenotype_df[,1]<-fit_phenotype[,1]
dim(adjust_phenotype_df)
校正表型
for(i in 3:ncol(fit_phenotype)){
# 为了方便后面合并，制造用于合并的表格
stat <- matrix(0,nrow(phenotype),1)
colnames(stat) <- c("sampleID")
stat[,1] <- rownames(fit_phenotype)
pheno<-colnames(fit_phenotype)[i]
lm_formula <- paste(pheno, paste(covars, collapse = " + "), sep = " ~ ")
lm_fit <- lm(formula = as.formula(lm_formula), data = fit_phenotype)
pheno_adjust<-as.data.frame(lm_fit[["residuals"]])
colnames(pheno_adjust)<-pheno
pheno_adjust$sampleID<-rownames(pheno_adjust)
pheno_adjust_2 <- merge(stat,pheno_adjust, by.x = "sampleID", all.x = TRUE)
adjust_phenotype_df[,i]<-pheno_adjust_2[,2]
gc()
}
adjust_phenotype_df_2 <- adjust_phenotype_df[,3:ncol(adjust_phenotype_df)]
write.table(adjust_phenotype_df_2,"adjust_phenotype.txt",quote = F)
