问题的由来

参考文献《Multi-view BLUP: a promising solution for post-omics data integrative prediction》
表型预测是加速植物育种的一种很有前途的策略，随着大数据时代的来临，对于同一种性状的表征往往可以采用不同的数据进行推断。例如从基因组层面，转录组层面，表观组层面都可以对表型进行推断。主要问题：如何整理不同来源的数据并进行表型推断？
作者在这里提出了Multi-view BLUP的框架，用于整合不同来源的数据（基因组层面，转录组层面，表观组层面记作不同来源的数据，在计算机中称之为多视图）并推断其表型

理论部分

差分进化算法的目的是，利用不同层次数据计算亲缘关系矩阵，并全局最优化亲缘关系矩阵，此时该亲缘关系矩阵融合了不同层次的数据信息

差分进化算法的步骤分为：初始化，变异，交叉，选择

1. 初始化

设解空间内存在NP个个体（即种群大小为NP），每个个体是 d 维向量，初始种群随机产生，可依靠任意函数生成，上述为一个例子

不同行代表不同种群，行向量（X₁~X₄₀代表种群个体的特征向量

2. 变异

差分进化算法使用种群中两个不同的特征向量来干扰一个现有向量，进行差分操作，来实现变异；g 代表种群代数

从当前群体中随机选择 3 个互不相同的种群个体 r₁，r₂，r₃；c 代表变异因子（缩放因子），g 代表种群代数；由此得到变异种群个体特征向量 V₁~V₄₀

3. 交叉

对于种群中每个个体和它所生成的子代变异向量进行交叉，即对每一个特征向量按照一定的概率选择子代变异向量（否则就是原种群个体的特征向量）来生成试验个体

CR（即 δ）为交叉概率因子；randb(j) 为随机自然数，可理解为阈值；大于该阈值保持原种群个体的特征向量，赋予对应索引的 X₁~X₄₀ 向量；小于该阈值则选择子代变异向量 V₁~V₄₀ 向量

4. 选择

通过以上的变异、交叉和选择操作，种群进化到下一代并反复循环，直到算法迭代次数达到预定最大次数，或种群最优解达到预定误差精度时算法结束

判断最优解需要设置目标函数，该例子的目标函数 f(X) 是混合线性回归模型；按照一定比例分train和test部分，利用train部分训练混合线性回归模型，用test部分得到预测值，计算test部分预测值与真实值的相关性，在一定迭代数中，满足前后两代的40个种群个体的平均准确率向量小于阈值时，停止计算，此时得到的亲缘关系矩阵即为融合了8个层次数据的亲缘关系矩阵

代码部分

整个流程：

source('./MVBLUP/R/GP.R')
source('./MVBLUP/R/MVBLUP.R')
source('./MVBLUP/R/plot.R')
source('./MVBLUP/R/similarity_matrix.R')

library(rrBLUP)

NP <- 40 
n <- 8
CR <- 0.5
Mu <- 0.5
s0 <- 0
s1 <- 1
thre <- 0.0001
iter <- 50
cv <- 5
trid <- 1
output_path <- "./MVBLUP/res/"

# read the multi-view dataset
Multi_view1 <- read.table("./MVBLUP/data/Multi_view1.txt", header=TRUE, sep="", stringsAsFactors=0, check.names=FALSE)
Multi_view2 <- read.table("./MVBLUP/data/Multi_view2.txt", header=TRUE, sep="", stringsAsFactors=0, check.names=FALSE)
Multi_view3 <- read.table("./MVBLUP/data/Multi_view3.txt", header=TRUE, sep="", stringsAsFactors=0, check.names=FALSE)
Multi_view4 <- read.table("./MVBLUP/data/Multi_view4.txt", header=TRUE, sep="", stringsAsFactors=0, check.names=FALSE)
Multi_view5 <- read.table("./MVBLUP/data/Multi_view5.txt", header=TRUE, sep="", stringsAsFactors=0, check.names=FALSE)
Multi_view6 <- read.table("./MVBLUP/data/Multi_view6.txt", header=TRUE, sep="", stringsAsFactors=0, check.names=FALSE)
Multi_view7 <- read.table("./MVBLUP/data/Multi_view7.txt", header=TRUE, sep="", stringsAsFactors=0, check.names=FALSE)
Multi_view8 <- read.table("./MVBLUP/data/Multi_view8.txt", header=TRUE, sep="", stringsAsFactors=0, check.names=FALSE)
Multi_dataset <- list(Multi_view1, Multi_view2, Multi_view3, Multi_view4, Multi_view5, Multi_view6, Multi_view7, Multi_view8)

