机器学习应用于基因组预测，以苜蓿为例

原创

生物信息与育种 2024-01-18 15:30:27 博主文章分类：GS ©著作权

文章标签 数据机器学习预处理 文章分类 Html/CSS 前端开发

©著作权归作者所有：来自51CTO博客作者生物信息与育种的原创作品，请联系作者获取转载授权，否则将追究法律责任

目的
示例代码

表型文件预处理
基因型文件预处理

机器学习建模预测
绘图
其他结果

利用代码复现一个机器学习应用于基因组预测的项目，张志武老师于2023年发表在《园艺研究》上的一篇文章。

机器学习应用于基因组预测，以苜蓿为例_预处理

目的

使用 GWAS 和 GP 结合重测序数据和从世界各地收集的 220 份紫花苜蓿种质的秋季休眠（Fall dormancy, FD），定义为秋剪后 30 天的植物再生高度：（i）比较不同GP方法的预测准确性;（ii）使用GWAS和GP的混合模型计算FD预测精度;（iii）扩展GP模型以进行新的性状预测。

示例代码

表型文件预处理

2017年温室种植，2018年移植大田，按三个重复的完全随机区组试验设计。2018年割茬三次，2019年割茬四次。秋季休眠是最后一次割茬后1个月的植物再生高度。使用广义线性模型（GLM）估计秋季休眠的方差：

机器学习应用于基因组预测，以苜蓿为例_数据_02

广义遗传力计算：

机器学习应用于基因组预测，以苜蓿为例_机器学习_03

表型文件整理格式如下：

机器学习应用于基因组预测，以苜蓿为例_预处理_04

基因型文件预处理

#过滤保留双等位基因
vcftools --vcf mdp_genotype.vcf --min-alleles 2 --max-alleles 2 --maf 0.05 --recode --recode-INFO-all --out snp_maf0.05_allele2.vcf 
#通过LD过滤
plink --vcf snp_maf0.05_allele2.vcf --indep-pairwise 100 50 0.2 --out snp_LD_filter.vcf --allow-extra-chr --make-bed 
#挑选LD保存的标记，生成过滤后标记
plink --vcf snp_LD_filter.vcf --extract snp_filter.vcf.prune.in --out prunData --recode 12 --allow-extra-chr 
#转换回vcf文件
perl /mnt/e/GWAS_model_test/tassel_v5/run_pipeline.pl -fork1 -plink -ped prunData.ped -map prunData.map -export snp_use_LD.vcf -exportType VCF -runfork1 
#利用beagle软件填充缺失值
java -Xmx5g -jar beagle.22Jul22.46e.jar nthreads=20 gt=snp_use_LD.vcf out=snp_use_impute_out 
#解压
gzip -d snp_use_impute_out.vcf.gz 
#改成数值型格式，把0|0改成0,0|1改成1,1|1改成2，直接sed也行
perl change_vcf_format_impute.pl snp_use_impute_out.vcf > snp_numeric_format.vcf 
# #去掉不需要的列
cut -f3,10- snp_numeric_format.vcf > R_use_format.vcf 
#去掉注释信息，生成模型构建需要用的文件
grep -v '^#' R_use_format.vcf >R_use_format1.vcf

预处理后基因型文件格式如下：

机器学习应用于基因组预测，以苜蓿为例_机器学习_05

机器学习建模预测

使用Python 包 sklearn 用于构建两类模型，一是SVM回归，包括线性核SVM linear和多核SVM poly；二是三种正则化模型，包括Lasso（L1正则化）、ridge regression（L2正则化）和 ElasticNet 模型（L1和L2正则化结合）。

选择使用的5种模型分别代表了控制表型的不同QTN模型，Lasso代表少量大效应的QTN模型，Ridge代表大量微效QTN, ElasticNet代表Lasso和Ridge的中间类型。SVM（support-vector machines）是一种有监督的机器学习模型，从多纬度层面评估不同QTN间关系，分别从线性和非线性层面模型构建。

交叉验证方法：80%个体建模，20%个体作为预测，重复100次。代码如下：

#加载数据处理模块
import numpy as np
import pandas as pd
import scipy as sp
import statsmodels.api as sm 

