From 03409db31c3b03117cfc06ad76364c9601af1ff2 Mon Sep 17 00:00:00 2001
From: "xavier.didelot" <xavier.didelot@gmail.com>
Date: Tue, 25 Mar 2014 07:59:35 +0000
Subject: [PATCH] First version

---
 glueTrees.m             | 32 ++++++++++++++
 hostFromFulltree.m      | 26 +++++++++++
 infer.m                 | 48 +++++++++++++++++++++
 intXtYtdtGivenTTree.m   | 17 ++++++++
 makeFullTreeFromPTree.m | 27 ++++++++++++
 makeTTree.m             | 35 +++++++++++++++
 plotBothTree.m          | 96 +++++++++++++++++++++++++++++++++++++++++
 plotPTree.m             | 11 +++++
 plotTTree.m             | 13 ++++++
 probPTreeGivenTTree.m   | 87 +++++++++++++++++++++++++++++++++++++
 probTTree.m             | 27 ++++++++++++
 proposal.m              | 61 ++++++++++++++++++++++++++
 ptreeFromFullTree.m     | 11 +++++
 readPTree.m             | 17 ++++++++
 simu.m                  | 31 +++++++++++++
 ttreeFromFullTree.m     | 14 ++++++
 tutorial.m              | 31 +++++++++++++
 withinhost.m            | 44 +++++++++++++++++++
 18 files changed, 628 insertions(+)
 create mode 100644 glueTrees.m
 create mode 100644 hostFromFulltree.m
 create mode 100644 infer.m
 create mode 100644 intXtYtdtGivenTTree.m
 create mode 100644 makeFullTreeFromPTree.m
 create mode 100644 makeTTree.m
 create mode 100644 plotBothTree.m
 create mode 100644 plotPTree.m
 create mode 100644 plotTTree.m
 create mode 100644 probPTreeGivenTTree.m
 create mode 100644 probTTree.m
 create mode 100644 proposal.m
 create mode 100644 ptreeFromFullTree.m
 create mode 100644 readPTree.m
 create mode 100644 simu.m
 create mode 100644 ttreeFromFullTree.m
 create mode 100644 tutorial.m
 create mode 100644 withinhost.m