# read phenotype
Phenotype <- read.table("./MVBLUP/data/Phenotype.txt", header = TRUE, sep = "", stringsAsFactors = 0, check.names = FALSE)
train_test_id <- read.table("./MVBLUP/data/train_test_id.txt", header = TRUE, sep = "", stringsAsFactors = 0, check.names = FALSE)

# obtain single-view inner product for calculating multi-view kinship matrix
A_N <- similarity_inner_product(Multi_dataset, n)

# optimize the weights of various single-view datasets by employing the MVBLUP model
MVBLUP_results <- MVBLUP_weights(NP, n, CR, Mu, s0, s1, 
                                 thre, iter, cv, trid, A_N,
                                 Phenotype, train_test_id)


# multi-view kinship matrix obtained from the MVBLUP model
kinship_test <- MVBLUP_results$kinship_test
# testop1_ex: true phenotypes of the test set
testop1_ex <- MVBLUP_results$testop1_ex
# trainop1_ex: true phenotypes of the train set
trainop1_ex <- MVBLUP_results$trainop1_ex

# We output a demo result of MVBLUP model, which includes the optimal weight values for different single-view datasets and the number of iterations required for convergence.
write.table(MVBLUP_results$MVBLUP_information, paste0(output_path, "trait", trid, "-MVBLUP_information.txt"), row.names = FALSE)

# test the MVBLUP model and store the results
MVBLUP_accuracy <- MVBLUP_PreTest(kinship_test, testop1_ex, trainop1_ex)
write.table(MVBLUP_accuracy, paste0(output_path, "trait", trid, "-MVBLUP_accuracy.txt"), row.names = TRUE)

1. 关于数据的输入

Multi_view1.txt - Multi_view8.txt，相当于不同层面的数据：

Multi_view1.txt

Multi_view8.txt

行代表不同的sample，列表示不同的特征，可以是SNP信息，可以是基因表达谱数据。Multi_view1.txt为 100（sample）× 3000（特征数量）的矩阵，Multi_view8.txt为 100（sample）× 1000（特征数量）的矩阵。各个来源数据的样本数量相同，特征的数量（和值）可以不同。

Phenotype.txt，代表表型数据值

Phenotype.txt

train_test_id.txt，划分sample的训练集和测试集

train_test_id.txt

2. 关于函数 similarity_inner_product()

similarity_inner_product() 解析：

#=============Computing the inner product of a single-view dataset==============
similarity_inner_product <- function(Multis,n) {
  A_N <- array(0, dim = c(nrow(Phenotype), nrow(Phenotype), n))
  for (n0 in 1:n) {
    # Multis 中含有 8 个不同层次来源的数据，Multis[[n0]][, -1] 相当于将每一个层次的数据矩阵取出来并删除第一列sample名
    # as.numeric(x) - mean(as.numeric(x))) / sd(as.numeric(x) 相当于对同一层次数据矩阵中的每一个特征（每一列）的向量进行标准化
    data_S <- apply(Multis[[n0]][, -1], 2, function(x) {
      (as.numeric(x) - mean(as.numeric(x))) / sd(as.numeric(x))
    })
   # 计算 sample 特征向量之间的乘积以便后续计算 kinship 矩阵（协方差） ，A_M 为100x100的矩阵（sample 一共 100 个）
    A_M <- data_S %*% t(data_S)
   # 对角线改为 1
    A_M[diag(A_M) == 0] <- 1
    A_N[, , n0] <- A_M
  }
  return(A_N)
}

# example，计算不同层次来源的数据（矩阵）的 similarity
A_N <- similarity_inner_product(Multi_dataset, n)

