R语言计算BLUP值
时间: 2024-01-25 11:10:13 浏览: 274
***UP值的示例代码:
```r
# 安装和加载blup包
install.packages("blup")
library(blup)
# 读取数据
data <- read.csv("data.csv")
# 构建基因型矩阵
G <- model.matrix(~ SNP1 + SNP2 + SNP3, data = data)
# 构建表型矩阵
Y <- as.matrix(data$Phenotype)
# 构建随机效应模型
model <- mixed.solve(Y = Y, X = G, Z = NULL, K = NULL, REML = TRUE)
# 提取加性效应的BLUP值
additive_blup <- model$u
# 提取显性效应的BLUP值
dominant_blup <- model$v
# 提取母体效应的BLUP值
maternal_blup <- model$w
# 提取BLUP值大于0.1的个体
selected_individuals <- data$ID[model$u > 0.1]
# 生成ID列和类型
ID <- data$ID
type <- data$Type
result <- data.frame(ID, type)
# 输出结果
print(additive_blup)
print(dominant_blup)
print(maternal_blup)
print(selected_individuals)
print(result)
```
相关问题
R语言计算blup值代码
在 R 语言中,计算 Blup 值需要使用 BLUP 程序包。以下是一个简单的 Blup 计算代码示例:
首先,安装 BLUP 程序包:
```
install.packages("BLUP")
```
然后,加载 BLUP 程序包:
```
library(BLUP)
```
接着,准备数据。假设我们有一个数据框 df,其中包括了目标动物的表现数据和家族信息:
```
df <- data.frame(
animal_id = c(1, 2, 3, 4, 5),
phenotype = c(120, 140, 110, 130, 150),
sire_id = c(NA, 1, 2, 1, 2),
dam_id = c(NA, 3, 4, 3, 4)
)
```
其中,animal_id 是动物的唯一标识符,phenotype 是表现数据,sire_id 和 dam_id 是亲代的 animal_id。
接下来,构建一个数据框,包括所有动物的 animal_id 和遗传联系矩阵(A 矩阵):
```
A <- makeAinv(df$animal_id, df$sire_id, df$dam_id)
A_df <- data.frame(animal_id = df$animal_id, A = A)
```
然后,使用 lmer 函数拟合一个线性模型:
```
model <- lmer(phenotype ~ (1 | animal_id), data = df)
```
最后,使用 blup 函数计算 Blup 值:
```
blup_values <- blup(model, A_df, "A")
```
运行完毕后,blup_values 将包含每个动物的 Blup 值。
注意:这只是一个简单的示例,实际应用中需要根据具体情况调整代码。
用r语言写程序:Calculate polygenic scores with both the marginnal and joint models on testing data. And compare their performance under low and high h^2
以下是一个用R语言写的计算多基因得分(polygenic scores)并分别采用边缘和联合模型在测试数据上评估其性能的示例程序。同时还会比较在低和高h^2情况下的表现,供您参考:
```R
# Load required packages
library(rrBLUP)
library(plyr)
# Load training and testing data
train_data <- read.csv("train_data.csv")
test_data <- read.csv("test_data.csv")
# Define the phenotype and genotypes
pheno <- train_data$pheno
geno <- train_data[,4:ncol(train_data)]
# Define the relationship matrix
A <- A.mat(geno, n.core=4)
# Fit the marginal and joint models
marginal_model <- kin.blup(pheno, A, method="average")
joint_model <- mixed.solve(pheno ~ A, random=1)
# Calculate the polygenic scores
marginal_scores <- as.numeric(marginal_model$pred[test_data[,1]])
joint_scores <- as.numeric(joint_model$u[test_data[,1]])
# Calculate the correlation between true and predicted phenotypes
marginal_corr <- cor(test_data$pheno, marginal_scores)
joint_corr <- cor(test_data$pheno, joint_scores)
# Compare the performance under low and high h^2
h2_range <- seq(0.1, 0.9, by=0.1)
results <- data.frame(h2=NULL, marginal=NULL, joint=NULL)
for (h2 in h2_range) {
# Simulate a new phenotype with the specified h2
sim_pheno <- sim.giv(h2, A=A, K=A, n=1)
# Calculate the polygenic scores
marginal_scores <- as.numeric(marginal_model$pred[sim_pheno$id])
joint_scores <- as.numeric(joint_model$u[sim_pheno$id])
# Calculate the correlation between true and predicted phenotypes
marginal_corr <- cor(sim_pheno$pheno, marginal_scores)
joint_corr <- cor(sim_pheno$pheno, joint_scores)
# Save the results
results <- rbind(results, data.frame(h2=h2, marginal=marginal_corr, joint=joint_corr))
}
# Plot the results
library(ggplot2)
ggplot(results, aes(x=h2)) +
geom_line(aes(y=marginal, color="Marginal")) +
geom_line(aes(y=joint, color="Joint")) +
scale_color_manual(values=c("blue", "red")) +
labs(x="h2", y="Correlation") +
theme_classic()
```
请注意:这只是一个示例程序,您需要根据您的具体数据和研究问题进行修改和调整。
阅读全文
相关推荐
