diff --git a/glueTrees.m b/glueTrees.m
new file mode 100644
index 0000000..ee256c8
--- /dev/null
+++ b/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)];
diff --git a/hostFromFulltree.m b/hostFromFulltree.m
new file mode 100644
index 0000000..bdb7ef2
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/infer.m b/infer.m
new file mode 100644
index 0000000..f932e35
--- /dev/null
+++ b/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));
\ No newline at end of file
diff --git a/intXtYtdtGivenTTree.m b/intXtYtdtGivenTTree.m
new file mode 100644
index 0000000..52a3ffb
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/makeFullTreeFromPTree.m b/makeFullTreeFromPTree.m
new file mode 100644
index 0000000..a595cbd
--- /dev/null
+++ b/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)];
\ No newline at end of file
diff --git a/makeTTree.m b/makeTTree.m
new file mode 100644
index 0000000..dc0ccd3
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/plotBothTree.m b/plotBothTree.m
new file mode 100644
index 0000000..7dc467e
--- /dev/null
+++ b/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]);
\ No newline at end of file
diff --git a/plotPTree.m b/plotPTree.m
new file mode 100644
index 0000000..796e0c6
--- /dev/null
+++ b/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'));
\ No newline at end of file
diff --git a/plotTTree.m b/plotTTree.m
new file mode 100644
index 0000000..51e4468
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/probPTreeGivenTTree.m b/probPTreeGivenTTree.m
new file mode 100644
index 0000000..f63464d
--- /dev/null
+++ b/probPTreeGivenTTree.m
@@ -0,0 +1,87 @@
+function prob=probPTreeGivenTTree(fulltree,neg)
+%To calculate the probability of a phylogenetic tree given a transmission
+%tree, split up the phylogenetic tree into individual within-host subtrees, and
+%calculate the product of their probabilities
+n=sum(fulltree(:,2)==0&fulltree(:,3)==0);
+prob=0;
+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);
+for i=1:n
+    subtree=extractSubtree(fulltree,i,fathers);
+    prob=prob+probSubtree(subtree,neg);
+end
+
+function subtree=extractSubtree(fulltree,which,fathers)
+%Take all nodes in host
+host=fulltree(:,4);
+ind=1:size(fulltree,1);
+ind=ind(host==which);
+%Add father of oldest node
+[~,j]=min(fulltree(ind));
+ind=[ind fathers(ind(j))];
+%Create subtree
+subtree=zeros(length(ind),2);
+invind=zeros(size(fulltree,1));invind(ind)=1:length(ind);
+for i=1:length(ind)
+    subtree(i,1)=fulltree(ind(i),1);
+    if i<length(ind)
+        subtree(i,2)=invind(fathers(ind(i)));
+    end
+end
+
+function p=probSubtree(tab,rate)
+%tab(:,1)=age at bottom;tab(:,2)=father;rate=coalescence rate
+%Return the log-prior probability of a subtree
+%This is an extension to Eq1 of Drummond et al (2002) Genetics
+%161:1307-1320 that accounts for condition TMRCA<INCUBATION_PERIOD
+p=0;
+tab(:,1)=max(tab(:,1))-tab(:,1);%convert times to ages
+isiso=ones(size(tab,1),1);
+isiso(tab(1:end-1,2))=0;
+iso=find(isiso);
+ex=zeros(size(tab,1),1);%Keep track of which nodes are active
+[s,ind]=sort(tab(:,1));
+cur=iso(1);while tab(cur,2)>0,ex(cur)=1;cur=tab(cur,2);end;ex(end)=1;%Activate all ancestors of first leave
+for l=2:length(iso)%For all others leaves in increasing order of age
+    isanc=zeros(size(tab,1),1);%Ancestors of the current node
+    anc=iso(l);while tab(anc,2)>0,isanc(anc)=1;anc=tab(anc,2);end
+    bra1=0;
+    bra2=0;
+    start=false;
+    found=false;
+    curage=0;
+    k=0;
+    for i=1:length(s)
+        if ind(i)==iso(l),start=true;curage=tab(iso(l),1);end
+        if ex(ind(i))==0,continue;end%Ignore non-existant nodes
+        if start
+            if found==false
+                bra1=bra1+k*(tab(ind(i),1)-curage);
+                if isanc(ind(i)),found=true;end
+            end
+            bra2=bra2+k*(tab(ind(i),1)-curage);
+        end
+        curage=tab(ind(i),1);
+        if isiso(ind(i)),k=k+1;else k=k-ex(ind(i))+1;end
+    end
+    p=p-log(rate)-bra1/rate-log(1-exp(-bra2/rate));
+    if k~=1,'errorHere',end
+    if bra1>=bra2,'errorThere',end
+    cur=iso(l);while tab(cur,2)>0,if ex(cur)==1,ex(cur)=2;break;end;ex(cur)=1;cur=tab(cur,2);end%Make all ancestors of current node active
+end
+
+function p=probSubtreeOLD(tab,rate)
+%Return the log-prior in Eq1 of Drummond et al (2002) Genetics 161:1307-1320
+%This is not quite right because of the truncation of the tree distribution TMRCA<INCUBATION_PERIOD
+n=size(tab,1)/2;
+p=-log(rate)*(n-1);
+tab(:,1)=max(tab(:,1))-tab(:,1);%convert times to ages
+[s,ind]=sort(tab(:,1));
+k=1;
+for i=2:length(s)-1
+    p=p-(k*(k-1)/(2*rate)*(s(i)-s(i-1)));
+    if isempty(find(tab(:,2)==ind(i))),k=k+1;else k=k-sum(tab(:,2)==ind(i))+1;end
+end
+if k~=1,'errorP',end
\ No newline at end of file
diff --git a/probTTree.m b/probTTree.m
new file mode 100644
index 0000000..c3326ab
--- /dev/null
+++ b/probTTree.m
@@ -0,0 +1,27 @@
+function prob=probTTree(ttree,popsize,beta,nu)
+%Calculate the log-probability of a transmission tree given popsize, beta and nu
+
+for i=1:size(ttree,1)
+    infector=ttree(i,3);
+    if infector==0,continue;end
+    if ttree(i,1)>ttree(infector,2)||ttree(i,1)<ttree(infector,1),prob=-inf;return;end
+end
+
+[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];
+prob=0;
+for i=2:length(times)
+    rate=nu*sir(2)+beta*sir(1)*sir(2);
+    prob=prob+log(rate)-rate*(times(i)-times(i-1));%;log(exppdf(times(i)-times(i-1),1/rate));
+    prob=prob+log(1/sir(2));
+    if types(i)==1
+        %Recovery
+        prob=prob+log(nu*sir(2)/rate);
+        sir(3)=sir(3)+1;sir(2)=sir(2)-1;
+    else
+        %Infection
+        prob=prob+log(1-nu*sir(2)/rate);
+        sir(1)=sir(1)-1;sir(2)=sir(2)+1;
+    end
+end
\ No newline at end of file
diff --git a/proposal.m b/proposal.m
new file mode 100644
index 0000000..5fb3bb5
--- /dev/null
+++ b/proposal.m
@@ -0,0 +1,61 @@
+function res=proposal(fulltree)
+n=sum(fulltree(:,2)==0&fulltree(:,3)==0);
+host=fulltree(:,4);
+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);
+
+%Choose a transmission event
+f=find(fulltree(:,2)>0&fulltree(:,3)==0);
+w=1+floor(rand*length(f));
+w=f(w);
+
+if w==size(fulltree,1)
+    %Move for the transmission to the index case
+    r=(rand-0.5)/1;
+    mini=-fulltree(end-1,1);
+    if r<mini,r=mini+(mini-r);end
+    fulltree(1:end-1,1)=fulltree(1:end-1,1)+r;
+else
+    %Move for the transmission to someone else
+    infector=host(w);
+    infected=host(fulltree(w,2));
+    %First remove the transmission node
+    f=fathers(w);
+    if fulltree(f,2)==w,fulltree(f,2)=fulltree(w,2);fathers(fulltree(w,2))=f;else fulltree(f,3)=fulltree(w,2);fathers(fulltree(w,2))=f;end
+    fulltree(w,:)=0;
+    %Second add it back on the path from infector to infected
+    islocs=zeros(size(fulltree,1),1);
+    path1=infector;islocs(path1)=1;
+    while path1<size(fulltree,1),path1=fathers(path1);islocs(path1)=1;end;
+    path2=infected;islocs(path2)=1;
+    while path2<size(fulltree,1),path2=fathers(path2);islocs(path2)=1-islocs(path2);end
+    locs=find(islocs);
+    lens=zeros(length(locs),1);
+    for i=1:length(locs)
+        lens(i)=fulltree(locs(i),1)-fulltree(fathers(locs(i)),1);
+    end
+    r=rand*sum(lens);
+    for i=1:length(locs)
+        if r>lens(i)
+            r=r-lens(i);
+        else
+            fulltree(w,1)=fulltree(locs(i),1)-r;
+            j=fathers(locs(i));
+            if fulltree(j,2)==locs(i),fulltree(j,2)=w;else fulltree(j,3)=w;end
+            fulltree(w,2)=locs(i);            
+            break;
+        end
+    end
+end
+
+%Reorder nodes chronologically
+[~,ind]=sort(fulltree(n+1:end,1),'descend');
+for i=n+1:size(fulltree,1)
+    for j=2:3
+        if fulltree(i,j)>n,fulltree(i,j)=n+find(ind==fulltree(i,j)-n);end
+    end
+end
+fulltree=fulltree([(1:n)';n+ind],:);
+res=[fulltree(:,1:3) hostFromFulltree(fulltree)];
\ No newline at end of file
diff --git a/ptreeFromFullTree.m b/ptreeFromFullTree.m
new file mode 100644
index 0000000..e9d9dd4
--- /dev/null
+++ b/ptreeFromFullTree.m
@@ -0,0 +1,11 @@
+function ptree=ptreeFromFullTree(tree)
+%Extract phylogenetic tree from a phylogenetic+transmission tree
+n=sum(tree(:,2)+tree(:,3)==0);
+tra=n+1;
+while tra<size(tree,1)
+    if tree(tra,3)~=0,tra=tra+1;continue;end
+    t=tree(:,2:3);f=find(t==tra);t(f)=tree(tra,2);f=find(t>tra);t(f)=t(f)-1;tree(:,2:3)=t;
+    tree=tree(setdiff(1:size(tree,1),tra),:);
+    tra=n+1;
+end
+ptree=tree(1:end-1,1:end-1);
\ No newline at end of file
diff --git a/readPTree.m b/readPTree.m
new file mode 100644
index 0000000..70917ad
--- /dev/null
+++ b/readPTree.m
@@ -0,0 +1,17 @@
+function ptree=readPTree(filename)
+
+ptree=phytreeread(filename);
+dist=get(ptree,'DISTANCES');
+poin=get(ptree,'POINTERS');
+n=get(ptree,'NUMLEAVES');
+ptree=zeros(n*2-1,3);
+ptree(n+1:end,2:3)=poin;
+for i=1:n*2-1
+    time=0;
+    j=i;
+    while j<2*n-1
+        time=time+dist(j);
+        j=n+find(poin(:,1)==j|poin(:,2)==j);
+    end
+    ptree(i,1)=time;
+end
diff --git a/simu.m b/simu.m
new file mode 100644
index 0000000..606e568
--- /dev/null
+++ b/simu.m
@@ -0,0 +1,31 @@
+function truth=simu(neg,popsize,beta,nu,ninf)
+%neg is Ne*g, ie the within-host effective population size (Ne) times  generation duration (g)
+%popsize is the size of the population
+%beta is the infectiveness
+%nu is the rate of recovery
+%ninf is the number of infected hosts, or 0 for any number
+
+if nargin<1,neg=100/365;end
+if nargin<2,popsize=1000;end
+if nargin<3,beta=0.001;end
+if nargin<4,nu=1;end
+if nargin<5,ninf=5;end
+
+%Create a transmission tree
+n=-1;
+while n~=ninf
+    ttree=makeTTree(popsize,beta,nu);
+    n=size(ttree,1);
+    if ninf==0,ninf=n;end
+end
+
+%Create a within-host phylogenetic tree for each infected host
+for i=1:n
+    times=[ttree(i,2);ttree(find(ttree(:,3)==i),1)]-ttree(i,1);
+    wtree{i}=withinhost(times,neg);
+end
+
+%Glue these trees together
+truth=glueTrees(ttree,wtree);
+truth(:,1)=truth(:,1)+2005;%Epidemic started in 2005
+timeLastRem=max(truth(:,1));
\ No newline at end of file
diff --git a/ttreeFromFullTree.m b/ttreeFromFullTree.m
new file mode 100644
index 0000000..f5119a0
--- /dev/null
+++ b/ttreeFromFullTree.m
@@ -0,0 +1,14 @@
+function ttree=ttreeFromFullTree(fulltree)
+%Takes in a fulltree, and extracts the ttree
+host=fulltree(:,4);
+ttree=fulltree(fulltree(:,2)==0&fulltree(:,3)==0,1);
+n=length(ttree);
+ttree=[zeros(length(ttree),1) ttree zeros(length(ttree),1)];
+for i=1:n
+    j=i;
+    while host(j)==i
+        j=find(fulltree(:,2)==j | fulltree(:,3)==j);
+    end
+    ttree(i,1)=fulltree(j,1);
+    ttree(i,3)=host(j);
+end
diff --git a/tutorial.m b/tutorial.m
new file mode 100644
index 0000000..b84da9f
--- /dev/null
+++ b/tutorial.m
@@ -0,0 +1,31 @@
+%Initialisation of random number generators
+rand('state',0);randn('state',0);
+
+%Simulation of an outbreak
+outbreak=simu(100/365,1000,0.001,1,10);
+
+%Plot the outbreak as a colored phylogeny
+plotBothTree(outbreak);
+
+%Extract and plot just the transmission tree
+plotTTree(ttreeFromFullTree(outbreak));
+
+%Extract and plot the phylogeny
+plotPTree(ptreeFromFullTree(outbreak));
+
+%Inference of transmission tree given a phylogeny
+ptree=ptreeFromFullTree(outbreak);
+timeLastSeq=max(outbreak(:,1));
+result=infer(ptree,10000,timeLastSeq);
+
+%Trace of log-posterior
+plot([result(:).pTTree]+[result(:).pPTree]);
+
+%Trace of parameter beta
+plot([result(:).beta]);
+
+%Histogram of probabilities to be the source case
+hist([result(end/2:end).source],1:10);
+
+%Plot configuration after first half of the MCMC
+plotBothTree(result(end/2).tree);
\ No newline at end of file
diff --git a/withinhost.m b/withinhost.m
new file mode 100644
index 0000000..793f8b9
--- /dev/null
+++ b/withinhost.m
@@ -0,0 +1,44 @@
+function [nodes,prob]=withinhost(times,neg)
+%Input=times at which N samples are taken (counted forward in time from infection time)
+%Output=array of size (2N)*3 where each row is a node, the first column indicate the date of the node and the last two
+%columns indicate the two children
+%This array has internal nodes sorted in order of most recent to most ancient node (and remains so during the algorithm)
+%The last node corresponds to infection time and only has one child
+%neg is the product of the within-host effective population size and the generation duration in days 
+prob=0;
+[tim,ind]=sort(times,'descend');
+n=length(tim);
+nodes=[0,ind(1),0];%Start with one node at time 0 and with the first isolate connected to it
+i=2;
+while i<=n%Graft branches one by one
+    r=-log(rand)*neg;
+    curt=tim(i);%Current time: start with date of isolate and go back in time until coalescence happens
+    fi=find(nodes(:,1)<curt,1);
+    for j=fi:size(nodes,1)
+        if r>(curt-nodes(j,1))*(i-j)
+            r=r-(curt-nodes(j,1))*(i-j);
+            curt=nodes(j,1);
+        else
+            curt=curt-r/(i-j);%Found the time for grafting
+            prob=prob+log(exppdf(r,neg));
+            r=0;break;
+        end
+    end
+    if r>0,continue;end%Rejection sampling forcing trees to be smaller than duration of infection
+    %Create new node
+    a=nodes(:,2:3);a(a>=j+n)=a(a>=j+n)+1;nodes(:,2:3)=a;%Renumbering according to table insertion in next line
+    nodes=[nodes(1:j-1,:);curt ind(i) 0;nodes(j:end,:)];
+    %Now choose on which branch to regraft amongst the branches alive at time curt
+    no=j;side=2;
+    prob=prob+log(1/(size(nodes,1)-j));
+    w=1+floor(rand*(size(nodes,1)-j));
+    while w>0
+        no=no+side-1;side=3-side;
+        if nodes(no,side+1)<=n || (nodes(no,side+1)>n && nodes(nodes(no,side+1)-n,1)>curt),w=w-1;end
+    end
+    nodes(j,3)=nodes(no,side+1);
+    nodes(no,side+1)=n+j;
+    i=i+1;
+end
+nodes=[zeros(n,3);nodes];
+nodes(1:n,1)=times;
\ No newline at end of file