具有遗传性疾病和性状的遗传位点分析

1、参赛密码 （由组委会填写）全“华为杯”第十三届全国研究生数学建模竞赛学 校江苏科技大学参赛队号队员姓名1.孙佳伟2.李袁3.李肇基参赛密码 （由组委会填写） “华为杯”第十三届全国研究生数学建模竞赛题 目 具有遗传性疾病和性状的遗传位点分析摘 要：本文设计了基于属性重要度的选择算法，并通过SVM分类器构建出预测模型，对不同的位点和基因进行分析，判断每个位点或基因对某种疾病的预测精度，从而判断是否为致病位点或者致病基因。最后利用该算法和模型，预测出十种性状相关的致病位点。问题1，针对每个位点有碱基对组成的性质，为了方面描述和分析，本文采用了十进制编码，每个位点的属性值可以通过0-9中的一个数进行

2、表示，具体的编码格式，文中给出了详细的编码表。问题 2，设计了基于属性重要度的特征选择算法，通过SVM分类器构建出预测模型，通过问题1中的特征表示方式，提取所有样本每列的特征，并对每列的特征属性进行重要度分析，进而判断该疾病与位点rs, rs,rs,rs,rs相关。 问题 3，每个基因是由不同位点组成的集合，则每个基因的所有特征属性即对应集合里位点特征属性的集合，利用问题2优化的模型，通过预测精度，对每个基因的对某疾病的重要度进一步分析得出，该疾病与致病基因gene_171相关。问题 4，利用本文提出的模型，对10种性状中的每种性状中，继续通过属性重要度分析，识别出在不同性状中最有可能的致病位

3、点，最后得出10个形状的相关致病位点分别为：rs, rs, rs, rs, rs, rs, rs, rs, rs, rs。本文亮点是，提出基于属性重要度的选择算法，通过SVM构建出预测模型，利用网格搜索进行寻优，判断每列属性的重要度，从而判断致病位点或者致病基因。 关键词：属性重要度；SVM分类器；优化模型；位点(SNPs)一 问题重述问题1，请用适当的方法，把genotype.dat中每个位点的碱基(A,T,C,G)编码方式转化成数值编码方式，便于进行数据分析。问题2，根据附录中1000个样本在某条有可能致病的染色体片段上的9445个位点1的编码信息(见genotype.dat)和样本患有遗

4、传疾病A的信息（见phenotype.txt文件）。设计或采用一个方法，找出某种疾病最有可能的一个或几个致病位点，并给出相关的理论依据。问题3，同上题中的样本患有遗传疾病A的信息（phenotype.txt文件）。现有300个基因，每个基因所包含的位点名称见文件夹gene_info中的300个dat文件，每个dat文件列出了对应基因所包含的位点(位点信息见文件genotype.dat)。由于可以把基因理解为若干个位点组成的集合，遗传疾病与基因的关联性可以由基因中包含的位点的全集或其子集合表现出来请找出与疾病最有可能相关的一个或几个基因，并说明理由。问题4，在问题二中，已知9445个位点，其编码

5、信息见genotype.dat文件。在实际的研究中，科研人员往往把相关的性状或疾病看成一个整体，然后来探寻与它们相关的位点或基因。试根据multi_phenos.txt文件给出的1000个样本的10个相关联性状的信息及其9445个位点的编码信息(见genotype.dat)，找出与multi_phenos.txt中10个性状有关联的位点。二 问题分析问题1的分析 问题1提供的1000个样本，每个样本均有9445个位点，每个位点是一个碱基对，本文通过十进制对所有碱基对2进行编码，利用此特征表示方法，每个样本得到9445个属性，每个属性通过0-9中的数值进行表示。问题2的分析问题2和问题3都是要利

6、用某种方法，进行致病位点和致病基因的检测，预测某种疾病的致病位点。预测该疾病相关的致病位点，其实就是判断不同位点对该疾病的影响程度，即判别每个位点的属性重要度。本文设计了基于属性重要度的选择算法，并利用SVM构建预测器，计算每个位点对疾病的预测精度，从而得到与该疾病相关的致病位点。问题3的分析问题3与问题2比较，问题2是识别与某疾病相关的致病位点，而某个基因是由多个位点组成的集合，问题3是识别某疾病相关的致病基因，为了识别致病基因，由于基因是由多个位点组成的集合，则基因的特征属性即多个位点特征属性组成的集合。利用每个基因的特征属性，通过问题2中设计的算法和构建的预测模型，从而识别与某疾病相关的

