1. Load in the aligned sequences from the TB cluster in FASTA format using the ‘ape’ package.

sequences <- read.dna("TB_Cluster.fasta",as.matrix = T,format = "fasta")

2. Use the dist.dna() function in the ‘ape’ package to calculate the pairwise SNP distance between sequences.

distance <- dist.dna(sequences,model = "N")

3. Define a SNP threshold as the maximum distance allowed to link sequences.

threshold <- 5  # maximum allowed distance for linking

4. We must now transform the distance matrix to a data frame with the columns: 1. Column name, 2. Row name, 3. Distance value. We can also ignore the diagonal (where the row name = column name).

distance_mat <- as.matrix(distance)

distance_transform <- as.data.frame(distance_mat) %>%
  mutate(rowname = rownames(distance_mat)) %>%
  pivot_longer(
    -rowname,
    names_to = "colname",
    values_to = "value"
  ) %>%
  filter(!is.na(value), rowname != colname) %>% 
  filter(match(rowname, rownames(distance_mat)) > match(colname, colnames(distance_mat))) # ignore diagonal

5. Extract the rows with a value <= the SNP threshold.

close_sequences <- distance_transform[which(distance_transform$value <= threshold),]

6. Build an igraph object using the row and column names (columns 1 and 2)

g <- graph_from_edgelist(as.matrix(close_sequences[,c(1,2)]), directed = FALSE)

7. Extract the clusters and build a results table with the cluster designations.

clusters <- components(g)

results <- data.frame(Names = names(clusters$membership), Cluster = as.numeric(clusters$membership))

How many clusters do you find? Are all sequences assigned a cluster?

Repeat this analysis with a higher SNP threshold, what differences do you find?

1. Load in the timed phylogeny produced from BEAST2 and the table containing the names and collection dates for each sequence.

tree <- read.nexus("TB_Cluster.tree")

dates <- read.table("TB_Cluster.txt",header = T)

TransPhylo requires the date of the last collected sample in the dataset to scale the branches of the tree. This has to be in the decimal date format.

2. The following command will find the maximum date of the sequences and convert it to a decimal with the ‘lubridate’ package.

date_last_sample <- max(decimal_date(as.Date(dates$Date)))

3. TransPhylo next requires the tree to be converted to a specialised ptree format:

ptree <- ptreeFromPhylo(tree, dateLastSample = date_last_sample) 

As TransPhylo employs a Bayesian approach, we can include prior epidemiological parameters in the model and choose to fix them or update them through the MCMC sampling. The full list of parameters that can be specified can be seen if you enter help(“inferTTree”).

Here we will specify prior parameters for the generation time distribution (the time between cases in transmission chain becoming infectious) and the sampling time distribution (the time between infection and sampling of individuals), which are modelled with a gamma distribution. This best models these distributions for an organism such as TB that can have a long and variable latency (non-infectious) period. A value between 0 and 1 for the proportion of cases we believe we have sampled in our dataset can also be specified with the pi parameter and we can also set the number of MCMC iterations to run our model.

4. The following lines of script with set the prior parameters:

w.shape = 1.3 # Shape of the generation time distribution
w.scale = 3.3333 # Scale of the generation time distribution
ws.shape = 1.1 # Shape of the sampling time distribution
ws.scale = 2.5 # Shape of the sampling time distribution

pi = 0.3 # A starting value for the sampling proportion (here we think we have ~30% of all cases in our dataset)

num_iterations <- 5000 # The number of MCMC iterations

5. We can now run the main function of TransPhylo that uses an MCMC sampling algorithm to infer a transmission tree from the phylogenetic tree - inferTTree. This can also include our starting values for the epidemiological parameters specified in step 4:

set.seed(1234) # Set a seed for reproducibility

record <- inferTTree(ptree, w.shape = w.shape, w.scale = w.scale,
                    ws.shape = ws.shape, ws.scale = ws.scale,
                    startPi = pi, mcmcIterations = num_iterations)

Once inferTTree has completed, we will have a posterior collection of estimates sampled at each of our MCMC iterations. We want to check for convergence of our analysis, which is when the MCMC sampling has reached a stationary distribution and results of our model are meaningful.

6. We can inspect the trace plots of a few key metrics with the following command:

plot(record)

Question: Do you think your analysis has converged and the results are meaningful?

mat = computeMatWIW(record)

lattice::levelplot(mat,xlab='',ylab='')

ttree <- consTTree(record, burnin = 0.2, minimum = 0.5) # Burnin of 20% and minimum probability of 0.5.

plotTTree2(ttree, maxTime = 2020)

Questions:

Can you determine the difference between the black and white dots?

How many direct transmission events between sampled cases can you identify (you may need to increase the size of the plot using the Zoom button or by exporting the plot)?