优雅的使用matlab机器学习大作业 疾病miRNA关联分析personalrank方法

优雅的使用matlab机器学习大作业 疾病miRNA关联分析 personalrank方法

摘要

已知部分疾病和miRNA之间的关联性，本文基于MATLAB通过Personalrank算法对关联性进行预测并计算AUC曲线，对结果进行分析和评价。此法得到的auc值只有0.82，小白能力有限，欢迎大家指导并指出错误。

题目描述

1、 数据集：数据集中提供了疾病和miRNA名称、相似性数据（已经整理好）和已知的miRNA-疾病关联

2、 读取数据，自学matlab文件读取方式

3、 构建关联矩阵或者样本矩阵

4、 采用一种预测方法：线性回归、随机游走、推荐算法预测所有未知关联有关联的可能性

5、 绘制ROC曲线，计算AUC值

6、 评分标准：计算出来的AUC值*100就是你的得分

解决思路

1、通过load读取数据集

2、构建miRNA和疾病的关联矩阵

3、采用10折交叉验证，即去除关联矩阵中的某些值并使用训练后的模型预测被去除的值

4、根据预测结果绘制ROC曲线并计算auc值

详细步骤

1、通过load读取数据集

%数据读取
mRNA = load('data\mirsim.txt');%495个mRNA之间的关联
diseaseim = load('data\diseasesim.txt');%383个疾病之间的关联
mi_disease = load('data\Known miRNA-disease association number.txt');%mRNA与疾病的关联

2、构建miRNA和疾病的关联矩阵

这里构建了一个大小为$(m+n)\ast (m+n)$的对称矩阵$A$，分为了四个部分

$\left[\begin{array}{cc}P& R\\ R^{\prime }& Q\end{array}\right]$

其中$PQ$分别为miRNA和miRNA、疾病和疾病之间的关联，而$R$是miRNA和疾病之间的关联，有关联为1无关联为0。举个栗子，第$i$个miRNA和第$j$个miRNA的关联性为$A(i,j)$或$A(j,i)$，第$i$个miRNA和第$j$个疾病之间的关联性为$A(i,m+j)$或$A(m+j,i)$，以此类推。

m = size(mRNA,2);%m为mRNA的数量，n为疾病的数量
n = size(diseaseim,2);
A = [mRNA,zeros(m,n);
    zeros(n,m),diseaseim];%构建关联矩阵

for i=1:size(mi_disease,1)%mRNA与疾病的关联
    row = mi_disease(i,1);
    col = mi_disease(i,2)+m;
    A(row,col) = 1;
    A(col,row) = 1;
end

3、随机游走的迭代

根据随机游走算法，迭代的时候关联矩阵$A$不是0和1，而是百分比（即一行中1占的比例），故将$A$的出度进行修改。这里参考用PersonalRank实现基于图的推荐算法_peoplerank算法实现过程_HarryHuang1990的博客-CSDN博客的公式

$\begin{array}{l}\mathrm{PR}(i)=(1-\mathrm{d})r_i+d{\sum }_{j\in in(i)}\frac{PR(j)}{|out(i)|}\end{array}$

