First version

glueTrees.m

@@ -0,0 +1,32 @@
+function fulltree=glueTrees(ttree,wtree)
+%Glue together the within-host trees using the transmission tree
+n=size(ttree,1);
+no=n+1;
+for i=n:-1:1%Consider all hosts
+    ni=size(wtree{i},1)/2;
+    for j=1:ni
+        lab(i,j)=no;no=no+1;
+    end
+    labh(i)=no-1;
+end
+leaves=[];
+intnodes=[];
+for i=1:n
+    tree=wtree{i};
+    ni=size(wtree{i},1)/2;
+    tree(:,1)=tree(:,1)+ttree(i,1);%Add infection time to all nodes
+    leaves=[leaves;tree(1,:)];
+    tree=tree(ni+1:end,:);
+    f=find(ttree(:,3)==i);%Infected by current nodes
+    for j=1:size(tree,1)
+        for k=2:3
+            if tree(j,k)==0,continue;end
+            if tree(j,k)==1,tree(j,k)=i;
+            elseif tree(j,k)<=ni,tree(j,k)=labh(f(tree(j,k)-1));
+            else tree(j,k)=lab(i,tree(j,k)-ni);end
+        end
+    end
+    intnodes=[tree;intnodes];
+end
+fulltree=[leaves;intnodes];
+fulltree=[fulltree hostFromFulltree(fulltree)];

hostFromFulltree.m

@@ -0,0 +1,26 @@
+function host=hostFromFulltree(fulltree)
+%Build vector 'host' indicating in which host each node of the fulltree is found
+fathers=zeros(size(fulltree,1)+1,1);
+fathers(fulltree(:,2)+1)=1:size(fulltree,1);
+fathers(fulltree(:,3)+1)=1:size(fulltree,1);
+fathers=fathers(2:end);
+host=zeros(size(fulltree,1),1);
+n=sum(fulltree(:,2)==0&fulltree(:,3)==0);
+for i=1:n
+    j=i;
+    while 1
+        host(j)=i;
+        j=fathers(j);
+        if fulltree(j,3)==0,break;end
+    end
+end
+f=n+find(fulltree(n+1:end-1,3)==0);%All transmission events
+for i=1:length(f)
+    j=f(i);
+    tocol=[];
+    while host(j)==0
+        tocol=[tocol j];
+        j=fathers(j);
+    end
+    host(tocol)=host(j);
+end

infer.m

@@ -0,0 +1,48 @@
+function record=infer(tree,mcmc,timeLastSeq)
+%Tree is a phylogeny
+%mcmc is the number of MCMC iterations to perform
+%timeLastSeq is the date of the last sequence
+if nargin<3,timeLastSeq=2010;end
+if nargin<2,mcmc=10000;end
+
+n=ceil(size(tree,1)/2);
+neg=100/365;%Within-host effective population size (Ne) times  generation duration (g)
+popsize=1000;%Size of population
+beta=0.001;%Infectiveness
+nu=1;%Rate of recovery
+
+%MCMC loop
+fulltree=makeFullTreeFromPTree(tree);%Starting point
+record(mcmc).tree=0;
+pTTree=probTTree(ttreeFromFullTree(fulltree),popsize,beta,nu);
+pPTree=probPTreeGivenTTree(fulltree,neg);
+for i=1:mcmc
+    %Record things
+    record(i).tree=absolute(fulltree,timeLastSeq);
+    record(i).pTTree=pTTree;
+    record(i).pPTree=pPTree;
+    record(i).neg=neg;
+    record(i).source=fulltree(end-1,4);
+    record(i).nu=nu;
+    record(i).beta=beta;
+    %Metropolis update for transmission tree
+    fulltree2=proposal(fulltree);
+    pTTree2=probTTree(ttreeFromFullTree(fulltree2),popsize,beta,nu);
+    pPTree2=probPTreeGivenTTree(fulltree2,neg);
+    if log(rand)<pTTree2+pPTree2-pTTree-pPTree,fulltree=fulltree2;pTTree=pTTree2;pPTree=pPTree2;end
+    %Metropolis update for Ne*g, assuming Exp(1) prior
+    neg2=neg+(rand-0.5)*record(1).neg;
+    if neg2<0,neg2=-neg2;end
+    pPTree2=probPTreeGivenTTree(fulltree,neg2);
+    if log(rand)<pPTree2-pPTree-neg2+neg,neg=neg2;pPTree=pPTree2;end
+    %Gibbs updates for beta and nu based on O'Neill and Roberts JRSSA 16:121-129 (1999) and assuming Exp(1) priors
+    ttree=ttreeFromFullTree(fulltree);
+    intXtYtdt=intXtYtdtGivenTTree(ttree,popsize);
+    beta=gamrnd(1+n-1,1/(1+intXtYtdt));
+    nu=gamrnd(1+n,1/(1+n*mean(ttree(:,2)-ttree(:,1))));
+end
+
+function ft=absolute(fulltree,timeLastRem)
+%Convert from relative to absolute time using timeLastRem which is the date of the last removal (ie the last sequence)
+ft=fulltree;
+ft(:,1)=ft(:,1)+timeLastRem-max(ft(:,1));