7、致病基因。问题4的分析 问题4中，人的某些疾病是和性状相关的，材料中提供了10种性状，要求判断与这10种性状的相关的致病位点。可以对每个性状分析，通过问题2种设计的选择算法和构建的预测模型，识别出某个性状相关的致病位点。三 模型假设和符号说明3.1 模型假设（1）假设给出的样本数据能分别代表整个的正常和患病群体。（2）忽略寻优时SVM分类器本身造成的偏差。（3）给出的位点和基因均是有效的位点和基因。3.2 符号说明（未说明）（1） A,C,G,T ：DNA中的四种碱基（2） Acc ：即accuracy，预测精度（3） PCn ：第n个位点的属性（4） S ：特征属性子集（5） ：通过十进制编

8、码得到位点子集的特征属性（6） ：第i个位点的属性，在第1组属性选择中进行融合（7） Sig ：即significant，属性重要度四、模型建立与求解4.1问题1：每个位点碱基编码方式转化4.1.1数据分析本文所使用位点测试数据集，来自1000个可能致病的染色体片段试验检测结果，标签分布为500个无病染色体使用0表示，500个患病染色体使用1表示，且每个致病染色体上有9445个碱基对，以此作为位点。4.1.2十进制编码本文采用十进制0,1,2,.9编码将每个碱基对转化成数据编码方式，以便于数据分析。“AA”为“0”；“AC”为“1”；“AG”为2；“AT”为3.“TT”为9，详见碱基对编码表1

9、（其中AC,CA；CG,GC；.碱基对表示方式相同）。表1碱基对编码ACGTA0136C1247G3458T6789另外，位点中出现字符I和D，根据说明，分别用T和C代替4.2 问题2：找出与疾病最有可能相关的一个或几个位点4.2.1 位点属性矩阵由于所有样本序列上的本一个二核苷酸位点代表了一个属性，本文总共有9445个位点即9445个不同的属性，这些属性由十进制表示（见附图1）。其中，属性列中PC1PCn 表示9445个不同的属性指标；AA,AC,AG,AT,.,TT表示16中不同的原始二核苷酸。图1 十进制编码编码碱基对4.2.2 实验测试方法和分类器设计1. 实验测试方法K折交叉验证、J

10、ackknife（留一法）测试和独立数据集测试是三种常用的实验测试方法。在本试验中采取了K折交叉验证，K折交叉验证是指经过初始采样后将初始样本分割成K个子样本，然后选择一个单独的子样本作为测试数据，剩下的K-1个样本用来训练模型。将这个交叉验证重复K次，使得每个子样本都验证一次，将K次结果取平均，得到一个估测结果。这个方法优势在于随机产生的子样本的训练和验证能够同时进行，每次的结果验证一次，本文采用了5折交叉验证。2. 分类器设计特征向量提取完成后便要面临对向量的分类问题，选取一个合适的分类器，将大大提高相关识别问题的精度，然而分类器研究至今出现了许多，例如贝叶斯分类器3，支持向量机（SVM）

11、4，5，6，7，K近邻分类器8，9,10以及DeepLearning11，12。在本篇论文中，我们采用支持向量机对样本进行分类。支持向量机是由Vapnik 等根据统计学习理论提出的，在考虑结构风险最小化原则和VC维理论的基础上，在有限的样本信息中，寻找模型复杂性和学习能力的最佳折中点。这个分类方法在解决模式分类与非线性映射问题中非常适用。例如图2便是在二维平面中一个分类，线AB上的点组成的向量就是这个样本的支持向量，然而在多数情况下数据样本将不仅仅局限于二维，通常都是在高维空间中的样本，如图3，这时的支持向量是以面的形式表现出来。图2 支持向量机二维原理图图3 支持向量机高维原理图在生物信息学

