理解这段代码需要对 MATLAB 编程语言、数学计算以及稀疏极大似然估计的相关概念有一定了解。由于代码相对长且包含复杂的数学计算，这里我将对其中的关键部分进行逐行解释。请注意，这是一个高度专业的代码，理解涉及到的数学概念可能需要深入的数学背景。
function[B,me]=sparse_maximum_likelihood_B(W,B,Y,K,Q,lambda_factor,lambda_factor_prev,sigma2,params)
if(nargin==8)
params=[];
end
if(isfield(params,'maxiter'))
maxiter=params.maxiter;
else
maxiter=1000;
end
```
- `function [B,me]=sparse_maximum_likelihood_B(W,B,Y,K,Q,lambda_factor,lambda_factor_prev,sigma2,params)`: 定义 MATLAB 函数 `sparse_maximum_likelihood_B`，该函数接受输入参数 `W, B, Y, K, Q, lambda_factor, lambda_factor_prev, sigma2, params`，并返回输出参数 `B, me`。
- `if(nargin==8) params=[]; end`: 如果输入参数的个数为 8，则将 `params` 初始化为空。
- `if(isfield(params,'maxiter')) maxiter=params.maxiter; else maxiter=1000; end`: 如果 `params` 结构体中包含字段 `maxiter`，则将其值赋给 `maxiter`；否则，默认设置 `maxiter` 为 1000。
```matlab
[M,N]=size(Y);
[Y,meanY]=centerY(Y,K);
```
- `[M,N]=size(Y)`: 获取矩阵 `Y` 的大小，`M` 为行数，`N` 为列数。
- `[Y,meanY]=centerY(Y,K)`: 调用 `centerY` 函数对矩阵 `Y` 进行中心化处理，并返回中心化后的矩阵 `Y` 和均值 `meanY`。
```matlab
IM=speye(M);
```
- `IM=speye(M)`: 创建大小为 MxM 的单位矩阵 `IM`。
```matlab
Ylamb=zeros(M);
for i=1:M
nontargs=K(i,:);
Yc=Y{i}(:,nontargs)*Y{i}(:,nontargs)';
Ylamb(i,:)=Yc(i,:);
end
```
- `Ylamb=zeros(M)`: 创建大小为 MxM 的零矩阵 `Ylamb`。
- `for i=1:M ... end`: 循环遍历 M 行，对 `Ylamb` 进行计算。
- `nontargs=K(i,:)`: 获取 `K` 矩阵的第 i 行。
- `Yc=Y{i}(:,nontargs)*Y{i}(:,nontargs)'`: 计算矩阵 `Y{i}(:,nontargs)` 的外积。
- `Ylamb(i,:)=Yc(i,:)`: 将计算结果的第 i 行赋给 `Ylamb`。
```matlab
lambda_max=max(max(abs((-sigma2*sum(K,2).*IM+Ylamb)./W)));
```
- `lambda_max=max(max(abs((-sigma2*sum(K,2).*IM+Ylamb)./W)))`: 计算正则化参数 `lambda` 的最大值。
```matlab
lambda=lambda_max*lambda_factor;
```
- `lambda=lambda_max*lambda_factor`: 计算正则化参数 `lambda`。
```matlab
Dlambda=(2*lambda_factor-lambda_factor_prev)*lambda_max*W;
S=(abs(Q)>=Dlambda);
```
- `Dlambda=(2*lambda_factor-lambda_factor_prev)*lambda_max*W`: 计算正则化参数的变化量。
- `S=(abs(Q)>=Dlambda)`: 创建一个与 `Q` 大小相同的矩阵 `S`，元素值为 `abs(Q) >= Dlambda` 的逻辑值。
```matlab
for i=1:M/2
S(i,i+M/2)=0;
end
```
- `for i=1:M/2 ... end`: 将 `S` 矩阵的前半部分对角线元素设为零。
```matlab
s=sum(S,2);
ei=zeros(M,1);
B=B.*S;
```
- `s=sum(S,2)`: 计算 `S` 矩阵每行元素的和。
- `ei=zeros(M,1)`: 创建大小为 Mx1 的零向量 `ei`。
- `B=B.*S`: 将 `B` 矩阵中对应于 `S` 中非零元素的位置上的元素保留，其他位置设为零。
```matlab
for k=1:maxiter
Bold=B;
for i=1:M
nontargs=K(i,:);
Nhat=sum(nontargs);
if (s(i))
ei(i)=1;
zi=(IM-B)\ei;
js_i=find(S(i,:)>0);
for j=js_i
m_ij=zi(j);
B_old=B(i,j);
if(j~=i)
lambdaW=lambda*W(i,j);
y_j=Y{i}(j,nontargs)';
BiT=B(i,:);
BiT(j)=0;
a_iT=(ei'-BiT)*Y{i}(:,nontargs);
r_ij=y_j'*y_j;
beta_ij=a_iT*y_j;
if(m_ij==0)
Bij=(beta_ij-lambdaW)/r_ij;
if Bij>0
B(i,j)=Bij;
else
Bij=(beta_ij+lambdaW)/r_ij;
if Bij<0
B(i,j)=Bij;
else
B(i,j)=0;
end
end
else
d_ij=1/m_ij+B(i,j);
theta_ijp=r_ij*d_ij+beta_ij-lambdaW;
k_ijp=d_ij*(beta_ij-lambdaW)-Nhat*sigma2;
q_ijp=theta_ijp^2-4*r_ij*k_ijp;
Bijpp=(1/(2*r_ij))*(theta_ijp+sqrt(q_ijp));
Bijpm=(1/(2*r_ij))*(theta_ijp-sqrt(q_ijp));
q_ijm=q_ijp +4*lambdaW*(beta_ij-r_ij*d_ij);
theta_ijm=theta_ijp+2*lambdaW;
Bijmm=(1/(
这段代码是一个用于网络推断的 MATLAB 脚本。以下是代码的主要步骤和解释：
```matlab
% Copyright 漏 2013 Cai et al; 漏 2019 North Carolina State University. All rights reserved.
% 鈥渓ignin_modelSML鈥� by Cai X, Bazerque JA, Giannakis GB and North Carolina State University is licensed under CC BY 3.0 (https://creativecommons.org/licenses/by/3.0/us/legalcode)
% rng_seed=1010;
% rng(rng_seed) %set seed for random number generator, for reproducibility
%%% Load transcript and protein abundances
load ../TranscriptProteinAbundances.mat % load lignin data
Ytable = Ymeanlines_imputed; % Mean of the experimental replicates
Xmask = Xmeanlines_mask; % indicates targeted transcripts
Y = Ytable{:,:}'; % Rows of Y are the transcripts and proteins, columns are the experiments
X = Xmask{:,:}';
[M, N] = size(Y); % M is the number of transcripts and proteins, N is the number of experiments
K = ~logical(X); % binary matrix size MxN, 1 when transcript/protein was not targeted in an experiment, 0 when it was
[Ysc, Asc] = scaleY(Y); % Scale each transcript and protein so that their abundances fall between 0 and 1
%% Set parameters
lambda_factors = 10.^(0:-.1:-3); % SML regularization parameters
params.Kcv = 5; % # of cross-validation folds
params.rho_factors = 10.^(-6:0.1:1); % Ridge regularization parameters
params.lambda_factors = lambda_factors;
params.maxiter = 1000; % max iterations for SML algorithm
params.cv_its = 5; % # of iterations of the kcv-fold cross-validations
%% Network inference
[ilambda_cv] = cross_validation_SML_B(Y, K, params); % Use cross-validation to determine SML regularization parameter
[rho_factor, sigma2] = cross_validation_ridge_B(Ysc, K, params); % Use cross-validation to determine ridge regularization parameter, and estimate sigma2 for SML
[BRhat, mueRhat] = constrained_ridge_B(Ysc, K, rho_factor); % Calculate ridge regression of Ysc using rho_factor from cross_validation
BR = diag(1./Asc) * BRhat * diag(Asc); % Un-scales calculated BRhat and mueRhat to use with Y
mueR = diag(1./Asc) * mueRhat;
W = 1./abs(BRhat); % Use BRhat to compute regularization weights to use with Ysc
Q = Inf; % Ignore discarding rules from SML
% Initialize variables for SML
lambda_factor_prev = 1;
BLhat = zeros(M);
for ilambda = 1:ilambda_cv
[BLhat, mueLhat] = sparse_maximum_likelihood_B(W, BLhat, Ysc, K, Q, lambda_factors(ilambda), lambda_factor_prev, sigma2, params);
Q = Inf;
lambda_factor_prev = lambda_factors(ilambda);
end % ilambda
SL = abs(sign(BLhat)); % Detected edges from SML algorithm
[BChat, mueChat] = constrained_ML_B(BLhat, SL, Ysc, K, sigma2, params); % Use constrainedML algorithm on just the identified edges. This solves for the weights of the identified edges without the l1-norm penalty
BC = diag(1./Asc) * BChat * diag(Asc); % Un-scale BChat and mueChat to use with Y
mueC = diag(1./Asc) * mueChat;
BCsparse = sparse(BC);
```
Lignin_modelSML.m
这段代码主要进行以下步骤：
1. **加载数据：** 从文件中加载转录本和蛋白质的丰度数据，以及指示目标转录本的掩码数据。
2. **数据处理：** 将加载的数据进行处理，获取转录本和蛋白质的丰度矩阵 `Y` 和目标掩码矩阵 `X`，其中 `K` 是非目标转录本的二进制矩阵。
3. **数据缩放：** 对 `Y` 进行归一化，使其丰度值在 [0, 1] 范围内。
4. **参数设置：** 设置网络推断所需的参数，包括 SML 的正则化参数、交叉验证折叠数、岭回归的正则化参数等。
5. **网络推断：** 利用 SML 和岭回归进行网络推断。通过交叉验证确定 SML 的正则化参数和岭回归的正则化参数，然后执行岭回归得到 `BRhat` 和 `mueRhat`。
6. **SML 算法：** 利用 SML 算法迭代计算 `BLhat` 和 `mueLhat`。
7. **网络边缘检测：** 根据 SML 算法得到的结果 `BLhat`，通过 l1-norm 惩罚将其稀疏化，得到 `SL`。
8. **Constrained ML 算法：** 利用 Constrained ML 算法在检测到的网络边缘上求解权重，得到 `BChat` 和 `mueChat`。
9. **结果缩放：** 对结果进行反缩放，得到原始数据空间中的 `BC`。
10. **稀疏表示：** 将 `BC` 转换为稀疏矩阵 `BCsparse`。
这段代码的主要目的是执行网络推断，找到相关性或连接性较强的转录本和蛋白质。
这个函数的目的是对输入矩阵 `Y` 进行缩放，使得每一行的元素范围都在 `[0, 1]` 之间。以下是函数的详细解释：
scaleY.m
function [Ysc, Ascaling] = scaleY(Y)
% Y is an M x N matrix with non-negative elements
% Scale Data so each row has the range [0 1]
% Find the maximum value along each row
Ymax = max(Y, [], 2);
% Calculate the scaling factor for each row
Ascaling = 1 ./ Ymax;
% Scale the data by multiplying each row by its corresponding scaling factor
Ysc = diag(Ascaling) * Y;
end
```
解释：
- `Y` 是一个大小为 M x N 的矩阵，其中的元素都是非负数。
- `Ymax` 计算了每一行的最大值。
- `Ascaling` 是一个向量，包含了每一行的缩放因子。每个缩放因子是 1 除以对应行的最大值。
- `diag(Ascaling)` 创建了一个对角矩阵，其中对角线上的元素是 `Ascaling` 中的元素。
- `Ysc` 是经过缩放的矩阵，通过将每一行乘以对应的缩放因子得到。
因此，函数的输出 `Ysc` 是一个每行元素范围在 `[0, 1]` 之间的矩阵，而 `Ascaling` 是用于还原缩放的缩放因子。
Sparse_maximum_likelihood_B.m
这是一个用于稀疏极大似然估计（Sparse Maximum Likelihood，SML）的 MATLAB 函数。下面是对代码的详细解释：
1. **输入参数解释**：
- `W`: 正则化惩罚项的权重矩阵。
- `B`: 待估计的稀疏系数矩阵。
- `Y`: 观测数据矩阵。
- `K`: 二进制矩阵，表示哪些观测是非目标观测。
- `Q`: 二进制矩阵，表示哪些系数需要被忽略。
- `lambda_factor`: 正则化参数。
- `lambda_factor_prev`: 先前的正则化参数。
- `sigma2`: 方差。
- `params`: 一些其他参数，是一个结构体，包含 `maxiter`（最大迭代次数）。
2. **函数主体解释**：
- `Y` 和 `meanY` 被标准化，其中 `meanY` 存储了每个观测的平均值。
- 计算了一些矩阵，如 `IM`（单位矩阵）和 `Ylamb`（根据观测数据生成的矩阵）。
- 计算了正则化参数 `lambda` 和 `Dlambda`。
- 通过迭代进行 SML 估计，其中更新了系数矩阵 `B`，并在每次迭代中更新了一些辅助矩阵和向量。
- 在每次迭代中，通过解线性或二次方程来更新 `B` 的值。
- 最后，计算 `delta_B`，如果小于阈值则停止迭代。
3. **输出解释**：
- `B`: 最终的系数矩阵。
- `me`: 与 `B` 相关的估计均值。
总体而言，这个函数实现了一种通过迭代求解 SML 问题的方法，其中在每次迭代中通过解线性或二次方程来估计系数矩阵。
centerY.m
这是一个用于中心化数据的 MATLAB 函数。下面是对代码的详细解释：
1. **输入参数解释**：
- `Y`：原始观测数据矩阵。
- `K`：二进制矩阵，表示哪些观测是非目标观测。
2. **函数主体解释**：
- `Y_bar` 是一个单元格数组，其中每个元素 `Y_bar{i}` 是对应观测的非目标观测的均值向量。
- `Yc` 是一个单元格数组，其中每个元素 `Yc{i}` 是对应观测减去其均值向量得到的中心化观测数据矩阵。
3. **中心化过程解释**：
- 对于每个观测 `i`，计算非目标观测的数量 `Nhat`。
- 计算 `Y_bar{i}`，即对应观测的非目标观测的均值向量。
- 计算 `Yc{i}`，即将原始观测数据矩阵 `Y` 减去对应观测的均值向量。
4. **输出解释**：
- `Yc`：中心化后的观测数据矩阵数组，每个元素对应一个观测的中心化数据。
- `Y_bar`：均值向量数组，每个元素对应一个观测的非目标观测的均值向量。
总体而言，这个函数实现了中心化数据的操作，以便在后续处理中更好地处理观测数据。
Constrained_ML_B.m
这是一个 MATLAB 函数，用于执行约束最大似然估计。下面是对该函数的解释：
```matlab
function [B,me]=constrained_ML_B(B,S,Y,K,sigma2,params)
% 检查是否提供了参数，如果没有，则创建一个空的参数结构
if(nargin==5)
params=[];
end
% 设置默认的最大迭代次数
if(isfield(params,'maxiter'))
maxiter=params.maxiter;
else
maxiter=1000;
end
% 获取输入矩阵的大小
[M,~]=size(Y);
% 创建单位矩阵
IM=speye(M);
% 对输入数据进行中心化处理
[Y,meanY]=centerY(Y,K);
% 初始化一些变量
ei=zeros(M,1);
s=sum(S,2);
I=speye(M);
% 迭代求解
for k=1:(maxiter)
Bold=B;
for i=1:M
% 获取未被目标的实验的索引
nontargs=K(i,:);
% 计算未被目标实验的数量
Nhat=sum(nontargs);
if (s(i))
ei(i)=1; % ei 是单位向量
% zi 用于计算 m_ij=cofactor_ij/det(I-B)
zi=(I-B)\ei;
js_i=find(S(i,:)>0);
for j=js_i
m_ij=zi(j);
B_old=B(i,j);
y_j=Y{i}(j,nontargs)';
% 移除 B(i,j) 并获取 B\ij 的第 i 行
BiT=B(i,:);
BiT(j)=0;
a_iT=(ei'-BiT)*Y{i}(:,nontargs);
r_ij=y_j'*y_j;
beta_ij=a_iT*y_j;
if (m_ij==0) % 转到线性方程
B(i,j)=beta_ij/r_ij;
else % 如果 m_ij ~= 0 转到二次方程
d_ij=1/m_ij+B(i,j);
theta_ij=r_ij*d_ij+beta_ij;
k_ij=d_ij*beta_ij-Nhat*sigma2;
q_ij=theta_ij^2-4*r_ij*k_ij;
Bijp=(1/(2*r_ij))*(theta_ij+sqrt(q_ij));
Bijm=(1/(2*r_ij))*(theta_ij-sqrt(q_ij));
candsBij=[Bijm, Bijp]';
Lss=sigma2*Nhat*log((d_ij-candsBij).^2)/2-0.5*sum((ones(length(candsBij),1)*a_iT-candsBij*Y{i}(j,nontargs)).^2,2);
[~,l_max]= max(Lss);
B(i,j)=candsBij(l_max);
end %m_ij
dB=B_old-B(i,j);
zi=zi/(1+dB*m_ij);
end%j
% 现在对扰动-基因连接矩阵 F 进行坐标上升
ei(i)=0;
end %s(i)
end %i
% 计算 B 的变化
delta_B=norm(Bold-B,'fro')/(norm(Bold,'fro')+1e-10);
% 如果 B 的变化小于 1e-2，则停止求解 B
if delta_B<1e-2
for i=1:M
me(i,1)=(IM(i,:)-B(i,:))*meanY{i}; % 计算与 B(i,:) 相关的 mue(i)
end
return
end
end %k
% 计算与 B(i,:) 相关的 mue(i)
for i=1:M
me(i,1)=(IM(i,:)-B(i,:))*meanY{i};
end
```
这个函数的目的是通过坐标上升法来执行约束最大似然估计，其中 `B` 是一个矩阵，`S` 是一个二值矩阵，`Y` 是一个细胞数组，`K` 是一个二值矩阵，`sigma2` 是方差参数，`params` 是可选参数。在迭代过程中，它尝试找到最大似然估计的解，同时满足一些约束条件。
Constrained_ridge_B.m
这是一个 MATLAB 函数，用于执行约束 Ridge 回归。下面是对该函数的解释：
```matlab
function [B,me]=constrained_ridge_B(Y,K,rho_factor)
% 获取输入数据的大小
[M,~]=size(Y);
% 将输入数据以中心值为0的方式进行中心化处理
[Y,meanY]=centerY(Y,K);
% 初始化 B 矩阵和单位矩阵
B=zeros(M);
IM=speye(M);
for i=1:M
% 如果是转录本，相关蛋白质不能进行调节
if i<=M/2
rows_rem=[i i+M/2];
I=speye(M-2);
else
rows_rem=i;
I=speye(M-1);
end
% 获取未被目标实验的索引
nontargs=K(i,:);
% 计算未被目标实验的数量
Nhat=sum(nontargs);
yi=Y{i}(i,nontargs)';
Yi=Y{i}(:,nontargs)';
% 移除相关的行
Yi(:,rows_rem)=[];
% 设置 Ridge 参数
Pi=eye(Nhat);
rho=rho_factor*norm(Yi'*Pi)^2;
% 计算 Ridge 回归系数
bi=(Yi'*Pi*Yi+rho*I)\(Yi'*(Pi*yi));
% 更新 B 矩阵
indexes=1:M;
indexes(rows_rem)=[];
B(i,indexes)=bi';
% 计算与 B(i,:) 相关的 mue(i)
me(i,1)=(IM(i,:)-B(i,:))*meanY{i};
end
```
这个函数的目的是通过 Ridge 回归执行约束最大似然估计，其中 `Y` 是一个细胞数组，`K` 是一个二值矩阵，`rho_factor` 是 Ridge 回归的参数。在迭代过程中，它尝试找到最大似然估计的解，同时满足一些约束条件。
Knockdown_visualizations.m
这个 MATLAB 脚本文件似乎用于可视化基因敲除实验结果。以下是一些关键部分的解释：
1. **输入参数:**
- `g`: 感兴趣的基因（基因名称）。
- `TargConst`: 感兴趣的构建（构建名称）。
- `Yact_table`: 包含实验观测数据的表格。
- `Xtable`: 与观测数据相关的其他信息表格。
- `Ypred`: 包含预测数据的结构体。
- `legend_flag`: 一个标志，表示是否显示图例。
2. **函数功能:**
- 该脚本函数的目的是生成多个图形，可视化基因敲除实验的结果。
- 第一个图形绘制了转录本的估计值，包括实验测量值、新模型估计值、旧模型（S1和S2）的估计值。
- 第二个图形绘制了蛋白质的估计值，也包括实验测量值、新模型估计值、旧模型的估计值。
- 第三个图形绘制了有针对性的基因敲除的条形图，显示了针对感兴趣的转录本和蛋白质的估计值。
3. **细节:**
- 图中的条形图显示了实验测量值、新模型估计值和旧模型估计值的比较。
- 有关转录本和蛋白质的测量值，图中还包括了相应的 WT (Wild Type) 条形，用虚线表示。
- 不同颜色的条形表示不同的数据源和估计模型。
- 图例（如果设置了 `legend_flag`）显示了不同颜色的含义。
总体而言，这个脚本通过条形图和线图可视化了基因敲除实验的结果，以便更好地理解不同模型和实验观测值之间的关系。
Model_predictions.m
这个 MATLAB 脚本文件似乎是用于生成基因敲除实验结果的预测。以下是一些关键部分的解释：
1. **版权声明:**
- 文件开头包含了版权声明，指出该程序是根据 GNU General Public License 的第二版发布的自由软件。它允许用户重新分发和/或修改软件，但在一定条件下，包括没有任何担保的情况下。
2. **函数功能:**
- 该脚本文件定义了一个函数 `Model_Predictions`，用于生成基因敲除实验的预测结果。
- 函数采用实验测量值 `Xact` 和针对性基因敲除的输入 `Xtarg`，然后使用已定义的模型参数进行预测。
- 如果提供了第三个输入 `Xfulltrans`，表示使用旧模型进行预测，并且该参数包含所有转录本的测量值。
3. **细节:**
- 函数首先加载了一些模型参数和实验数据。
- 针对旧模型（`prevhalf_flag` 为真），函数使用实验的全转录本测量值进行预测。
- 针对新模型，函数使用实验的部分转录本测量值，同时根据已知的蛋白质测量值进行调整。
- 函数返回一个结构体 `Ypred`，其中包含了对实验数据的三种预测：`prevhalf`、`prevfull` 和 `full`。
- 预测值通过矩阵求解和模型参数进行计算。
总体而言，这个脚本用于根据已知的实验测量值和模型参数生成基因敲除实验的预测结果。
TranscriptProteinModelPredictions.m
这段脚本文件用于模拟基因敲除实验，并通过可视化展示新模型和旧模型的预测结果。以下是一些关键部分的解释：
1. **版权声明:**
- 文件开头包含了版权声明，指出该程序是根据 GNU General Public License 的第二版发布的自由软件。
2. **清理工作区:**
- `clear` 和 `clc` 命令用于清除 MATLAB 工作区和命令窗口。
3. **加载实验数据和模型预测:**
- 使用 `load` 命令加载实验数据，包括转录本和蛋白质的丰度数据。
- 调用 `Model_Predictions` 函数，传入实验测量值和模型参数，生成基因敲除实验的预测结果。
4. **构建绘图数据:**
- 将实验数据进行进一步处理，生成用于可视化的数据，包括目标基因敲除的输入 `Xtarg` 和旧模型中使用的所有转录本的数据 `Xfulltrans`。
5. **可视化预测结果:**
- 使用 `Knockdown_Visualizations` 函数创建条形图，展示目标基因和目标基因敲除的实验数据、新模型预测和旧模型预测。
- 用户可以选择关注的基因 (`gene_of_interest`) 和目标基因敲除 (`knockdown_of_interest`)，并通过 `legend_flag` 控制是否显示图例。
总体而言，该脚本文件通过可视化展示新模型和旧模型在特定基因敲除实验条件下的预测结果，以便用户更好地理解模型性能。
扫一扫
609
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