Multi_dataset 代表不同层次来源的数据（矩阵，Multi_view1-Multi_view8矩阵的 list），n 代表有 8 个不同层次来源的数据

补充说明：

其中 x_i 代表第 i 个sample的特征向量；x_j 代表第 j 个sample的特征向量，上式计算的是计算 sample 之间的 similarity（协方差）

A_N <- similarity_inner_product(Multi_dataset, n)中 A_N 代表 8 个不同层次来源数据中，每一个不同层次来源数据（矩阵中）各个 sample 之间的 similarity

3. 关于函数 MVBLUP_weights()

MVBLUP_weights() 函数

MVBLUP_weights <- function(NP, n, CR, Mu, s0, s1, thre, iter, cv, trid, A_N, Phenotype, train_test_id) {
  # 设置 5 x 交叉验证的 sample
  num_ex <- floor(nrow(Phenotype) / cv)
  # 区分训练集和测试集，并将他们的 index 取出来
  index_sample <- match(train_test_id[, 1], Phenotype[, 1])
  # 根据 id 将表型值取出来赋予 dt_ex 
  dt_ex <- Phenotype[index_sample, -1]
  samplenames <- train_test_id[, 1]
  # trid 代表需要选择的性状，trid = 1 代表选择 Phenotype 矩阵中的第一个性状
  # Phenotype 矩阵中的第一个性状值赋予 dt_in
  dt_in <- dt_ex[, trid]
  # 如果 dt_in 矩阵中的元素有NA值，则取第一个性状值的均值覆盖该 NA 值
  dt_in[is.na(dt_in)] <- mean(dt_in, na.rm = TRUE)
  # 区分 5 x 交叉验证的训练集和测试集
  trainop1_ex <- dt_in[-c(1:num_ex)]
  testop1_ex <- dt_in[1:num_ex]
  
  # 将各个来源数据 sample 之间的 similarity 矩阵 A_N（list）赋予 AA
  AA <- array(0, dim = c(nrow(train_test_id), nrow(train_test_id), n))
  for (n0 in 1:n) {
    AA[,, n0] <- A_N[,, n0][index_sample, index_sample]
  }
  
  # 创建两个空的 dataframe，一个讨论 Training Accuracy 另一个讨论 Training Weight
  Training_Accuracy <- data.frame(iter = rep(1:iter, each = NP), group = rep(1:NP, iter), Accuracy = NA)
  Training_Weight <- data.frame(iter = rep(1:iter, each = n), Parameter = rep(1:n, iter), Weight = NA)
  num_test <- floor(length(trainop1_ex) / cv)
  
  # 初始化部分 F0 代群体
  # 随机生成个 8 x 40 的矩阵（0-1）之间，40 相当于种群群体
  BM <- matrix(runif(n * NP, 0, 1), nrow = n)
  F0 <- NULL
  # 利用差分进化算法表征 sample 间的亲缘关系
  ## 初始化 8 个不同来源（层次）的 sample 间的矩阵信息，初始化 kinship 矩阵
  ## 1. 初始化 40 个群体
  for (k in 1:NP) {
    b <- BM[1:n, k]
    # 将 8 个不同来源（层次）的 sample 间的矩阵分别乘上一个微扰动值 b[n0]^2，然后全部累加起来（Reduce）
    # 此处可理解为融合了 8 个不同层次数据制作亲缘关系矩阵的过程，见补充说明 8
    # AA[,,n0] 为每一层次数据的 similarity 矩阵
    A_M <- Reduce("+", lapply(1:n, function(n0) (b[n0]^2) * AA[,,n0]))
    rownames(A_M) <- samplenames
    colnames(A_M) <- samplenames
    ## 初始化 8 个不同来源（层次）的 sample 间的矩阵信息，得到 kinship 矩阵
    A <- similarity_kinship(A_M)
    A_train <- A[-c(1:num_ex), -c(1:num_ex)]
    # 用混合线性模型 training
    # 详见以下补充说明 2
    res_M1 <- Pre_ER5(trainop1_ex, A_train, num_test)
    # F0[k] <- res_M1[length(res_M1)] 长度为 1；F0 （1 × 40）为初始化时 40 个种群的平均准确率，见补充说明 7
    F0[k] <- res_M1[length(res_M1)]
  }
  # 提取每一次测试的准确率，统称为F0代测试的准确率，见补充说明 4
  Training_Accuracy[1:NP, ncol(Training_Accuracy)] <- F0
  # 提取准确率最高时，对应的种群矩阵，作为训练权重，统称为F0训练权重，见补充说明 5
  Training_Weight[1:n, ncol(Training_Weight)] <- BM[, which.max(F0)]
  F_N0 <- matrix(0, nrow = NP, ncol = iter)
  F_N0[, 1] <- F0
  
  # 2. Mutation 部分
  # 选出最优的种群权重值后，进行变异部分计算
  # iter 为迭代次数，可以理解为代数，从 2 开始是因为初始化部分相当于第一代
  for (j in 2:iter) {
    # 初始化 V 和 U 矩阵（8 × 40），40 相当于种群群体，参考 BM 矩阵
    V <- matrix(0, nrow = n, ncol = NP)
    U <- matrix(0, nrow = n, ncol = NP)
    # 初始化 U 矩阵，即 BM矩阵 8 x 40 的矩阵（0-1）之间，40 相当于种群群体
    U <- BM
    # 变异算法：
    # 遍历 40 个群体，见补充说明 1
    for (l in 1:NP) {
      # 产生1~40的自然数
      Index_NP <- c(1:NP)
      # 1~40中随机选三个数
      r <- sample(Index_NP[-l])[1:3]
      # 执行上述遗传算法 Mu = 0.5
      V[, l] <- BM[, r[1]] + Mu * (BM[, r[2]] - BM[, r[3]])
    }

    # 3. Crossover 部分
    # 阈值从标准正态分布中选取其中一个值
    randrow <- runif(n)
    # CR = 0.5 代表 δ 为交叉概率，见补充说明 2
    # V 和 U 矩阵为 8 × 40 的矩阵
    # 如果满足 randrow 最小值大于 CR 
    # 之前的 U 矩阵有 U <- BM
    if (min(randrow) > CR) {
      # 随机换一个变异的特征向量
      rid <- sample(nrow(BM))[1]
      U[rid, ] <- V[rid, ]
    # 如果不满足，选择小于 CR 的元素对应的向量索引，利用该索引将变异矩阵 V 的行向量赋予对应的 U 矩阵，其余部分保持 BM 的原有向量
    } else {
      U[which(randrow <= CR), ] <- V[which(randrow <= CR), ]
     }
    
    # 4. Selection 部分
    for (k in 1:NP) {
      # 将超过 1 或者不足 0 的 U 矩阵元素分别替换为 1，0
      U[U[, k] <= s0, k] <- s0
      U[U[, k] > s1, k] <- s1
      if (max(U[, k]) == 0) {
        U[, k] <- BM[, k]
      }
      # 类似于初始化部分，利用混合线性模型进行推理，然后计算测试集预测值与真实值的相关性，
      b <- U[1:n, k]
      A_M <- Reduce("+", lapply(1:n, function(n0) (b[n0]^2) * AA[,,n0]))
      rownames(A_M) <- samplenames
      colnames(A_M) <- samplenames
      A <- similarity_kinship(A_M)
      A_train <- A[-c(1:num_ex), -c(1:num_ex)]
      res_M1 <- Pre_ER5(trainop1_ex, A_train, num_test)
      # j 为迭代次数（可以理解为代数），第 j 代，k 为第 k 个种群
      # F_N0[k,j] 代表第 j 代的第 k 个个体的平均准确率，见补充说明 3
      F_N0[k, j] <- res_M1[length(res_M1)]
      if (F_N0[k, j] < F0[k]) {
        U[, k] <- BM[, k]
        F_N0[k, j] <- F0[k]
      }
    }
    BM <- U
    # 当满足前后两代的40个种群个体的平均准确率向量小于阈值时，停止计算，见补充说明 3
    if (max(F_N0[, j] - F_N0[, (j - 1)]) < thre) {
      F0 <- F_N0[, j]
      Training_Accuracy[((j - 1) * NP + 1):(j * NP), ncol(Training_Accuracy)] <- F0
      Training_Weight[((j - 1) * n + 1):(j * n), ncol(Training_Weight)] <- U[, which.max(F0)]
      break
    }
    F0 <- F_N0[, j]
    Training_Accuracy[((j - 1) * NP + 1):(j * NP), ncol(Training_Accuracy)] <- F0
    Training_Weight[((j - 1) * n + 1):(j * n), ncol(Training_Weight)] <- U[, which.max(F0)]
  }
  MVBLUP_information <- data.frame(parameters = c(paste0("View", 1:n, "_W"), "iter", "Trait"),
                                   values = c(BM[, which.max(F0)], j, colnames(dt_ex)[trid]))
  
  results <- list(Training_Accuracy = Training_Accuracy, 
                  Training_Weight = Training_Weight, 
                  MVBLUP_information = MVBLUP_information, 
                  kinship_test = A,
                  trainop1_ex = trainop1_ex,
                  testop1_ex = testop1_ex)
  return(results)
}