12、中，非线性以及高维数据的处理十分常见，所以支持向量机在生物信息学的领域被广泛使用。因此本篇文章中，libsvm-3.18包作为支持向量机的实现。由于在非线性分类过程的有效性和速度，因为径向基核函数（RBF）利用网格搜索法优化了正则化参数和核参数，所以我们使用这个函数对样本进行预测。从支持向量机获得的概率得分被用来预测本次实验的最终结果。4.2.3 评价标准为了对我们的预测方法进行合理的评价,我们在考虑结果的准确性上，加入了相关系数，用这些数据对预测结果进行一系列评价。本文使用下列度量模型的性能进行了评价：准确度（Acc），它可以表示为 (1)其中，TP表示被正确判断为正样本的正样本数；FP表示

13、被错误判断为正样本的负样本数；TN表示被正确判断为负样本的负样本数；FN表示被错误判断为负样本的正样本数。将公式（1）化简，并定义新的字符表示，化为如下形式： (2)用代替TP，代替FN，代替TN，代替FP，显然，表示所有的正样本数。同理，表示所有的负样本数。根据此式，我们可以很容易的得出以下的结构。当=0时，即所有的样本都被正确判断时，；当且时，即所有的样本都被错误判断时，；Acc是对整个数据集的评价指标。所以我们在比较预测器性能的时候，主要对比Acc的值。4.2.4 基于属性重要度的选择算法设计1. 属性重要度问题1中通过十进制编码方法，对样本中9445个位点进行特征表示，得到9445个特

14、征属性，且通过实验发现，多个属性的组合来进行该模型的预测可能比单个属性预测的结果要更高，但是并不是所有的属性都与该疾病有必然的联系，可能存在冗余的属性，即可能存在冗余的位点，本文设计一种基于属性重要的选择算法，通过属性选择的方式可以有效的去除冗余属性，提高判断的准确性，从而保障了预测的可靠性。2. 属性选择算法假设第n列的属性记为PCn,则全部9445个位点的属性集合为PC1，PC2.PCnPC9445，选取位点属性子集为S，通过问题1中的编码方式编码样本序列特征记为。 （3）公式(3)中的和分别表示，采用位点属性子集和性选择算法得到的子集，利用SVM构建的预测器的预测精度。判断标准：（1）当

15、Sig(S,PCn)0，表明增加PCn这列属性对模型预测性能的提高有积极的作用，其可以融合为多重属性并利用预测器进行下一步判断。且Sig(S,PCn)越大表明效果越好，属性重要度越高。（2）当Sig(S,PCn)ACC ACC=max_ACC(i).acc; row=i; endendsave select_result2_3 ACC row2. 功能函数 random_5fold.m 代码：function Acc = random_5fold( Pdata, Ndata, data_atrr)Pnum=size(Pdata,1);Nnum=size(Ndata,1);indices1 =c

16、rossvalind(Kfold,Pnum,5); indices2 =crossvalind(Kfold,Nnum,5); Acc=;for i=1:5 Ptest = (indices1 = i); Ptrain = Ptest; Ntest = (indices2 = i); Ntrain = Ntest; Ztrain=Ptrain;Ntrain; Ztest=Ptest;Ntest; % 选出数据集训练属性和测试属性 train_attr=; test_attr=; train_attr =train_attr;data_atrr(find(Ztrain=1),:); test_at

