DNNGP中文使用文档

wheat1.tsv 制表符分割的表型文件，第一列是id,第2列是表型值
wheat599_pc95.tsv 制表符分割的主成分矩阵文件，第一列是id,后面是pc1到pc251列
wheat599_pc95.pkl是模型可读取的主成分矩阵文件

从vcf文件获取PCA矩阵

plink2 --threads 30 --vcf test.vcf --pca 200 --out pca200 --allow-extra-chr
cat pca200.eigenvec|sed 's/#IID/ID/' >pca200.tsv
#这里是使用前300行作为训练的样本
head -301 pca200.tsv >pca200.train.tsv

#最后的20行作为要预测的样本
cat <(head -1 pca200.tsv) <(tail -20 pca200.tsv) >pca200.predict.tsv
tsv2pkl.py pca200.train.tsv pca200.train.pkl
tsv2pkl.py pca200.predict.tsv pca200.predict.pkl

结果会生成两个文件，后缀名分别为.eigenval 和.eigenvec，eigenval 显示了每个
PC 所占的比重，各个 PC 的比重/比重和为每个 PC 的解释度。eigenval 为我们需
要使用的 PCA 矩阵

DNNGP解决的一个核心就是SNP文件的高维度问题，使用PCA把snp降低维度，示例中降低到95个维度（95个主成分PCA），然后使用这些作为基因型文件，进行预测。因为是使用PCA的维度作为训练的值，所以在开始训练之前，必须把训练集和要预测的基因型放到一起计算获得PCA,然后再分开成训练的tsv和预测的tsv,使用tsv2pkl.py脚本把这2个文件分别转为pkl文件。

原作者没有提供脚本tsv2pkl.py，我自己编写了一个tsv2pkl.py.可以实现tsv文件和pkl文件的互相转换。注意后缀名必须是对应的格式tsv或pkl,脚本是通过判断后缀名，完成转换格式的。
pkl转为tsv用法
python3 tsv2pkl.py wheat599_pc95.pkl wheat599_pc95.tsv
tsv转为pkl用法
python3 tsv2pkl.py wheat599_pc95.tsv wheat599_pc95.pkl

预处理步骤：

1.使用plink2把vcf文件转为PCA的主成分,获取tsv文件。注意：一定要把训练的和要预测的基因型合并到一个vcf里面，然后在获得PCA的tsv文件后，再分为训练数据的tsv和要预测的tsv
2.使用脚本tsv2pkl.py把PCA的TSV文件转为pkl文件。

训练数据集

mkdir Test_files
cd ~/AI/DNNGP/DNNGP-main/
python3 Scripts/dnngp_runner.py --batch_size 28 --lr 0.001 --epoch 100 --dropout1 0.5 --dropout2 0.3 --patience 5 --seed 123 --cv 10 --part 1 --earlystopping 10 --snp "Input_files/wheat599_pc95.pkl" --pheno "Input_files/wheat1.tsv" --output ./Test_files/

参数说明：
--batch_size 训练模型所调用的样本量
--lr 初始学习率
--epoch 迭代次数
--dropout1 第一次特征抛弃
--dropout2 第2次特征抛弃
--patience 无提升则减小学习率阈值
--seed 随机种子
--cv 交叉验证折数
--part 设定选取第几折数据作为验证集
--earlystopping 无提升则停止训练阈值，
--snp 训练集的基因型文件，格式是pkl
--pheno 表型文件，第一列是id,第2列是表型值
--output 输出文件夹

输出结果如下：
('Corr obs vs pred =', PearsonRResult(statistic=0.5732735932580789, pvalue=1.690124919853689e-06)) Pearson相关系数是0.57，pvalue值是1.69e-6.即通过基因型预测表型的相关性达到0.57.
('Prediction validation save in:', './Test_files/Prediction.validation.csv') 验证集的预测结果
('Model save in:', './Test_files/training.model.h5') 训练得到的模型文件
('Model history save in:', './Test_files/Modelhistory.csv') 训练的历史记录

• Prediction.validation.csv：DNNGP 模型对验证集的预测结果（第一列的序号代表预测值个体在原数据集中的名称）。
• training.model.h5：训练好的模型文件，用于下一步对育种群体表型性状预测。
• Modelhistory.csv：记录了训练过程中，各项数值的变化情况，打开后可以看到各项指标的变化情况，从这里可以选择决断，什么时间终止训练，随着训练层数的增加，各项指标不再继续变优就是决断层数
  自己的数据要根据实际情况调整参数，尽量获得更高的相关系数和更低的Pvalue值
  我试着通过降低--lr，可以一定程度提升模型的Pearson值。

使用模型预测

python3 ~/AI/DNNGP/DNNGP-main/Scripts/Pre_runner.py --Model "~/AI/DNNGP/DNNGP-main/Test_files/training.model.h5" --SNP ~/AI/DNNGP/DNNGP-main/Input_files/wheat599_pc95.pkl" --output "~/AI/DNNGP/DNNGP-main/Test_predict/

参数说明：
--Model 上一步训练获得的模型文件，格式是*.h5
--SNP 指定要预测的基因型的pkl文件。
--output 指定输出路径，注意最后面一定要有/
输出结果：
Prediction.ALL.csv
文件有2列，第1列是ID,第2列是预测的值。

该算法存在的问题，当你拥有某一个物种的新的群体的基因型时，想预测该群体的表型，你需要把你的基因型和已知表型的基因型合并到一起进行PCA,然后再分开，进行建模预测。因为预测使用的是降维之后的PCA,所以每次有新的基因型都得重复这么做一遍，这是不友好的，特别是如果你训练的数据集的样本数量比较多时，每次都需要重新训练模型。
理想的模型应该是，你只要是使用的同一个参考基因组比对获得的基因型，那就直接调用模型进行预测即可，即不需要在重新建模。但是这种是需要你拥有基因型覆盖足够全面，而且样本数量足够大，材料代表性足够强，这种群体的规模至少是要数千个材料才行，否则群体不具有代表性，预测出的准确性很低。
每一个算法都不一定是你的数据的最优算法模型，机器学习的结果的可重复性并不高，一定要根据实际情况，调整参数，选择最优的结果。可能论文里的准确性很高，你实际运行后，没有获得那么高的R2，调参后只要接近即可。