intXtYtdtGivenTTree.m

@@ -0,0 +1,17 @@
+function intXtYtdt=intXtYtdtGivenTTree(ttree,popsize)
+%return the integral of Xt*Yt*dt
+
+[times,ind]=sort([ttree(:,1);ttree(:,2)]);%Event times
+types=ind>size(ttree,1);%0 for infection, 1 for recovery
+sir=[popsize-1 1 0];
+intXtYtdt=0;
+for i=2:length(times)
+    intXtYtdt=intXtYtdt+sir(1)*sir(2)*(times(i)-times(i-1));
+    if types(i)==1
+        %Recovery
+        sir(3)=sir(3)+1;sir(2)=sir(2)-1;
+    else
+        %Infection
+        sir(1)=sir(1)-1;sir(2)=sir(2)+1;
+    end
+end

makeFullTreeFromPTree.m

@@ -0,0 +1,27 @@
+function fulltree=makeFullTreeFromPTree(tree)
+%Return a non-zero probability phylogenetic+transmission tree for a given phylogenetic tree
+n=ceil(size(tree,1)/2);
+tree=[tree;zeros(n,3)];
+tree(end,1)=min(tree(1:(2*n-1),1))-1;
+tree(end,2)=2*n-1;
+[~,source]=max(tree(1:n,1));%The source is the last leaf
+notsource=setdiff(1:n,source);
+i2=0;
+for i=notsource
+    i2=i2+1;
+    f=find(tree(:,2)==i|tree(:,3)==i);
+    if tree(f,2)==i,tree(f,2)=2*n-1+i2;else tree(f,3)=2*n-1+i2;end
+    tree(2*n-1+i2,2)=i;
+    tree(2*n-1+i2,1)=(tree(f,1)+tree(i,1))/2;
+end
+tree(:,1)=tree(:,1)-min(tree(:,1));
+
+%Reorder nodes chronologically
+[~,ind]=sort(tree(n+1:end,1),'descend');
+for i=n+1:size(tree,1)
+    for j=2:3
+        if tree(i,j)>n,tree(i,j)=n+find(ind==tree(i,j)-n);end
+    end
+end
+tree=tree([(1:n)';n+ind],:);
+fulltree=[tree hostFromFulltree(tree)];

makeTTree.m

@@ -0,0 +1,35 @@
+function [ttree,prob]=makeTTree(popsize,beta,nu)
+%Creates a transmission tree, and returns the result as a N*3 matrix in the
+%following format:
+%One row per infected host
+%First column is time of infection
+%Second column is time of recovery
+%Third column is infector
+%nu is the rate of recovery
+%beta is the rate of infection
+%popsize is the total size of population
+ttree=[0 0 0];
+sir=[popsize-1 1 0];
+curt=0;
+prob=0;
+while 1
+    rate=nu*sir(2)+beta*sir(1)*sir(2);
+    if rate==0,break;end
+    t=exprnd(1/rate);
+    prob=prob+log(exppdf(t,1/rate));
+    curt=curt+t;
+    w=find(ttree(:,2)==0);%All infected
+    prob=prob+log(1/length(w));
+    w=w(1+floor(length(w)*rand));%Pick one at random
+    if rand<nu*sir(2)/rate
+        %Recovery
+        prob=prob+log(nu*sir(2)/rate);
+        ttree(w,2)=curt;
+        sir(3)=sir(3)+1;sir(2)=sir(2)-1;
+    else
+        %Infection
+        prob=prob+log(1-nu*sir(2)/rate);
+        ttree=[ttree;curt 0 w];
+        sir(1)=sir(1)-1;sir(2)=sir(2)+1;
+    end
+end

plotBothTree.m

@@ -0,0 +1,96 @@
+function plotBothTree(tree,method,names)
+%Plot both phylogenetic and transmission trees on the same figure
+%The second parameter indicates whether to show stars for transmission
+%events (method=0), color branches according to host (method=1) or use
+%boxes to represent hosts (method=2)
+n=sum(tree(:,2)+tree(:,3)==0);
+if nargin<3,for i=1:n,names{i}=sprintf('%d',i);end;end
+if nargin<2,method=1;end
+hold on
+host=tree(:,4);
+cmap=colormap;
+
+%Assign a unique color to each host
+cols=zeros(n,3);
+for i=1:n
+    cols(i,:)=cmap(min(size(cmap,1),1+floor(size(cmap,1)*i/n)),:);
+end
+
+%Assign alternating colors to hosts
+%basiccols=cmap(1:round(size(cmap,1)/6):size(cmap,1),:);
+%basiccols=[1 0 0;0 1 0;0 0 1];
+%co=zeros(size(tree,1),1);
+%td=size(tree,1);
+%co(td)=0;
+%while ~isempty(td)
+%    chi=setdiff(tree(td(1),2:3),0);
+%    if length(chi)==1,co(chi)=1+mod(co(td(1)),size(basiccols,1));else co(chi)=co(td(1));end
+%    td=[td(2:end) chi];
+%end
+%cols=basiccols(co(1:n),:);
+
+%Determine ys
+ys=zeros(n,1);
+todo=[size(tree,1) 0 0.5 1];%Matrix of nodes to do, with associated starting x and y coordinates and scale
+while size(todo,1)>0
+    w=todo(1,1);x=todo(1,2);y=todo(1,3);scale=todo(1,4);
+    if tree(w,2)==0 && tree(w,3)==0
+        %Leaf node
+        ys(w)=y;
+    elseif tree(w,3)==0
+        %Transmission node
+        todo=[todo;tree(w,2) tree(w,1) y scale];
+    else
+        %Binary node
+        todo=[todo;tree(w,2) tree(w,1) y+scale/2 scale/2;tree(w,3) tree(w,1) y-scale/2 scale/2];
+    end
+    todo=todo(2:end,:);%Remove first node
+end
+[~,ys]=sort(ys);
+ys(ys)=1:length(ys);
+ylims=[ys ys];
+for i=n+1:size(tree,1)
+    f=1+find(tree(i,2:3)>0);
+    ylims(i,1)=min(ylims(tree(i,f),1));
+    ylims(i,2)=max(ylims(tree(i,f),2));
+    ys(i)=mean(ys(tree(i,f)));
+end
+
+todo=[size(tree,1) tree(end,1)];%Matrix of nodes to do, with associated starting x
+while size(todo,1)>0
+    w=todo(1,1);x=todo(1,2);y=ys(w);
+    if method==1&&host(w)>0,col=cols(host(w),:);else
+    if tree(w,2)>0,col=cols(host(tree(w,2)),:);else col=[0 0 0];end;end
+    if tree(w,2)==0 && tree(w,3)==0
+        %Leaf node
+        if method==2
+           patch([tree(w,1) tree(w,1) max(tree(:,1)) max(tree(:,1))],[ylims(w,1)-0.5 ylims(w,2)+0.5 ylims(w,2)+0.5 ylims(w,1)-0.5],[1 1 1],'LineStyle','none');
+        end
+        p=plot([x tree(w,1)],[y y],'Color',col*(method==1),'LineWidth',2);
+        uistack(p,'bottom');
+        text(tree(w,1)+(max(tree(:,1))-min(tree(:,1)))/100,y,names{w},'FontSize',12);
+    elseif tree(w,3)==0
+        %Transmission node
+        p=plot([x tree(w,1)],[y y],'Color',col*(method==1),'LineWidth',2);
+        uistack(p,'bottom');
+        if method==0||method==1
+            plot(tree(w,1),y,'r*','MarkerSize',10);
+        elseif method==2
+            patch([tree(w,1) tree(w,1) max(tree(:,1)) max(tree(:,1))],[ylims(w,1)-0.5 ylims(w,2)+0.5 ylims(w,2)+0.5 ylims(w,1)-0.5],col,'LineStyle','none');
+%            patch([tree(fi,1) tree(fi,1) max(tree(:,1)) max(tree(:,1))],[y-scale y+scale y+scale y-scale],[1 1 1],'LineStyle','none');
+        end
+        todo=[todo;tree(w,2) tree(w,1)];
+    else
+        %Binary node
+        p=plot([x tree(w,1)],[y y],'Color',col*(method==1),'LineWidth',2);
+        uistack(p,'bottom');
+        p=plot([tree(w,1) tree(w,1)],[ys(tree(w,2)) ys(tree(w,3))],'Color',col*(method==1),'LineWidth',2);
+        uistack(p,'bottom');
+        todo=[todo;tree(w,2) tree(w,1);tree(w,3) tree(w,1)];
+    end
+    todo=todo(2:end,:);%Remove first node
+end
+set(gcf,'Color','w');
+set(gca,'YTick',[]);
+set(gca,'YColor','w');
+ylim([0 n+1]);

plotPTree.m

@@ -0,0 +1,11 @@
+function plotPTree(tree)
+%This function takes a phylogenetic tree as input
+
+n=ceil(size(tree,1)/2);
+br=tree(n+1:end,2:3);
+for i=1:n*2-2
+    f=find(tree(:,2)==i|tree(:,3)==i);
+    d(i)=tree(i,1)-tree(f,1);
+end
+d=[d 0];
+plot(phytree(br,d'));

plotTTree.m

@@ -0,0 +1,13 @@
+function plotTTree(ttree)
+%Plots a transmission tree
+n=size(ttree,1);
+
+cm=zeros(n+1,n+1);
+ids{1}='SOURCE';
+for i=1:n
+    cm(ttree(i,3)+1,i+1)=1;
+    ids{i+1}=sprintf('%d [%.2f;%.2f]',i,ttree(i,1),ttree(i,2));
+end
+bg=biograph(cm,ids,'LayoutScale',0.5);%Create biograph
+set(bg.nodes(:),'FontSize',12);
+g=biograph.bggui(bg);%View biograph