17、tr =test_attr;data_atrr(find(Ztrain=0),:); % 训练集和测试集中的正样本 Ptrain_data=; Ptest_data=; Ptrain_data=Ptrain_data; Pdata(find(Ptrain=1); Ptest_data=Ptest_data; Pdata(find(Ptrain=0); % 训练集和测试集中的负样本 Ntrain_data=; Ntest_data=; Ntrain_data =Ntrain_data; Ndata(find(Ntrain=1); Ntest_data =Ntest_data; Ndata(fin

18、d(Ntrain=0); % 训练集的标签 train_label = vertcat( Ptrain_data, Ntrain_data); % 测试集的标签 test_label = vertcat( Ptest_data, Ntest_data); % 通过SVM构建的预测器，训练集训练模型 model = train(train_label,sparse(train_attr); % 通过SVM构建的预测器，测试集预测精度predict_label,accuracy,values=predict(test_label,sparse(test_attr),model); Acc=Acc,

19、accuracy(2,1);endend附件2：问题3的代码1.问题3中主函数problem2.m的MATLAB代码：close all;clear;clc;% 加载十进制编码9445个位点的所有属性load nowenary_encoding_feature% 读取对应文件下的300个基因path=dir(F:数模比赛Bproblem3gene_info*.dat); paths, ind = sort(path(:).datenum, ascend);path = path(ind);Pdata=zeros(500,1);Pnum=size(Pdata,1);Ndata=ones(500,

20、1);Nnum=size(Ndata,1);L=size(path,1);count=1;% 循环每个基因，利用每列属性进行该疾病的预测，得到预测精度for i=1:L filename = path(i).name; name = importdata(filename); num=size(name,1); data_attr1=all_feature(:,count:count+num-1); data_attr2=horzcat(data_attr,data_attr1); % 通过计算每个基因预测该疾病的准确率accuracy，取3次，取平均值 K=3; org_accuracy=z

21、eros(K,1); for j=1:K Acc = random_5fold( Pdata, Ndata,data_attr2); org_accuracy(j,1)=mean(Acc); end accuracy=mean(org_accuracy); class(i).acc=accuracy; class(i).name=filename; count=count+num;endsave liblinear_result_gene1_2 class% 选出最大的预测精度，及该最大预测精度所对应的基因ACC=-1;L=size(class,2);for i=1:L if class(i)

22、.accACC ACC=class(i).acc; row=i; endendsave select_result_gene1_2 ACC row2. 功能函数 random_5fold.m 代码：function Acc = random_5fold( Pdata, Ndata, data_atrr)Pnum=size(Pdata,1);Nnum=size(Ndata,1);indices1 =crossvalind(Kfold,Pnum,5); indices2 =crossvalind(Kfold,Nnum,5); Acc=;for i=1:5 Ptest = (indices1 = i

23、); Ptrain = Ptest; Ntest = (indices2 = i); Ntrain = Ntest; Ztrain=Ptrain;Ntrain; Ztest=Ptest;Ntest; % 选出数据集训练属性和测试属性 train_attr=; test_attr=; train_attr =train_attr;data_atrr(find(Ztrain=1),:); test_attr =test_attr;data_atrr(find(Ztrain=0),:); % 训练集和测试集中的正样本 Ptrain_data=; Ptest_data=; Ptrain_data=Pt

24、rain_data; Pdata(find(Ptrain=1); Ptest_data=Ptest_data; Pdata(find(Ptrain=0); % 训练集和测试集中的负样本 Ntrain_data=; Ntest_data=; Ntrain_data =Ntrain_data; Ndata(find(Ntrain=1); Ntest_data =Ntest_data; Ndata(find(Ntrain=0); % 训练集的标签 train_label = vertcat( Ptrain_data, Ntrain_data); % 测试集的标签 test_label = vertc

25、at( Ptest_data, Ntest_data); % 通过SVM构建的预测器，训练集训练模型 model = train(train_label,sparse(train_attr); % 通过SVM构建的预测器，测试集预测精度predict_label,accuracy,values=predict(test_label,sparse(test_attr),model); Acc=Acc,accuracy(2,1);endend附件3：问题4的代码1.main.m主程序的代码close all;clear;clc;% 加载所有位点的名称load SNPs_name% 加载 multi

26、z_phenos.txt 中的混合标签filename = .multi_phenos.txt;data=textread(filename);str=size(data,2);% 加载十进制编码9445个位点的所有属性load nowenary_encoding_featuretrain_attr=all_feature;% 属性数据归一化dataset_scale,ps = mapminmax(train_attr,0,1);data_attr1 = dataset_scale;W=size(data_attr1,2);for s=1:str% 加载所有样本对应的标签 Data=data(

27、:,s);Pdata=Data(1:500,:);Pnum=size(Pdata,1);Ndata=Data(1:500,:);Nnum=size(Ndata,1);for w=1:W K=5; org_accuracy=zeros(K,1); for i=1:K Acc = random_5fold( Pdata, Ndata,data_attr1(:,w); org_accuracy(i,1)=mean(Acc); end accuracy=mean(org_accuracy); max_ACC(w).acc=accuracy;endACC=-1;for i=1:W if max_ACC(i).accACC ACC=max_ACC(i).acc; row=i; endend multis.ACC=ACC;multis.row=row;=gene(row);endsave multi_result multi2. 功能函数 random_5fold.m 代码：function Acc = random_5fold( Pdata, Ndata, data_atrr)Pnum=size(Pdata,1);Nnum=size(Ndata,1);indi