函数 Pre_ER5() 和 PR() 的解析：

# 定义用于测试的样本个数为 num_test =16
Pre_ER5<-function(dt,A_train,num_test){
  A <-A_train
  # 随机选择测试集，多次测试
  # 利用性状数据与初始化的亲缘关系矩阵 A_train 建立混合线性回归关系
  res1<-PR(dt,A,1:num_test)
  res2<-PR(dt,A,(num_test+1):(num_test*2))
  res3<-PR(dt,A,(num_test*2+1):(num_test*3))
  res4<-PR(dt,A,(num_test*3+1):(num_test*4))
  res5<-PR(dt,A,(num_test*4+1):nrow(A))
  res_Pcor<-c(res1$pcor_te,res2$pcor_te,res3$pcor_te,res4$pcor_te,res5$pcor_te)
  res_Pl<-c(res1$pl_te,res2$pl_te,res3$pl_te,res4$pl_te,res5$pl_te)
  # 计算测试集部分对应的预测数据与 pl_te 的相关性的均值，作为 training 准确性的指标，相关性越高，准确性越强；见补充说明 4
  Pcor5<-sum(res_Pcor)/5
  res=c(res_Pl,res_Pcor,Pcor5)
  return(res)
}

nut = num_test
PR<-function(dt,A,nut){
  dt_M<-dt
  # 将选作为测试集的数据进行 NA 化，见补充说明 3
  dt_M[nut]<-NA
  data <- data.frame(y=dt_M,gid=rownames(A))
  # 利用性状数据 dt_M 与初始化的亲缘关系矩阵 A_train 建立混合线性回归关系
  ans <- kin.blup(data=data,geno="gid",pheno="y",K=A)
  # 取出测试集部分对应的预测数据（未被NA），计算去除测试集部分后对应预测数据的均值
  # 上述两部分相加得到 pl_te；见补充说明 6
  pl_te<-ans$g[nut]+mean(dt[-nut]) # ans$g[nut] 为测试集部分对应的预测数据（未被NA）；mean(dt[-nut]) 为去除测试集部分后对应预测数据的均值
  # 计算测试集部分对应的预测数据与 pl_te 的相关性
  pcor_te<-cor(as.numeric(dt[nut]),pl_te)
  Model=list(pcor_te=pcor_te,pl_te=pl_te) 
  return(Model)
}

