文献—Emergent simplicity in microbial community assembly——中使用的交叉互养模型的代码分析
2024-06-02 11:22:52
本文对Emergent simplicity in microbial community assembly——中使用的交叉互养模型的代码分析
原始文献Goldford J E , Lu N , Bajic D , et al. Emergent Simplicity in Microbial Community Assembly[J]. Science, 2018, 361(6401):469-474.
原始代码链接https://github.com/jgoldford/mcrm.
文章说明和系列内容
本系列更新主要是关于2018年一篇群落结构组成的文献进行阅读的分析和总结，将全文的过程进行分解。主要分为下面几个篇章。本系列的全部内容都为原创，由于个人水平有限，整理过程可能对一些理解不到位甚至错误，请辩证的阅读。
本文字数共（11703）大于阅读5分钟
- 论文全过程详细阅读整理
- 全文分析
- 文献模型过程推导
- 代码整理过程注释
本文代码运行说明
matlab2018b（原始代码环境matlab2015a）本文对原始代码进行整理，加入了一下自己的思考。适合初学者了解文献中实现模型的过程。对原始文献部分地方进行了调整（具体见代码部分）
代码过程
下图是对文献中的代码流程的一个汇总绘制，其中关于涉及的具体公式参考文献内部和本系列数学推导过程
对文献中实现过程的思考
- 持续更新 2020年11月25日
文献在对群落中代谢过程的分析可以从采样模型进行改进，采取傅里叶变化的方式对代谢规律进行分析，观察群落中不同代谢物质在群落持续变化对资源的利用效果。分析在不同科水平中加入一下生化相关的信息。
代码结果展示
文献具体图在上一篇文章有上传
结果输出：
代码过程
主函数部分
%% 时间 2020年11月25日08:55:26
%% 版本 matlab 2018b
%% 注释 尚帝
%% 原始代码地址 https://github.com/jgoldford/mcrm.
%% 整理邮箱: 尚帝@163.com
%% %%
%% 设置群落的参数属性
% define hyperparamters
% the total number of species in your "world"Total_Number_Of_Species = 1000;% the number of species in your ecosystem (i.e. inocula)
numSpecies = 100;% the total number of metabolic families
numGroups = 5;% the total number of resources (must equal or exceed numGroups)
numResource = 10; %资源总数
% x 获取资源数目
% dirichlet hypoerparameter (controls the variablity of the family
% metabolism)
Total = 100;
%% 生成分泌代谢矩阵
% get a secretion matrix for everyone %分泌矩阵
D = getMetabolism(numResource,'rand',1/numResource);% determine how much of the community you want to be specialists vs.
% generalists
% more generalist, set to 1/M, and more specialist set to 0.9
specialist = 0.5;specialistVariation = 0.01;
T=10; %外源物质转化效率/摄取资源量
% 每个科水平中 每个物质的消耗水平
priors = arrayfun(@(x) getConsumerPriors(T,x,specialist,specialistVariation,numResource)',1:numGroups,'uni',0);
%每个菌属的代谢矩阵
[out,~]=makePhyloConsumers(Total_Number_Of_Species,numResource,numGroups,priors,Total,true);% if you don't want community struture, use the function
% getRandomConsumerMatrix% randomly sample a sub population:
k = randsample(1:Total_Number_Of_Species,numSpecies);
% construct an ecosystem
params = mcrm_params(1,numSpecies,out.C(k,:),D,'eye','');% run ODE solver for ecosystem
r = run_mcrm(params);
figure();
subplot(1,3,1);
plot(r.species);
ylabel('species abundance');g = out.group(k((r.communityStruct > 1e-4)));
% x = r.communityStruct(r.communityStruct > 1e-4);
%
% subplot(1,3,2);
% pie(x)
%% 探究群落中丰富度对稳定性的影响，分析群落中物种的丢失率
% compute the coarse grained structure of the community
x = coarseGrainCommunityStructure(r.communityStruct,out.group(k),numGroups);
%figure();
subplot(1,3,2)
x(x<1e-10) = 0;
pie(x+1e-20)
legend(arrayfun(@(x) ['Family ',num2str(x)],1:numGroups,'uni',0));% how stable is the community to species loss?
% this function reduces the parameter structure to only include species
% that survived in the prior simulation
p_new = removeExtinctSpecies(params,r,1e-4);% knockout a species (i.e. set it's consumption rates to zero)I = zeros(p_new.num_species);for ko_species = 1:p_new.num_speciesparams_knockout = knockoutSpecies(p_new,ko_species);ko_sim = run_mcrm(params_knockout);% compute the number of species that went co-extinctextinct = find(ko_sim.communityStruct < 1e-4);%寻找物种丰富度低的I(ko_species,extinct) = 1;%extinct = setdiff(extinct,ko_species);end%compute average species loss
avgLoss = (sum(sum(I)) - length(I))/length(I);
fractionalLoss = avgLoss/p_new.num_species;
disp(strcat('平均损失的物种数目',num2str(fractionalLoss)))
%基于傅里叶变化分析代谢的种类差异 非原始代码
subplot(1,3,3)
Y = fft(r.resources,[],2);
P2 = abs(Y/(size(Y,1)*size(Y,2)));
P1 = P2(1:size(r.resources,2)); %以样本长度读取代谢信号
P1(2:end-1) = 2*P1(2:end-1); %去两倍信号延拓
plot([1:length(P1)],P1)
title('代谢的能力的傅里叶变化')
xlabel('10种资源')
ylabel('资源代谢频率分布')
% surf(r.species)
功能实现函数(非运行顺序)
function [coarsedGrainedAbundanceVec] = coarseGrainCommunityStructure(abundanceVector,groupVector,maxGroups)
% function sums the total abundance within each group (groupVector conains the group id for each species)
coarsedGrainedAbundanceVec = arrayfun(@(y) sum(abundanceVector(groupVector==y)),1:maxGroups);
end
function r = drchrnd(a,n)
p = length(a);
r = gamrnd(repmat(a,n,1),1,n,p);
r = r ./ repmat(sum(r,2),1,p);
end
微生物群落随机部分更新
function [params_out] = evolve_member(params,percentChange)
%EVOLVE_MEMBER takes a parameter structure, copies a random member and
%"mutates" the consumer vector by a max. percent change. This new member
%is appended into the mcrm structure% Detailed explanation goes here
species = randsample(1:params.num_species,1);
total = sum(params.C(species,:));
uptake = params.C(species,:)./total;
dir = binornd(1,[0.5]);
resource = randsample(1:params.num_resources,1);if dir > 0uptake(resource) = uptake(resource) + percentChange;
elseuptake(resource) = uptake(resource) - percentChange;
end
uptake(uptake<0) = 0;
uptake = uptake./sum(uptake);% now there is a small chance of changing the total enzyme capacity through
% some global mutation
p_change = 0;
q = binornd(1,p_change);
if q == 1total = total + total.*normrnd(0,0.01);
endc = uptake.*total;d = 1e-2.*params.x0(params.varIdx.species(species));% append a new member to the community structure
params_out = appendMember(params,c,d,0);end设定群落中微生物的能量代谢古城，相关代谢矩阵的属性
function [R] = getConsumerPriors(T,r,specialist,specialistVar,numResources)
% function generates a single consumer vector R, that represents the total
% uptake capacity for a group of organisms% the specialist paramaters controls the propotion of total uptake capacity
% (T) that is dedicated to uptake resource "r"% input
% T = total uptake capacity (e.g. controlled by fixed amount of enzymes, transporters)
% specialist, specialistVar = mean and std to randomly draw the
% fraction of uptake capcity dedicated to consumer resource "r"
% r = index of resource for which a consumer will be a "specialist" on
% numResources = total number of resources % initialize resource consumer vector
R = zeros(numResources,1);
% pick a dominate resource
if isempty(r)r = randsample(1:numResources,1);
end
%draw proportion of fixed resources 固定资源获取比例
f = normrnd(specialist,specialistVar); %正态分布的一随机数，（平均获取能力，均差）
% ensure f falls between 0 and 1
if f<0f = 0;
elseif f>1f = 1;
end% choose other resource consumption proportions from a uniform
% distribution
R = rand(numResources,1);
% (T-f*T) 剩余资源 除了自身以外的消耗速率百分比R(setdiff(1:numResources,r))./sum(R(setdiff(1:numResources,r)))'
R(setdiff(1:numResources,r)) = (T-f*T).*R(setdiff(1:numResources,r))./sum(R(setdiff(1:numResources,r)))';% compute fraction of total uptake capacity dedicated to consuming resource
% "r"
R(r) = f*T;end
初始物种丰富度
function [params] = knockoutSpecies(params,speciesIdx)
% Set species abundance to zero
%
params.C(speciesIdx,:) = zeros(length(speciesIdx),params.num_resources);end
按科水平进行对微生物的代谢资源种类
function [out,priors] = makePhyloConsumers(num_species,num_resources,num_groups,priors,Total,nonNormal)
% function to make sets of consumer matricies drawn from priors,
% representing "families" of consumers.% if no prior family-level consumer vectors, make some
if isempty(priors)
priors = arrayfun(@(x) Total.*rand(1,num_resources),1:num_groups,'uni',0);
else
end% draw consumption vectors from dirichlet distribution
% 每次运行返回一个群落数目{（物种数目×资源数目）}的元细胞
C = cellfun(@(y) drchrnd(y,num_species), priors,'uni',0);% assign wich consumers are part of which family
group_labels = arrayfun(@(y) y.*ones(1,num_species), 1:num_groups,'uni',0);
group_labels = cell2mat(group_labels);C_total = cell2mat(C');
% randomly sample a sub-population
idx = randsample(1:size(C_total,1),num_species);
C = C_total(idx,:);
group_labels = group_labels(idx);% breaks the constraint that there has to be a "fixed enzyme budget per
% species"if nonNormal %运行有随机的浮动代谢% allow for a ~2 percent change in the total uptake capacityu = normrnd(1,0.01,[num_species,1]);C = diag(u)*C;% for i = 1:max(group_labels)
% u = normrnd(1,0.01);
% C(group_labels==i,:) = C(group_labels==i,:).*u;
% uptakeFactor(i) = u;
% end% out.uptakeFactor = uptakeFactor;end% construct output structure
out.C = C;
out.group = group_labels;
[~,out.k]=sort(out.group);end
物种随机抽样
function [C] = getRandomConsumerMatrix(numSpecies,numResources,ctype)
% function to greate random consumer matricies
% if ctype is normalized, then make sure sum_alpha c_i,alpha = 1if ~isfield(ctype,'normalized')normalizedBool = true;
elsenormalizedBool = false;
end% randomly sample resource utilization matrix (C) from uniform distribution
C = rand(numSpecies,numResources);
% normalize distribution for each species
if normalizedBoolC = cell2mat(cellfun(@(x) x./sum(x), num2cell(C,2),'uni',0));
endend
设定代谢物质矩阵的生成方式（文献中为随机）
function [ D ] = getMetabolism(numResources,flag,p)
%GETMETABOLISM Function used to generate random stoichiometric matriciesswitch flagcase 'thermo'w = [numResources:-1:1];% construct D matrixd = fliplr(w)./numResources;D = cell2mat(arrayfun(@(x) [zeros(1,x),d(1:end-x)],1:numResources,'uni',0)');D = D';% rescale D s.t. W > D'WD = 0.3*D;case 'sparse'w = [numResources:-1:1];% construct D matrixd = fliplr(w)./numResources;D = cell2mat(arrayfun(@(x) [zeros(1,x),d(1:end-x)],1:numResources,'uni',0)');D = D';% rescale D s.t. W > D'WD = 0.3*D;M = binornd(1,p,[numResources,numResources]);D = M.*D;case 'rand'D = rand(numResources);scale = p;D = D.*scale;case'zeros'D = ones(numResources);
endend研究丰富度对群落的影响，去除不能生存的微生物（代谢的物质不足，或者丰富度低于预设值）
function [ params ] = removeExtinctSpecies(params,simulation,minAbundance)
%REVMOVEEXTINCTSPECIES removes species at the end of the simulation that
%are zero
% x = simulation.species(end,:);
x(x<0) = 0;
k = x > minAbundance;
params.num_species = sum(k);
params.varIdx.species = 1:(params.num_species);
params.varIdx.resources = (params.num_species+1):(params.num_resources + params.num_species);
x0 = [x(k)';simulation.resources(end,:)'];
params.x0 = x0;
params.C = params.C(k,:);
params.death_rate = params.death_rate(k);end资源竞争模型计算模拟
%% 资源竞争模型计算模拟
function [output] = run_mcrm(params)% parse input
num_species = params.num_species;
num_resources = params.num_resources;
varIdx = params.varIdx;
C = params.C;
D = params.D;
%Q = params.Q;
W = params.W;
B = params.B;
death_rate = params.death_rate;
mu = params.mu;
T = params.T;
tau = params.tau;
x0 = params.x0;
timeSteps = params.timeStep;
alpha = params.alpha;% create a function handle for dynamics
y = @(t,x) population_dynamics(t,x,num_species,C,D,W,T,mu,alpha,death_rate,B,tau,varIdx);% solve ODE
[T,Y] = ode15s(y,1:timeSteps,x0);% parse species, resource abundanes and consumer_matrix into output structure
output.species = Y(:,varIdx.species);
output.resources = Y(:,varIdx.resources);
%output.consumer_matrix = C;
%output.weights = W;
%output.production_matrix = D;
output.communityStruct = output.species(end,:);
output.environmentalStruct = output.resources(end,:);end
%%求解微分方程
function [dx] = population_dynamics(t,x,num_species,C,D,W,T,mu,alpha,death_rate,B,tau,varIdx)dx = zeros(sum(length(varIdx.resources)+num_species),1);% uptake, consumption matrix%C = update_consumption_matrix(C,x(varIdx.resources),diag(W));% calculate growth rate multiplier ([sum_j w_j * c_ij * r_j - T])growth_rate_multiplier = C*W*x(varIdx.resources) - T;% calculate consumption rate of resourcesconsumption = (C*diag(x(varIdx.resources)))'*x(varIdx.species);% compute production rate production = D*consumption;% calculate differentialdx(varIdx.species) = x(varIdx.species)*diag(mu).*growth_rate_multiplier - death_rate.*(x(varIdx.species));dx(varIdx.resources) = (alpha-x(varIdx.resources))./tau - consumption + production + B.*sum(death_rate.*x(varIdx.species));%dx = dx';end群落属性更新（忘了具体干啥的了。。。。）
function [params] = appendMember(params,c,abundance,death_rate)
%APPENDMEMBER takes a parmater structure and appends a new member (consumer vector) onto this structure
% inputs:
% params : mcrm paramater structure
% c : consumer vector (1XM dim)
% abundance : initial abundance of this consumer
% dealth rate: rate that this species dies relative to per capita growth (usually set to 0)params.x0 = [params.x0(params.varIdx.species);abundance;params.x0(params.varIdx.resources)];
params.varIdx.species = [params.varIdx.species,params.varIdx.species(end)+1];
N = length(params.varIdx.species);
params.num_species = N;
params.varIdx.resources = (N+1):N+length(params.varIdx.resources);params.C = [params.C;c];
params.death_rate = [params.death_rate;death_rate];end
总结
2020年11月25日
。。。。下次再说
文献—Emergent simplicity in microbial community assembly——中使用的交叉互养模型的代码分析相关推荐
- 文献—Emergent simplicity in microbial community assembly--论文全过程详细阅读整理与翻译
本文对Emergent simplicity in microbial community assembly-- 论文全过程详细阅读整理 原始文献Goldford J E , Lu N , Bajic ...
- 文献—Emergent simplicity in microbial community assembly-- 全文框架分析
本文对Emergent simplicity in microbial community assembly-- 全文框架分析 原始文献Goldford J E , Lu N , Bajic D , ...
- 目标检测，FFmpeg中第一个基于深度学习模型的视频分析功能
2021年4月,终于把目标检测(object detection)加到FFmpeg upstream了,有maintainer身份加持,还是交互了将近100封邮件,花了两个多月才完成upstream, ...
- 微生物恒化器中的进化压力----进化模型推导与分析
微生物恒化器中的进化压力 博文导读 内容说明 模型背景与假设 模型的建立 关键模型公式 完整推导过程 文献内容 基础代谢模型 文献延申 博文导读 微生物在进化过程中有多样性的代谢策略,通常我们观察微生 ...
- 文献阅读：神经网络提取生物医学文本中的关系
文献阅读:神经网络提取生物医学文本中的关系 题目 1 背景 2 相关工作 3 方法 3.1 生物医学关系提取 3.2 Dependency graphs and SDPs 3.3 句子嵌入表示 3.4 ...
- 爬取知网博硕士文献及中国专利存到mysql数据库中的代码及其注意事项
今天因为需要做了一个爬取知网博硕士论文及中国专利的爬虫,在制作的过程中遇到了不少坑,在网上查资料时都是很老的资源,在现在知网的反爬虫下不起作用,所以我来写这篇文章来供大家参考.(这篇文章主要介绍通过改 ...
- 文献记录(part79)--光学影像序列中基于多视角聚类的群组行为分析
学习笔记,仅供参考,有错必纠 光学影像序列中基于多视角聚类的群组行为分析 摘要 群组行为分析是光学影像序列分析中的一项重要课题, 在近年来引起了人工智能领域研究人员的广泛关注. 与行人个体相比, 群组 ...
- c# typescript_在任何IDE中从C＃，Java或Python代码获取TypeScript接口的简单方法
c# typescript by Leonardo Carreiro 莱昂纳多·卡雷罗(Leonardo Carreiro) 在任何IDE中从C#,Java或Python代码获取TypeScript接 ...
- NC：自体免疫水泡皮肤病中鉴定基因与微生物组互作(微生物组关联分析MWAS)
之前我们平台分享过<Nature:宏基因组关联分析>综述,让大家系统的了解这一领域.同时还分享过一篇 <GigaScience:谷子产量与微生物组关联分析>是植物领域中的优秀工 ...
最新文章
- 枚举保存到数据库中为序号解决方案
- NET130署名错误一事，改正也着实迅速
- ASP.NET企业开发框架IsLine FrameWork系列之七--AppLogProvider日志框架(上)
- 云计算教程学习入门视频课件：云计算架构参考模型
- 计算机云客户端技术指标,云服务器技术指标
- 如何正确使用广告素材、优化Facebook广告
- 北京大学生物信息学(转录组)
- jQuery - 选择器（五）
- iOS应用程序审核机制
- 用ie浏览器签章后保存在桌面显示不出文件
- sam卡和sim卡区别_PSAM卡、SAM卡与SIM卡
- 一个微信关联管理多个腾讯云账号
- Python——第一天的Suger Rush
- 【JZOJ4847】【NOIP2016提高A组集训第5场11.2】夕阳
- Android开发人才前景分析及建议
- 处理反走样常用的四种技术
- ShuZu冒泡排序选择排序
- 利普希茨【NOIP2017模拟8.7A组】
- 【转】iOS右滑返回手势全解和最佳实施方案
- latex行内公式和行间公式