#加载机器学习模块
from sklearn import datasets, linear_model 
from sklearn.linear_model import Ridge
from sklearn.kernel_ridge import KernelRidge
from sklearn.linear_model import ElasticNet
from sklearn.cross_decomposition import PLSRegression
from sklearn.metrics import mean_squared_error, r2_score 
from sklearn.svm import SVR
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import ElasticNetCV 

#导入基因型数据
data = pd.read_csv('R_use_format1.vcf',sep="\t",index_col=0)
#行列转置
data1=data.T 
#导入表型数据
phe = pd.read_csv('pheno_all_ind_mean.txt',sep="\t",index_col=0) 
#提取个体名称
index=phe.index 
#匹配表型和基因型数据的个体信息
data1 = data1.reindex(index = index) 
#挑选需要的表型
use_phe=phe["dormancy_ind_mean_all"] 
#去掉表型为NA的个体
use_phe=use_phe[use_phe.notna()] 
#把基因型数据也去掉对应的个体
GD1=data1.loc[use_phe.index] 
#定义一个导出数据的文件
data_results={} 
for i in range(100):
        train_y = use_phe.sample(frac =.8) #拆分出80%个体为建模群体
        train_x = GD1.loc[train_y.index] #得到对应的80%基因型数据
        test_y = use_phe.drop(train_y.index,axis=0) #剩下20%个体作为验证群体
        test_x = GD1.drop(train_y.index,axis=0) #20%个体对应的基因型数据
        ##以下模型可更换，去掉注释，选择其一
        regr = Ridge(alpha=0.1) #岭回归模型
        #regr = linear_model.Lasso(alpha=0.1) #lasso模型
        #regr = ElasticNet(random_state=0) #ElasticNet模型
        #regr = KernelRidge(alpha=1.0) #KernelRidge模型
        #regr = PLSRegression(n_components=2) #PLSRegression模型
        #regr = ElasticNetCV(cv=5, random_state=0) #ElasticNetCV模型
        #regr = SVR(kernel="linear", C=100, gamma="auto") #SVR-linear模型
        #regr = SVR(kernel="poly", C=100, gamma="auto") #SVR-poly模型
        #regr = linear_model.LinearRegression() #线性模型

        regr.fit(train_x, train_y.values.ravel()) #构建模型
        y_pred = regr.predict(test_x) #预测表型
        y_pred = pd.DataFrame(y_pred) #修改表型格式
        test_all=test_y.reset_index() #修改验证表型的个体顺序
        test_all['Pred_y']=y_pred #合并预测表型和验证表型
        test_all.columns = ['taxa', 'Real_Y', 'Predicted_Y'] #添加表头
        corr = test_all['Real_Y'].corr(test_all['Predicted_Y']) ##calculate correlation info ##预测数据和真实数据相关性
        #coef = sm.OLS(test_all['Predicted_Y'], test_all['Real_Y']).fit().params[0] ##回归系数
        data_results[i]=corr #保存结果
#保存数据
test_all=pd.DataFrame.from_dict(data_results, orient='index')
#导出csv格式
test_all.to_csv('alfalfa_100rep_ind_mean_SVR_poly_coefficient.csv', sep=',')

得到结果文件如下：

机器学习应用于基因组预测，以苜蓿为例_预处理_06

绘图

根据结果绘图。

library(ggplot2)
mydata<-read.table("combine_five_method_results.txt",header=TRUE,sep="\t")
mydata$Accuracy=abs(mydata$Accuracy) #相关系数有正负，准确率只考虑正值

p1<-ggplot(mydata, aes(x=Phenotype, y=Accuracy, fill=Method))+labs(x="Phenotype", y="Prediction accuracy") +
  geom_boxplot(width=0.8)+ylim(0,1)+
  theme_classic()+theme(
    legend.position="right",legend.title = element_text(color = "black", size = 10), legend.text = element_text(color = "black", size = 10),axis.text = element_text(color="black",size=10),axis.title.x=element_text(margin = margin(t = 10, r = 0, b = 0, l = 0),size=10),axis.title.y=element_text(size=10))+
  scale_fill_manual(values=c("#DD3226","#5EBCD5","#F88301","#07B294","#F1AC00"))

ggsave(plot=p1,"GS_accuracy.pdf", width = 18, height = 8, units = "cm")

机器学习应用于基因组预测，以苜蓿为例_预处理_07