？为什么要加平均：pl_te<-ans$g[nut]+mean(dt[-nut])
初始化部分：

补充说明：

NP 为差分进化算法的群体数量，例子中定为40，通过模拟种群中的遗传变异来筛选出优质种群

res_M1 <- Pre_ER5(trainop1_ex, A_train, num_test) 是利用线性混合模型进行训练并且

dt_M 的数据展示为：

用作测试的部分替换为NA，利用性状数据 dt_M 与初始化的亲缘关系矩阵 A_train 建立混合线性回归关系 A <-A_train; ans <- kin.blup(data=data,geno="gid",pheno="y",K=A)。A_train 的数据展示为：

然后通过混合线性回归模型可以预测测试集 NA 的部分，如下：

原来NA部分被预测出来

对应NA部分现已被预测出数值

随机选择五次测试集做训练，计算每次训练测试集部分对应的预测数据（未被NA）与 pl_te 的相关性：pcor_te<-cor(as.numeric(dt[nut]),pl_te)；然后取五次训练的均值：Pcor5<-sum(res_Pcor)/5 作为训练准确性的指标，称为F0代测试的准确率

初始的种群群体矩阵 BM：

BM 为 8 x 40 的矩阵，NP=40相当于40次采样，8对应8个不同来源（层次）数据的权重；当训练准确性达到最大时，BM矩阵对应的特征向量为最佳的群体权重向量，作为F0训练权重，F0 训练权重向量（长度为40），每个元素代表每个 NP（1~40）对应的 Pcor5 值（详见Pre_ER5()函数）

