Tag Snp：在给定的 SNP 中，连锁不平衡值（LD）最高的 SNP;

1. 先下载测试数据

https://github.com/YuetongXU/CropGBM-Tutorial-data

2.安装cropgbm

pip install --user cropgbm

3. 开始使用cropgbm

3.1基因型数据预处理，统计并展示基因型数据的缺失率和杂合率等

#conda activate cropgbm
cropgbm -o gbm_result -pg all --fileprefix genofile --fileformat ped --snpmaxmiss 0.10 --samplemaxmiss 0.10 --maf-max 0.05 --r2-cutoff 0.7

输出的pdf文件内容如下：

样本缺失率
SNP缺失率
genofile_maf.pdf

MAF（次要等位基因频率）的分布图，可以看到基因型的大概分布。

基因型预处理的主要内容

• 根据 缺失率、MAF 筛除不符合要求的样本个体及 snp 位点
• 根据高频基因型填补 snp 的缺失数据
• 根据 r^2 去除冗余 snp
• 基因型数据重编码（基因型 -> 012）
• 根据 snp 编号或样本编号提取特定样本的特定基因型数据
• 基因型数据的杂合率、缺失率和次等位基因频率分布情况的可视化

3.2 表型数据预处理

cropgbm -o gbm_result -pp --phe-plot --phefile-path phefile.txt --phefile-sep ',' --phe-name DTT --ppgroupfile-path phefile.txt --ppgroupfile-sep ',' --ppgroupid-name paternal_line --phe-norm

表型数据需要包含2列，第一列为样本编号。
输出的结果文件中的pdf是对表型的分布的箱线图可视化

phefile_scatter.pdf

图中的每个点代表一个样本，横坐标是分组，纵坐标是表型值。

表型预处理的主要内容

• 提取指定 样本ID、snpID 的表型数据
• 表型归一化 --phe-norm参数
• 表型重编码等。
  表型文件phefile.txt的内容如下：

seqname,paternal_line,pop_class,DTT
CUBIC_72_X_Tester_25,Tester_25,Reid,0.5408204825240001
CUBIC_74_X_Tester_25,Tester_25,Reid,0.176112246932
CUBIC_75_X_Tester_25,Tester_25,Reid,2.0793823313099997

参数说明：
-o 输出文件夹
-pp 调用表型预处理模块
--phe-plot 绘制箱线图
--phefile-path 表型文件的路径
--phefile-seq 分割的字符串
--phe-name 表型所在列的名称
--ppgroupfile-path 分组文件路径， 示例：此处是和表型文件在一起
--ppgroupid-name 分组变量的列名
--phe-norm 对表型进行归一化，使用的是z-score

3.3 群体结构分析

基于基因型数据分析样本的群体结构
使用t-SNE或PCA方法对基因型数据进行降维，
使用OPTICS或Kmeans方法降维
以散点图形式展示样本的群体结构。

cropgbm -o gbm_result/ -s --genofile-path genofile_filter.geno --structure-plot --redim-mode pca --pca-explained-var 0.98 --cluster-mode kmeans --n-clusters 30

参数说明：
-s 调用群体结构分析模块
--genofile-path 基因型文件，是预处理后的结果
--structure-plot 绘图
--redim-mode 降维算法，二选一：tsne或pca
--pca-explainded-var 参数值大于1的整数或0-1之间的小数。整数表示pca降维后数据的维度，参数值小数表示降维后数据可解释的变量值。该参数只在使用pca降维模式时有用
--cluster-mode 聚类算法 二选一：kmeans或optics
--n-clusters 形成的簇数以及生成的核心数，只在kmeans模式下有用。

输出的文件有：

• filename.redim 降维结果的输出文件
• filename.cluster 聚类结果的输出文件
• filename_redim.pdf 散点图展示降维结果

• filename_cluster.pdf散点图展示聚类结果

  
  redim.pdf

reachability.pdf genofile_filter_cluster.pdf

聚类的图

3.4构建模型与特征选择

建模之前先进行交叉验证
cropgbm -o gbm_result/ -e -cv --traingeno train.geno --trainphe train.phe --cv-nfold 5 --min-detal 0.5
-e 调用snp选择与表型预测模型
-cv 调用交叉验证功能
--traingeno 训练集基因型文件路径
--trainphe 训练集样本表型数据文件路径
--cv-nfold 交叉验证的折叠数，默认值是5
--min-detal 每次训练对模型精确度提升的最小百分值。默认是0.05

输入：基因型数据、表型数据
输出：filename.lgb_model

cropgbm -o gbm_result -e -t --traingeno train.geno --trainphe train.phe --validgeno valid.geno --validphe valid.phe --num-threads 48

3.5 SNP选择

cropgbm -o gbm_result -e -t -sf --bygain-boxplot --traingeno train.geno --trainphe train.phe --min-gain 0.05 --max-colorbar 0.6 --cv-times 6 --num-threads 48

-sf 调用snp选择功能，与-t训练模型连用
--bygain_boxplot 绘制随snp加入模型误差不断变化的箱线图
--min-gain 每个snp在模型中的总增益
--max-colorbar 每个snp在各树中的最大增益值
--cv-times 交叉验证的重复次数

输出文件：
train_bygain.pdf 散点图形式展示随着snp的加入，模型的预测误差的变化情况。x轴是每次模型新加入的snp的编号。
train_random.pdf 以箱线图形式展示随snp加入模型误差变化情况。x轴是随机抽选的snp
train_heatmap.pdf 热图形式展示各snp在不同决策树中的增益值及变化规律

train_bygain.pdf

横坐标是使用的增益SNP,可以看出，随着snp数量的增加，MSE的值逐渐降低到接近0.4之后，就趋于平缓。再增加snp的数量，也不会降低误差。这个可以帮助我们决定所需的最低的snp数量。

train_random.pdf
这是随机抽取的snp对模型的误差的贡献率，看出也是类似的曲线，到一定数量后，误差就不会继续降低了。但是看此时的稳定的误差的值是0.6左右。和前面一个图对比，说明增益SNP用于预测的准确度更高。
train_heatmap.pdf
热图形式展示各个snp在不同决策树中的增益值及变化规律，x轴为snp编号，按照featureGain值从大到小排列，y轴为决策树编号。tree20表示第20次训练所生成的决策树。显然拍在前面的snp的增益值较高，对模型的贡献率较大，结合前面的MSE的降低情况，前面的这些snp肯定就是与模型训练的表型紧密相关的位点。

3.6 表型预测

cropgbm -o gbm_result -e -p --testgeno test.geno --modelfile-path train.lgb_model --num-threads 48

-p 调用表型预测功能
--testgeno 测试集的基因型文件，该文件第一行是snp id,第1列是样本编号信息。
--modelfile-path 前面生成的模型的路径
输入文件是基因型文件，要求 分隔符是逗号，而且基因型已经经过预处理转为0,1,2这种数字矩阵
输出文件是filename.predict 为预测结果的输出文件

我使用了自己的数据测试该程序，一堆报错，程序似乎没人维护了。