$r_i=\{\begin{array}{ll}1& i=u\\ 0& i\ne u\end{array}$

虽然我没看懂但我还是得随便扯一下，有机会一定补

下面是作者的原话：

---

我们发现PersonalRank跟PageRank的区别只是用替换了1/N，也就是说从不同点开始的概率不同。u表示我们推荐的目标用户，这样使用上式计算的就是所有顶点相对于顶点u的相关度。

与PageRank随机选择一个点开始游走（也就是说从每个点开始的概率都是相同的）不同，如果我们要计算所有节点相对于用户u的相关度，则PersonalRank从用户u对应的节点开始游走，每到一个节点都以1-d的概率停止游走并从u重新开始，或者以d的概率继续游走，从当前节点指向的节点中按照均匀分布随机选择一个节点往下游走。这样经过很多轮游走之后，每个顶点被访问到的概率也会收敛趋于稳定，这个时候我们就可以用概率来进行排名了。

在执行算法之前，我们需要初始化每个节点的初始概率值。如果我们对用户u进行推荐，则令u对应的节点的初始访问概率为1，其他节点的初始访问概率为0，然后再使用迭代公式计算。而对于pageRank来说，由于每个节点的初始访问概率相同，所以所有节点的初始访问概率都是1/N （N是节点总数）。
————————————————
版权声明：本文为CSDN博主「HarryHuang1990」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/HarryHuang1990/article/details/10048383

---

说白了，对于第一个miRNA，我们构建
$v_0={\left[\begin{array}{c}1\\ 0\\ 0\\ ⋮\\ 0\end{array}\right]}_{(1,m+n)}$
并用我们刚刚得到的$A$矩阵进行迭代
$v_1=v_0$
$v_1=\alpha A{v}_1+(1-\alpha )v_0$
直至$v_1$收敛，这样我们就得到了第一个miRNA和其他miRNA以及疾病的关联度$v_1$。这里$\alpha$即重启概率。
同理对第二第三直到最后一个miRNA和疾病重复以上操作……
利用MATLAB的矩阵运算简化过程，就像代码里的一样，构建一个$v_0$单位矩阵，这样$v_0$和$v_1$的每一列就是上面迭代过程的一个向量。迭代完成后，得到的$v_1$就是新的关联矩阵了。之所以还有一个$v_2$是为了保证$v_1$已经收敛。最后得到如下代码

orig_A=A;%保留原来的关联矩阵
for i=1:size(A,1)
    A(i,:)=A(i,:)/sum(A(i,:));%设置出度
end
v0=eye(m+n);
v1=v0;
alpha=0.2;
%========================
%这里用k折交叉验证去除部分已知的数据后，A矩阵变成了t
%========================
v2=alpha*t*v1+(1-alpha)*v0;
while sum(sum(abs(v2-v1)))>1e-15
    v1=v2;
    v2=alpha*t*v1+(1-alpha)*v0;
end

4、k折交叉验证

k折交叉验证即搞个大的循环把上边的迭代过程重复多次，这里展示一下k折交叉验证中的去除部分数据的代码。先看这段代码

kf=10;
indices = crossvalind('Kfold',(m+n)*(m+n),kf);
for i=1:kf
	t=A(:);%先变成行向量
    t(indices==i)=0;%把部分数据设置成未知的
end

k折交叉验证会把数据随机分成10份，t(indices==1)=0就是把第一份的数据全部变为0，理解了这里后，下面是交叉验证的迭代代码

kf=10;
indices = crossvalind('Kfold',(m+n)*(m+n),kf);
for i=1:kf
    t=A(:);
    t(indices==i)=0;%把部分数据设置成未知的
    t=reshape(t,m+n,m+n);
    v0=eye(m+n);
    v1=v0;
    alpha=0.2;
    v2=alpha*t*v1+(1-alpha)*v0;
    while sum(sum(abs(v2-v1)))>1e-15
        v1=v2;
        v2=alpha*t*v1+(1-alpha)*v0;
    end
   t=v2;
   t=t(:);
   predicts(indices==i)=t(indices==i);
   %去除部分的数据才是模型预测的，故将这部分数据保存在predicts里
end

k折交叉验证每一折预测未知的$1/10$的数据（折数不一定要是10，也可以更多）最后得到完整的关联矩阵。对于得到的结果，我们只用抽取$\left[\begin{array}{cc}P& R\\ R^{\prime }& Q\end{array}\right]$矩阵的$R$部分就是miRNAh和疾病的关联矩阵了，然后计算auc，绘制图像想偷懒就不写了。

5、完整代码

clear
clc
tic%该代码运行大约需要5分钟

%数据读取
mRNA = load('data\mirsim.txt');%495个mRNA之间的关联
diseaseim = load('data\diseasesim.txt');%383个疾病之间的关联
mi_disease = load('data\Known miRNA-disease association number.txt');%mRNA与疾病的关联

kf=10;
m = size(mRNA,2);%m为mRNA的数量，n为疾病的数量
n = size(diseaseim,2);
indices = crossvalind('Kfold',(m+n)*(m+n),kf);    

A = [mRNA,zeros(m,n);
    zeros(n,m),diseaseim];%构建关联矩阵

for i=1:size(mi_disease,1)%mRNA与疾病的关联
    row = mi_disease(i,1);
    col = mi_disease(i,2)+m;
    A(row,col) = 1;
    A(col,row) = 1;
end
orig_A=A;
for i=1:size(A,1)
    A(i,:)=A(i,:)/sum(A(i,:));%设置出度
end
predicts=zeros(m+n,m+n);
predicts=predicts(:);

%k折交叉验证
for i=1:kf
    t=A(:);
    t(indices==i)=0;%把部分数据设置成未知的
    t=reshape(t,m+n,m+n);
    v0=eye(m+n);
    v1=v0;
    alpha=0.2;
    v2=alpha*t*v1+(1-alpha)*v0;
    while sum(sum(abs(v2-v1)))>1e-15
        v1=v2;
        v2=alpha*t*v1+(1-alpha)*v0;
    end
   t=v2;
   t=t(:);
   predicts(indices==i)=t(indices==i);
end
fin_pred=reshape(predicts,m+n,m+n);
fin_A=orig_A(m+1:end,1:m);
fin_pred=fin_pred(m+1:end,1:m);
[PX, PY, Auc] = calculate_roc(fin_pred(:),fin_A(:));
plot(PX,PY)
disp(Auc)
toc

%这个函数用来绘制auc曲线，可略过
function [PX,PY,Auc] = calculate_roc(predict, ground_truth)
 
    pos_num = sum(ground_truth == 1);
    neg_num = sum(ground_truth == 0);

    m = length(ground_truth);
    [~, index] = sort(predict);
    ground_truth = ground_truth(index);
    PX = zeros(m+1,1);
    PY = zeros(m+1,1);
    Auc = 0;
    PX(1) = 1; PY(1) = 1;

    for i = 2:m
        TP = sum(ground_truth(i:m)==1);
        FP = sum(ground_truth(i:m)==0);
        PX(i) = FP/neg_num;
        PY(i) = TP/pos_num;
        Auc = Auc + (PY(i)+PY(i-1))*(PX(i-1)-PX(i))/2;     % 梯形面积：（上底+下底）*高/2
    end
    PX(m+1) = 0;
    PY(m+1) = 0;
    Auc = Auc + PY(m)*PX(m)/2;
end

6、结果

auc为0.8243