关注我们：
【GS专栏】11-全基因组选择实战之RF和SVM
知识/政策规划
2021-04-12
1. 前言
前面我们已经演练了基因组选择常用的R包，包括对BLUP和Bayes这两类模型的实战。今天我们进入机器学习应用于GS的实战环节：随机森林(RandomForest，RF)和支持向量机(Support Vector Machine，SVM)。理论部分我们已经提到，RF和SVM除了可以解决常见的分类问题，还可以处理回归问题。因此，对于大量数量性状的预测，机器学习也就有了用武之地。
2. 回归预测
大部分农艺性状都属于数量性状，也就是说表型数值是连续的。在机器学习中将对这种数据建模预测称之回归。因此，对于基因组选择而言，我们在育种应用中需要处理的大部分其实也是回归问题。本次演练的示例数据仍然来自CIMMYT的小麦数据GYSS[1]。为节省运行时间，以其中242个小麦样本的1000个标记为例。
# 加载数据
data(GYSS)
# 基因型数据
dim(Markers)
Markers[1:5, 1:5]
# 表型数据
head(phenotype)
浏览表型和基因型数据。
> dim(Markers)
[1] 242 1000
> Markers[1:5, 1:5]
PZB00451.1 PZB00416.2 PZB00579.1 PZB00257.1 PZB00236.2
DT10 0 0 0 0 2
DT100 0 2 0 0 0
DT101 0 0 2 0 0
DT102 0 0 0 0 0
DT103 2 0 0 0 0
> head(phenotype)
[1] 0.7045494 0.3886150 0.1165665 0.3747609 0.4136186 0.6908267
2.1 RFR
随机森林的回归，即Random Forest Regression，简称RFR。在R语言中，随机森林算法常用randomForest包来实现。这个R包使用起来非常简单，同名主函数虽然参数众多，但需要设置的并不多，大部分默认即可，除非你很清楚要那么做。
# 加载包
library(randomForest)
# 建模
rf_fit <- randomForest(
x = Markers, #编码后的基因型
y = phenotype, #表型
ntree = 500, #树的数目
importance = F, #是否计算各个变量在模型中的重要性
nodesize = NULL #树的节点大小
)
# 预测
rf_pred <- predict(rf_fit, Markers)
# 预测能力
rfr <- cor(rf_pred, phenotype)
注：这里为了简化说明，并没有划分数据集。模型的验证也是采用的全部预测值和表型计算相关性。下面同理。
2.2 SVR
支持向量机的回归，即Support Vector Regression，简称SVR。我们采用R包e1071来分析。
# 加载包
library(e1071)
# 建模
svm_fit <- svm(
x = Markers, #编码后的基因型
y = phenotype, #表型值
kernel = "linear", #核函数
cost = 2^(-9), #惩罚因子，可与任意核函数搭配
gamma = 1, #数据映射到新特征空间后的分布，值越大，支持向量越少
probability = TRUE #模型是否允许概率预测
)
# 预测
svm_pred <- predict(svm_fit, Markers)
# 预测能力
svr <- cor(svm_pred, phenotype)
RFR和SVR适用于数量性状预测，如以上的株高性状。
3. 分类预测
对于质量性状或分类变量，我们自然地会想到用RF和SVM来解决。那么能否将株高这样的数量性状也当成分类问题来处理呢?思路上是可以这么处理的：将株高性状排序，将之归类为高、中、低三组，然后视为分类问题来解决。首先，我们可以编写一个用于分类的函数，如按比例来进行样本分类。
sampleClassify <- function(phenotype,
posPercentage = 0.4,
BestIndividuals = c("top", "middle", "bottom")) {
trainSampleSize <- length(phenotype)
posSampleSize <- round(posPercentage*trainSampleSize)
if (BestIndividuals == "top") {
trainPhenOrder <- order(phenotype, decreasing = TRUE)
} else if (BestIndividuals == "bottom") {
trainPhenOrder <- order(phenotype, decreasing = FALSE)
} else if (BestIndividuals == "middle") {
trainPhenOrder <- order(abs(phenotype), decreasing = FALSE)
}
posIdx <- trainPhenOrder[1:posSampleSize]
negIdx <- setdiff(c(1:trainSampleSize), posIdx)
res <- list(posSampleIndex = posIdx, negSampleIndex = negIdx)
return res
}
然后，将表型和基因型分别分类。如按前40%进行划分，高的个体的表型值归为“1”组，低的个体归为“0”组。
posNegSampleList <- sampleClassify(
phenotype = phenotype,
posPercentage = 0.4,
BestIndividuals = "top"
)
markers <- Markers[c(posNegSampleList$posSampleIndex,
posNegSampleList$negSampleIndex), ]
pheVal <-as.factor(c(rep("1",length(posNegSampleList$posSampleIndex), rep("0", length(posNegSampleList$negSampleIndex)
)))
数据划分好后，就可以用随机森林和支持向量机的分类器来进行处理了。
3.1 RFC
随机森林分类，即Random Forest Classifier，简称RFC。同样用R包randomForest来实现，用法是一样的。
# 建模
rfc_fit <- randomForest(
x = markers, #分类后的标记
y = pheVal, #分类后的表型
ntree = 500,
importance = F,
nodesize = NULL
)
# 预测
rfc_pred <- predict(rfc_fit, markers, type = "vote")[,"1"]
# 预测能力
rfc <- cor(rfc_pred, phenotype)
3.2 SVC
支持向量机的分类，即Support Vector Classifier，简称SVC。同样用R包e1071实现，用法是一样的。
# 建模
svc_fit <- svm(
x = markers, #分类后的标记
y = pheVal, #分类后的表型
kernel = c("linear"),
cost = 2^(-9),
gamma = 1,
probability = TRUE
)
# 预测
svc_pred <- predict(svc_fit, Markers, probability = TRUE)
svc_res <- as.matrix(attr(svc_pred, "probabilities")[, "1"])
# 预测能力
svc <- cor(svc_res, phenotype)
我这里仅仅作为演示，就不展示具体的结果了。各位感兴趣的话，可以用自己的数据测试一下，看看分类和回归的效果如何。
4. 后记
至此，我们基本上已经演练了一遍GS中常见的三类模型——BLUP、Bayes和机器学习。首先，恭喜各位老师同学顺利通过学习!但需要说明的是，本期及往期实战只是简单演示了GS数据分析中的核心部分，而非完整的流程，分析过程也不如真实数据那么严谨。在实际的数据分析过程中，应包括前期数据的预处理和统计分析、数据抽样划分、特征选取、模型选择与比较、交叉验证、模型性能评估、预测个体选择等环节。每一步都非常重要，涉及的细节也很多，需要大家利用实际数据多加练习应用。后续若有机会，这些内容也将逐步分享给大家。若有任何疑问，欢迎垂询交流。本期内容就到这里啦，谢谢大家!
相关推荐
百奥繁育 | 良种选育的信息管理系统
09-06
2023
8月25日，江西省水稻种质资源创新交流及智能育种技术研讨会在江西省农科院顺利召开。江西农业大学、江西师范大学、广东省农科院、江西省农科院、中国水稻所江西早稻中心、南昌市农科院、赣州市农科所有关专家代表以及省内种企代表共约40人参与了本次交流会议。
08-28
2023
有人说，植物育种的过程就像工厂流水线。挑选的种质即原材料，根据市场需求进行设计与加工，送入田间生产线后，一代代优中选优，还得经过严格的产量、品质、抗性测试，过五关斩六将，最终拿到审定编号，成为可以推广的成熟产品。 虽然流程相似，但育种工作远不及工厂生产那样标准。田间环境气候复杂、真实性状判断困难、水肥条件难以统一……种种困难下，选育良种成了概率事件，育种家们都有种“尽人事听天命”的无奈感。
08-22
2023
育种4.0时代，比拼的是科技创新，关键是通过数字化、信息化、智能化让育种过程缩短周期、提升效率、降低风险，而这也是长沙百奥云数据科技有限公司（百奥云）的创业初衷。
08-22
2023
1. 利用PopGenomics快速对基因组实现高质量组装 基因组组装是将全基因组测序的小片段（reads长度100 bp-100 kb）通过算法拼接成尽量长的片段（contig 和scaffold，长度几十kb 到Mb 不等）或者整条染色体的过程。获得包含基因组全序列的参考基因组是对动植物进行基因组学研究和育种利用的前提[1]。 由于植物基因组具有非常丰富的多样性，参考已发表的少数物种组装新的物种，有时却无法达到理想的组装效果。测序技术发展提供了短序列测序、单分子测序、光学图谱、Hi-C 图谱等多种测序技术及其组合的组装方案，到目前为止，已经有上千个植物基因组被组装（图1）。然而，如何以最低成本快速获得满足需求的高质量基因组，仍然是科研人员和育种数据分析人员普遍面临的一个问题。
08-22
2023
7月26日，顶着炎炎夏日 在公司CSO桂進老师的带领下 举行了一场别开生面的劳动体验。
08-04
2023