pl_te 的解释：pl_te<-ans$g[nut]+mean(dt[-nut]) ans$g[nut] 为测试集部分对应的预测数据（未被NA），mean(dt[-nut]) 为去除测试集部分后对应预测数据的均值；取出测试集部分对应的预测数据（未被NA），计算去除测试集部分后对应预测数据的均值；上述两部分相加得到 pl_te

F0[k] <- res_M1[length(res_M1)] 长度为 1；F0 （1 × 40）为初始化时 40 个种群的平均准确率

此处可理解为融合了 8 个不同层次数据制作亲缘关系矩阵的过程：
# 此处可理解为融合了 8 个不同层次数据制作亲缘关系矩阵的过程
# AA[,,n0] 为每一层次数据的 similarity 矩阵
A_M <- Reduce("+", lapply(1:n, function(n0) (b[n0]^2) * AA[,,n0]))
rownames(A_M) <- samplenames
colnames(A_M) <- samplenames
## 初始化 8 个不同来源（层次）的 sample 间的矩阵信息，得到 kinship 矩阵
A <- similarity_kinship(A_M)

变异部分，Crossover 部分和 Selection 部分：

补充说明：

变异算法：

X^g_r1（向量长度为8）为BM[, r[1]]；X^g_r2（向量长度为8）为BM[, r[2]]；X^g_r3（向量长度为8）为BM[, r[3]]；Mu 为 c scaling factor。V[, l] <- BM[, r[1]] + Mu * (BM[, r[2]] - BM[, r[3]]) 相当于

差分进化算法使用种群中两个不同的特征向量来干扰一个现有向量，进行差分操作，来实现变异
# 遍历 40 个群体
for (l in 1:NP) {
   # 产生1~40的自然数
   Index_NP <- c(1:NP)
   # 1~40中随机选三个数
   r <- sample(Index_NP[-l])[1:3]
   # 执行上述遗传算法 Mu = 0.5
   V[, l] <- BM[, r[1]] + Mu * (BM[, r[2]] - BM[, r[3]])
}
Crossover 部分
# 从标准正态分布中选取 8 个值，赋予向量 randrow 
# CR = 0.5 代表 δ 为交叉概率
randrow <- runif(n)
# 如果满足 randrow 最小值大于 CR 
# 之前的 U 矩阵有 U <- BM
if (min(randrow) > CR) {
   # ？随机换一个变异的特征向量
   rid <- sample(nrow(BM))[1]
   U[rid, ] <- V[rid, ]
# 如果不满足，选择小于 CR（交叉概率） 的元素对应的向量索引，利用该索引将变异矩阵 V 的行向量赋予对应的 U 矩阵
} else {
   U[which(randrow <= CR), ] <- V[which(randrow <= CR), ]
}
对于种群中每个个体和它所生成的子代变异向量进行交叉，即对每一个特征向量按照一定的概率选择子代变异向量（否则就保持原种群个体的特征向量）来生成试验个体

Selection 部分：F_N0[k,j] （40 × iter，这里的 iter 代表代数）代表第 j 代的第 k 个个体的平均准确率；当满足前后两代的40个种群个体的平均准确率向量小于阈值时，停止计算

最优化后，此时推导出的亲缘关系矩阵即为融合了8个层次数据的亲缘关系矩阵，里面提取了不同层次的亲缘关系信息

Multi-view BLUP 原理介绍