这是模式识别课上的一个作业，本人花了一点时间做了一下，发现分类器的学习结果勉强还能接受，在这里分享一下。
关于CpG岛：我们知道，生物的基因组是由A,T,G,C四种核苷酸组成的序列。人的一套基因组单链是由约30亿个A,T,G,C组成的一个超长序列，可以把它看成是一本由四个字母写成的天书。这样一个超长的字符串或其中一个子串，并不是由这四个字母随机组成的，而是遵循着很多特殊、复杂的规律。比如,由于一些生物化学机制的作用,基因组同一条链上出现相连的C
和G的概率要比随机情况小很多。把相连的C和G叫做一个CpG双核苷酸。人们已经观察到,CpG在基因组上出现的平均率比根据C、G各自出现的频率估计的组合出现频率小很多,但是,这些有限的CpG在基因组上分布的位置不是均匀的,而是倾向于集中在相对较短的一些片段上,这种CpG相对富集的区域被称作CpG岛,就像大海上的小岛一样。CpG岛在基因组上有重要的功能,研究CpG岛的识别是非常有意义的。CpG岛识别的基本问题是,给定一段DNA序列,判断它是否来自CpG岛,即判断这段序列样本是否是CpG岛的一部分。这可以看作是一个两类的分类问题,两类分别是CpG岛和非CpG岛,我们所能利用的特征就是序列本身有很多方法来预测或确定CpG岛,这里我们用一种简单的基于马尔可夫模型( MarkovModel)的方法来说明在这种场景下用统计决策来进行分类的一种思想。
当我们独立地考虑DNA序列每个位置上的核苷酸时,我们可以把它当作一个有四种取值的离散随机变量x={A,T,G,C},面前面例子中我们考虑的都是x∈R的连续情况,在连续情况下,我们用概率密度函数来表示变量取值的分布,而在离散情况下,我们则用变量取各个值的概率来表示P(x),显然,P(x=A)+P(x=T)+P(x=G)+P(x=C)=1。在连续情况下,如果我们有多个特征,则用随机向量x来表示,特征之间的关系反映在随机向量的联合概率密度上。在离散情况下,如果序列的每个位置上核苷酸的分布是独立同分布的,那么每个位置就是随机变量的一次实现。但是,当我们考虑在连续的位置上出现的CpG双核苷酸时,这种模型就不适用了,因为两个位置不再独立。我们可以用马尔可夫模型来表示这种相邻位置之间的依賴关系。有很多专门的教材和文献讲述马尔可夫模型,我们这里只是给出一种直观的简单描述。我们把DNA序列看作符号组成的时间序列x1,x2, ..., xL,每一时刻(位置)上的取值为xi={A,T,G,C},如果第i时刻上的取值依赖于且仅依赖于第i-1时刻的取值,即P(xi|xi-1,xi-2, ... ,x1)=P(xi|xi-1)。则把这个串称作一个一阶马尔可夫链或一阶马尔可夫模型。如果第i时刻上的取值依赖且仅依赖于第i-k到到第i-1时刻的取值,则这个串就是k阶马尔可夫链。在考虑CpG岛时,我们只需同时考虑相邻两个位置上的核苷酸,因此可以用一阶马尔可夫链来描述。马尔可夫链可以用条件概率模型来描述。（这些概念可以参考张学工教授的《模式识别》第三版，由于一些数学公式不方便打出来，所以这些概念就只介绍到这里，现在就假设你已经知道了一些概念，如果你还不知道这些概念，相信你看完我这篇文章之后依旧还是懵逼状态）
其中有一些名词你要知道它的含义：转移概率，马尔科夫链，状态转移矩阵，正样本，负样本，测试集，结果集，学习器，分类器，对数几率比，对数似然比矩阵（可不做要求），还有到后面分析用的：NP决策，阈值，ROC曲线。现在也假设这些概念和所用到的数学公式你都已经在书上看懂了。好了，前面科普了一大堆，真怕有人骂我发这种科普类的文章，现在就进入正题：
现在我们要做的事情：给你一个正样本数据集（CpG岛），和一个负样本数据集（非CpG岛），我后面会提供这些较为庞大的数据集的文件（正样本为2000多个数据集，负样本约为8万个数据集），然后让你通过马尔科夫模型编程写出你的学习器的算法。其中，你要算出正负样本的状态转移矩阵，然后计算出对应灵敏度sn和特异度sp（这里的阈值取0，为什么？自己想（对数几率比）），然后分别取100个正负样本，测试这样样本的分类结果如何，然后，你再根据测试结果集等区间划分阈值（当然你也可以选择其他的方法划分阈值），然后在以每个阈值对这100个正负样本进行测试，每一次测的结果集都把相应的sn和sp算出来，最后得到的是一系列的sn,sp一一对应的值，然后以sn为纵坐标，以1-sp为横坐标，画出曲线，得到的就是传说中的ROC曲线。你可以根据roc曲线对这个学习器的好坏进行评价。在这里可以比较难以理解的是对划分阈值后在进行测试的理解，还有就是roc曲线的本质，这些东西如果我全部在一篇文章里面陈述那就太庞大了，可以参考其他文章，然后你才可能看懂我写的是什么意思。我这里是用C++编写的源程序，对于C++不太熟悉的人来说，读懂源代码可能有一些难度，但是希望你能把一些逻辑和思路弄懂，而不要纠结在代码的细节上。当然我希望大家使用MATLAB或者Python。
最后想说一点，这个东西确实很有趣，也很强大，如果你是模式识别或者机器学习的一个入门者，请你把我以上提到的这些东西都看懂，这样会对你有很多的帮助。（再一次声明：这篇文章所涉及的知识量有点大，你要提前知道这些知识再来阅读！）
下面送上C++源程序代码：（可以拷贝运行看结果！）
// Text_C++.cpp : 定义控制台应用程序的入口点。
#include "stdafx.h"
#include<cmath>
#include<iostream>
#include <fstream>
#include <string>
#include<sstream>
#include<stdlib.h>
#include<vector>
#include<algorithm>
using namespace std;
inline void NEXT(const string&T, vector<int>&next){//按模式串生成vector,next(T.size())
next[0] = -1;
for (int i = 1; i<T.size(); i++){
int j = next[i - 1];
while (j >= 0 && T[i - 1] != T[j]) j = next[j];//递推计算
if (j >= 0 && T[i - 1] == T[j]) next[i] = j + 1;
else next[i] = 0;
}
}
inline string::size_type COUNT_KMP(const string&S, const string&T){
//利用模式串T的next函数求T在主串S中的个数count的KMP算法
//其中T非空，
vector<int>next(T.size());
NEXT(T, next);
string::size_type index, count = 0;
for (index = 0; index<S.size(); ++index){
int pos = 0;
string::size_type iter = index;
while (pos<T.size() && iter<S.size()){
if (S[iter] == T[pos]){ ++iter; ++pos; }
else{
if (pos == 0) ++iter;
else pos = next[pos - 1] + 1;
}
}
if (pos == T.size() && (iter - index) == T.size()) ++count;
}
return count;
}
void funt(double text_t[],double yuzhi,double *tp,double *fn){
// int j=0;
for(int i=0;i<100;i++)
{
if(text_t[i]>=yuzhi)
*tp=*tp+1;
else
*fn=*fn+1;
}
}
void funf(double text_f[],double yuzhi,double *tn,double *fp){
//int j=0;
for(int i=0;i<100;i++)
{
if(text_f[i]<=yuzhi)
*tn=*tn+1;
else
*fp=*fp+1;
}
} //其实上面的两个函数是一样的！只是传参数不同
int main(int argc, char*argv[])
{
char strue[100][18];
char sfalse[100][18];
//char aaa='dda';
double sum=0;
double aa,ca,ga,ta,ac,cc,gc,tc;
double ag,cg,gg,tg,at,ct,gt,tt;
double suma,sumc,sumg,sumt;
double pa,pc,pg,pt;
double sumag=0,sumcg=0,sumgg=0,sum_tg=0;
double sumta=0,sumtc=0,sumtg_=0,sumtt=0;
double sumag1=0,sumcg1=0,sumgg1=0,sum_tg1=0;
double sumta1=0,sumtc1=0,sumtg_1=0,sumtt1=0;
//已经通过本程序算出来的负样本状态转移矩阵，在此直接使用计算值
double aa1=0.298,ca1=0.264,ga1=0.231,ta1=0.216;
double ac1=0.226,cc1=0.268,gc1=0.229,tc1=0.278;
double ag1=0.341,cg1=0.078,gg1=0.312,tg1=0.270;
double at1=0.221,ct1=0.233,gt1=0.191,tt1=0.356;
double two=2.00;
double text_t[100]={0}; //记录测试结果
double text_f[100]={0};
double text_add[200];
double yuzhi[21];
double tph=0,fph=0,fnh=0,tnh=0,snh=0,sph=0;
double tp[21]={0},fp[21]={0},fn[21]={0},tn[21]={0},sp[21]={0},sn[21]={0};
//定义模式串
string Taa="aa",Tca="ca",Tga="ga",Tta="ta",Tac="ac",Tcc="cc",Tgc="gc",Ttc="tc";
string Tag="ag",Tcg="cg",Tgg="gg",Ttg="tg",Tat="at",Tct="ct",Tgt="gt",Ttt="tt";
//定义数量
string::size_type countaa,countca,countga,countta,countac,countcc,countgc,counttc;
string::size_type countag,countcg,countgg,counttg,countat,countct,countgt,counttt;
string sline,sline2;//每一行
string ss="",sf="";
string s1,s2;
int j=0,m=0;
ifstream inft,inff;
//打开文件
inft.open("D://truesamples.txt",ios::in);
inff.open("D://falsesamples.txt",ios::in);
//遍历文件
while(getline(inft,sline))
{
j=j+1;
int i=0;
istringstream sin(sline);
sin>>s1;//取出行给s1赋值
if(j<=100)
{
for(int k=0;k<18;k++){
strue[j-1][k]=s1.data()[k]; //存入二维数组
}
}
string::iterator it;
for (it =s1.begin(); it != s1.end(); ++it)
{
i=i+1;
if(i==9)
{
if (*it =='a')
sumag++;
if(*it=='c')
sumcg++;
if(*it=='g')
sumgg++;
else
sum_tg++;
}
if(i==12)
{
if (*it =='a')
sumta++;
if(*it=='c')
sumtc++;
if(*it=='g')
sumtg_++;
else
sumtt++;
}
}
ss=ss+s1; //把s1加到ss上
sum++;
}
while(getline(inff,sline2))
{
istringstream sin(sline2);
sin>>s2;//取出行给s1赋值
m=m+1;
if(m<=100)
{
for(int k=0;k<18;k++){
sfalse[m-1][k]=s2.data()[k]; //存入二维数组
}
}
if(m>100)
break;
}
//分别调用KMP函数求组合数量
countaa = COUNT_KMP(ss,Taa),countca=COUNT_KMP(ss,Tca),countga=COUNT_KMP(ss,Tga);
countta = COUNT_KMP(ss,Tta)-sumta,countac=COUNT_KMP(ss,Tac),countcc=COUNT_KMP(ss,Tcc);
countgc = COUNT_KMP(ss,Tgc),counttc=COUNT_KMP(ss,Ttc)-sumtc,countag=COUNT_KMP(ss,Tag)-sumag;
countcg = COUNT_KMP(ss,Tcg)-sumcg,countgg=COUNT_KMP(ss,Tgg)-sumgg,counttg=COUNT_KMP(ss,Ttg)-sum_tg-sumtg_;
countat = COUNT_KMP(ss,Tat),countct=COUNT_KMP(ss,Tct),countgt=COUNT_KMP(ss,Tgt)-sum;
counttt=COUNT_KMP(ss,Ttt);
//求分别以a,c,g,t为后一个元素的概率
suma=countaa+countca+countga+countta;
sumc=countac+countcc+countgc+counttc;
sumg=countag+countcg+countgg+counttg;
sumt=countat+countct+countgt+counttt;
//矩阵求值
aa=countaa/suma,ca=countca/suma,ga=countga/suma,ta=countta/suma;
ac=countac/sumc,cc=countcc/sumc,gc=countgc/sumc,tc=counttc/sumc;
ag=countag/sumg,cg=countcg/sumg,gg=countgg/sumg,tg=counttg/sumg;
at=countat/sumt,ct=countct/sumt,gt=countgt/sumt,tt=counttt/sumt;
//矩阵输出
cout<<aa<<" "<<ca<<" "<<ga<<" "<<ta<<endl;
cout<<ac<<" "<<cc<<" "<<gc<<" "<<tc<<endl;
cout<<ag<<" "<<cg<<" "<<gg<<" "<<tg<<endl;
cout<<at<<" "<<ct<<" "<<gt<<" "<<tt<<endl;
//验证每一行的概率之和是否为1
pa=aa+ca+ga+ta;
pc=ac+cc+gc+tc;
pg=ag+cg+gg+tg;
pt=at+ct+gt+tt;
cout<<pa<<" "<<pc<<" "<<pg<<" "<<pt<<endl; //验证概率是否正确
for(int i=0;i<100;i++){
for(int j=1;j<9;j++){
if(strue[i][j]=='a'){
if(strue[i][j-1]=='a'){
text_t[i]+=log10(aa/aa1)/log10(two);
}
if(strue[i][j-1]=='c'){
text_t[i]+=log10(ca/ca1)/log10(two);
}
if(strue[i][j-1]=='g'){
text_t[i]+=log10(ga/ga1)/log10(two);
}
if(strue[i][j-1]=='t'){
text_t[i]+=log10(ta/ta1)/log10(two);
}
}
if(strue[i][j]=='c'){
if(strue[i][j-1]=='a'){
text_t[i]+=log10(ac/ac1)/log10(two);
}
if(strue[i][j-1]=='c'){
text_t[i]+=log10(cc/cc1)/log10(two);
}
if(strue[i][j-1]=='g'){
text_t[i]+=log10(gc/gc1)/log10(two);
}
if(strue[i][j-1]=='t'){
text_t[i]+=log10(tc/tc1)/log10(two);
}
}
if(strue[i][j]=='g'){
if(strue[i][j-1]=='a'){
text_t[i]+=log10(ag/ag1)/log10(two);
}
if(strue[i][j-1]=='c'){
text_t[i]+=log10(cg/cg1)/log10(two);
}
if(strue[i][j-1]=='g'){
text_t[i]+=log10(gg/gg1)/log10(two);
}
if(strue[i][j-1]=='t'){
text_t[i]+=log10(tg/tg1)/log10(two);
}
}
if(strue[i][j]=='t'){
if(strue[i][j-1]=='a'){
text_t[i]+=log10(at/at1)/log10(two);
}
if(strue[i][j-1]=='c'){
text_t[i]+=log10(ct/ct1)/log10(two);
}
if(strue[i][j-1]=='g'){
text_t[i]+=log10(gt/gt1)/log10(two);
}
if(strue[i][j-1]=='t'){
text_t[i]+=log10(tt/tt1)/log10(two);
}
}
}
for(j=12;j<18;j++){
if(strue[i][j]=='a'){
if(strue[i][j-1]=='a'){
text_t[i]+=log10(aa/aa1)/log10(two);
}
if(strue[i][j-1]=='c'){
text_t[i]+=log10(ca/ca1)/log10(two);
}
if(strue[i][j-1]=='g'){
text_t[i]+=log10(ga/ga1)/log10(two);
}
if(strue[i][j-1]=='t'){
text_t[i]+=log10(ta/ta1)/log10(two);
}
}
if(strue[i][j]=='c'){
if(strue[i][j-1]=='a'){
text_t[i]+=log10(ac/ac1)/log10(two);
}
if(strue[i][j-1]=='c'){
text_t[i]+=log10(cc/cc1)/log10(two);
}
if(strue[i][j-1]=='g'){
text_t[i]+=log10(gc/gc1)/log10(two);
}
if(strue[i][j-1]=='t'){
text_t[i]+=log10(tc/tc1)/log10(two);
}
}
if(strue[i][j]=='g'){
if(strue[i][j-1]=='a'){
text_t[i]+=log10(ag/ag1)/log10(two);
}
if(strue[i][j-1]=='c'){
text_t[i]+=log10(cg/cg1)/log10(two);
}
if(strue[i][j-1]=='g'){
text_t[i]+=log10(gg/gg1)/log10(two);
}
if(strue[i][j-1]=='t'){
text_t[i]+=log10(tg/tg1)/log10(two);
}
}
if(strue[i][j]=='t'){
if(strue[i][j-1]=='a'){
text_t[i]+=log10(at/at1)/log10(two);
}
if(strue[i][j-1]=='c'){
text_t[i]+=log10(ct/ct1)/log10(two);
}
if(strue[i][j-1]=='g'){
text_t[i]+=log10(gt/gt1)/log10(two);
}
if(strue[i][j-1]=='t'){
text_t[i]+=log10(tt/tt1)/log10(two);
}
}
}
}
for(int i=0;i<100;i++){
for(int j=1;j<9;j++){
if(sfalse[i][j]=='a'){
if(sfalse[i][j-1]=='a'){
text_f[i]+=log10(aa/aa1)/log10(two);
}
if(sfalse[i][j-1]=='c'){
text_f[i]+=log10(ca/ca1)/log10(two);
}
if(sfalse[i][j-1]=='g'){
text_f[i]+=log10(ga/ga1)/log10(two);
}
if(sfalse[i][j-1]=='t'){
text_f[i]+=log10(ta/ta1)/log10(two);
}
}
if(sfalse[i][j]=='c'){
if(sfalse[i][j-1]=='a'){
text_f[i]+=log10(ac/ac1)/log10(two);
}
if(sfalse[i][j-1]=='c'){
text_f[i]+=log10(cc/cc1)/log10(two);
}
if(sfalse[i][j-1]=='g'){
text_f[i]+=log10(gc/gc1)/log10(two);
}
if(sfalse[i][j-1]=='t'){
text_f[i]+=log10(tc/tc1)/log10(two);
}
}
if(sfalse[i][j]=='g'){
if(sfalse[i][j-1]=='a'){
text_f[i]+=log10(ag/ag1)/log10(two);
}
if(sfalse[i][j-1]=='c'){
text_f[i]+=log10(cg/cg1)/log10(two);
}
if(sfalse[i][j-1]=='g'){
text_f[i]+=log10(gg/gg1)/log10(two);
}
if(sfalse[i][j-1]=='t'){
text_f[i]+=log10(tg/tg1)/log10(two);
}
}
if(sfalse[i][j]=='t'){
if(sfalse[i][j-1]=='a'){
text_f[i]+=log10(at/at1)/log10(two);
}
if(sfalse[i][j-1]=='c'){
text_f[i]+=log10(ct/ct1)/log10(two);
}
if(sfalse[i][j-1]=='g'){
text_f[i]+=log10(gt/gt1)/log10(two);
}
if(sfalse[i][j-1]=='t'){
text_f[i]+=log10(tt/tt1)/log10(two);
}
}
}
for(j=12;j<18;j++){
if(sfalse[i][j]=='a'){
if(sfalse[i][j-1]=='a'){
text_f[i]+=log10(aa/aa1)/log10(two);
}
if(sfalse[i][j-1]=='c'){
text_f[i]+=log10(ca/ca1)/log10(two);
}
if(sfalse[i][j-1]=='g'){
text_f[i]+=log10(ga/ga1)/log10(two);
}
if(sfalse[i][j-1]=='t'){
text_f[i]+=log10(ta/ta1)/log10(two);
}
}
if(sfalse[i][j]=='c'){
if(sfalse[i][j-1]=='a'){
text_f[i]+=log10(ac/ac1)/log10(two);
}
if(sfalse[i][j-1]=='c'){
text_f[i]+=log10(cc/cc1)/log10(two);
}
if(sfalse[i][j-1]=='g'){
text_f[i]+=log10(gc/gc1)/log10(two);
}
if(sfalse[i][j-1]=='t'){
text_f[i]+=log10(tc/tc1)/log10(two);
}
}
if(sfalse[i][j]=='g'){
if(sfalse[i][j-1]=='a'){
text_f[i]+=log10(ag/ag1)/log10(two);
}
if(sfalse[i][j-1]=='c'){
text_f[i]+=log10(cg/cg1)/log10(two);
}
if(sfalse[i][j-1]=='g'){
text_f[i]+=log10(gg/gg1)/log10(two);
}
if(sfalse[i][j-1]=='t'){
text_f[i]+=log10(tg/tg1)/log10(two);
}
}
if(sfalse[i][j]=='t'){
if(sfalse[i][j-1]=='a'){
text_f[i]+=log10(at/at1)/log10(two);
}
if(sfalse[i][j-1]=='c'){
text_f[i]+=log10(ct/ct1)/log10(two);
}
if(sfalse[i][j-1]=='g'){
text_f[i]+=log10(gt/gt1)/log10(two);
}
if(sfalse[i][j-1]=='t'){
text_f[i]+=log10(tt/tt1)/log10(two);
}
}
}
}
cout<<endl<<endl<<"正样本预测:"<<endl;
for(int i=0;i<100;i++){
cout<<text_t[i]<<endl;
if(text_t[i]>=0)
tph++;
else
fnh++;
}
cout<<endl<<endl<<"负样本预测:"<<endl;
for(int i=0;i<100;i++){
cout<<text_f[i]<<endl;
if(text_f[i]>=0)
fph++;
else
tnh++;
}
cout<<endl<<endl;
cout<<"TP:"<<tph<<endl;
cout<<"FN:"<<fnh<<endl;
cout<<"FP:"<<fph<<endl;
cout<<"TN:"<<tnh<<endl;
cout<<endl;
snh=tph/(tph+fnh);
sph=tnh/(tnh+fph);
cout<<"Sn:"<<snh<<endl<<"Sp:"<<sph<<endl;
for(int i=0;i<100;i++){
text_add[i]=text_t[i];
text_add[i+100]=text_f[i];
}
sort(text_add,text_add+200);
//等排列区间划分样本测试集阈值
int jj=0;
for(int i=0;i<21;i++){
yuzhi[i]=text_add[jj];
jj=jj+10;
if(jj>199)
jj=jj-1;
}
for(int i=0;i<21;i++){
funt(text_t,yuzhi[i],&tp[i],&fn[i]); //调用函数求值
funf(text_f,yuzhi[i],&tn[i],&fp[i]);
}
//求灵敏度与特异度
for(int i=0;i<21;i++){
sn[i]=tp[i]/(tp[i]+fn[i]);
sp[i]=tn[i]/(tn[i]+fp[i]);
}
cout<<endl<<endl;
//输出结果
cout<<"1-sp"<<" "<<"sn"<<endl;
for(int i=0;i<21;i++){
cout<<1-sp[i]<<" "<<sn[i]<<endl;
}
inft.close();//关闭文件
inff.close();
return 0;
}
很显然，还能在做优化，但是我实在是不想在弄了，时间和精力，以及水平都有限，谢谢大家的理解。
样本链接：https://pan.baidu.com/s/1C6kICZhK8A3uNO7qMgJbhg
密码：pf9d
209
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
