Issue about Figure S1

[ArtPoon]

A technical question regarding the Figure S1 simulation: If you simulate under INDELible (which uses a slightly different substitution model, e.g., GY94 with possibly different equilibrium frequencies) and infer with FUBAR, what exactly are you comparing via RMSE? FUBAR was not really designed to estimate true point values of dN and dS at a site super-accurately; it is meant to estimate the posterior probability of a site belonging to a specific rate class on a grid. A much more sensible approach would be to calculate the Wasserstein distance between the simulated distribution of dN/dS and the inferred FUBAR grid density.

I could not find any statement in paper associated with FUBAR to the effect that it is "not really designed to estimate true point values of dN and dS at a site" with accuracy. The HyPhy 2.5 paper cites FUBAR as one of the available methods for measuring site-specific selection. In addition, Murrell et al. 2016 stated:

FUBAR was originally designed for identifying sites under positive selection, and the accuracy of θ̂ is relatively unimportant for site-specific inference, which governed some of the design choices behind FUBAR.

FUBAR also remains the recommended option on the Datamonkey webserver for detecting pervasive (non-episodic) site-specific selection.

I also think this statement is inaccurate: "[FUBAR] is meant to estimate the posterior probability of a site belonging to a specific rate class on a grid". The grid in FUBAR is a means to an end. It makes the estimation of site-specific selection more efficient by constraining dN and dS values to a fixed grid of combined values. In contrast, the general bivariate method used in Pond et al. needs to determine the dN and dS values associated with each rate classes, as well as its probability.

We could take the posterior distribution (grid values) from the FUBAR analyses and compare them to the "true" distributions from the model simulation. Since INDELible only simulates variation in dN/dS (omega), this distribution will be constant with respect to dS.

Alternatively, we could redo this analysis using SLAC on the same simulated alignments. (If we still have these.)

[horaciobam]

I agree that the statement about FUBAR not being good at estimating site-specific dN and dS could be misleading. The documentation for FUBAR indicates that the method reports mean posterior estimates for α (dS) and β (dN) at each codon site.

From our current implementation:

• From the FUBAR results, we compute the site-specific dN/dS ratio by taking the ratio between the posterior mean of beta and alpha for each codon site:
  

  surfaces/scripts/indelible_vs_fubar.R

  Lines 13 to 33 in 5e304a5

  	##################################
  	# Read FUBAR results
  	##################################
  	setwd("/home/laura/Projects/surfaces/data_find_tree_length/march21_all")
  	files <- Sys.glob("*.FUBAR.json")
  	names(files) <- files # so master list has filenames as names
  	master <- lapply(files, function(f) read_json(f, simplifyVector=TRUE))
  	names(master)
  	# Returns MLE estimates for each dataset as list of data frames
  	# Data frame res[[1]] includes alpha, beta, alpha=beta, LRT, P-value,
  	# Total branch length, codon.pos and dnds
  	res <- lapply(master, function(d) {
  	stats <- d$MLE
  	site.stats <- as.data.frame(stats$content[[1]][,1:6])
  	names(site.stats) <- stats$headers[,1]
  	site.stats$codon.pos<-row.names(site.stats)
  	site.stats$dnds<- site.stats$beta/site.stats$alpha
  	return(site.stats)
  	})

• For each simulation generated with INDELible, we then compare the true omega values per codon site to the inferred dnds ratios computed above:
  

  surfaces/scripts/indelible_vs_fubar.R

  Lines 65 to 92 in 5e304a5

  	#####################################
  	# Calculate RSME for all tree lengths
  	#####################################
  	out <- matrix(ncol = 3, nrow=length(res))
  	# From:
  	for(i in 1:length(res)){
  	name <- names(res)[[i]]
  	sel <- res[[i]]
  	sim.rates <- rates[[i]]
  	# --- get tree length from name ----
  	spl <- strsplit(name, "_")[[1]]
  	tree.l <- spl[2]
  	# --- calculate RSME ----
  	rsme <- sqrt(mean(sim.rates$omegas - sel$dnds)^2)
  	# --- store info ----
  	out[i,1] <-as.numeric(tree.l)
  	out[i,2] <-as.numeric(rsme)
  	out[i,3] <-as.numeric(spl[3])
  	}
  	out <- as.data.frame(out)
  	colnames(out) <- c("tree.length", "rmse", "rep")
  	out <- out[order(out$tree.length),]

My understanding of your approach 1 would be:

• Use the posterior distributions of FUBAR rather than the posterior mean estimates
• Estimate the distribution of omega on each INDELible's simulation
• Compare the omega distribution from INDELible to the posterior distributions of dN and dS from FUBAR

What is not clear to me is how to compare INDELible's omega against the grid? Can we assume than in iNDELible dS=1?

[ArtPoon]

The thing about comparing posterior distribution (fingerprint) from FUBAR to simulations is that we have shown that the granularity (smoothness) of the fingerprint is highly dependent on sample size, e.g., tree length. This makes it very difficult to use this method to determine what tree length is necessary to estimate fingerprints.

I would much much rather use SLAC to compute site-specific dN and dS values, given the reviewer's issue with using FUBAR for this purpose, and make the same comparison.

[horaciobam]

New results running SLAC on the simulated data:

out <- read.csv("march21_slac_rsme.csv")

plot(out$tree.length, out$rmse, log = "x",
     pch=19, las=1
     )

---

• Note: I believe we had a typo in our previous RSME estimates, as we were calculating it as sqrt(mean(sim.rates$omegas - sel$dnds)^2) (absolute mean error) instead of sqrt(mean((sim.rates$omegas - sel$dnds)^2)) (root mean squared error):
  

  surfaces/scripts/indelible_vs_fubar.R

  Lines 68 to 92 in 5e304a5

  	out <- matrix(ncol = 3, nrow=length(res))
  	# From:
  	for(i in 1:length(res)){
  	name <- names(res)[[i]]
  	sel <- res[[i]]
  	sim.rates <- rates[[i]]
  	# --- get tree length from name ----
  	spl <- strsplit(name, "_")[[1]]
  	tree.l <- spl[2]
  	# --- calculate RSME ----
  	rsme <- sqrt(mean(sim.rates$omegas - sel$dnds)^2)
  	# --- store info ----
  	out[i,1] <-as.numeric(tree.l)
  	out[i,2] <-as.numeric(rsme)
  	out[i,3] <-as.numeric(spl[3])
  	}
  	out <- as.data.frame(out)
  	colnames(out) <- c("tree.length", "rmse", "rep")
  	out <- out[order(out$tree.length),]

• I have addressed this problem in the script that I used for comparing SLAC versus INDELible
  

  surfaces/scripts/indelible_vs_slac.py

  Lines 90 to 101 in 13bd217

  	def compute_rmse(truth, estimate):
  	"""Compute true RMSE over non-missing estimated values."""
  	sqerr = []
  	for true_value, est_value in zip(truth, estimate):
  	if est_value is None:
  	continue
  	sqerr.append((true_value - est_value) ** 2)
  	if not sqerr:
  	return None, 0
  	return math.sqrt(sum(sqerr) / len(sqerr)), len(sqerr)

• Note 2: I am excluding sites where dS=0 and dN>0, which sometimes happens in short branches
  

  surfaces/scripts/indelible_vs_slac.py

  Lines 46 to 73 in 13bd217

  	def load_slac_sitewise(path):
  	"""
  	Extract sitewise dN/dS values from a HyPhy SLAC JSON file.
  	Rules:
  	- if dS > 0, use dN / dS
  	- if dS == 0 and dN == 0, use 0
  	- if dS == 0 and dN > 0, use None and drop from RMSE
  	"""
  	with open(path) as handle:
  	data = json.load(handle)
  	rows = data["MLE"]["content"]["0"]["by-site"]["AVERAGED"]
  	out = []
  	for row in rows:
  	ds = row[5]
  	dn = row[6]
  	if ds is None or dn is None:
  	out.append(None)
  	continue
  	if ds > 0:
  	out.append(dn / ds)
  	elif ds == 0 and dn == 0:
  	out.append(0.0)
  	else:
  	out.append(None)
  	return out

[horaciobam]

Comparing FUBAR, SLAC, and FEL

---

Calculating omega from FEL and SLAC outputs

• I think FEL is problematic for small trees. It tries to estimate site-specific synonymous and nonsynonymous rates directly by maximum likelihood, and small trees often do not contain enough substitutions for those sitewise estimates to be stable.

  • Short branches means fewer observed substitutions across the alignment
  • FEL works site by site, for each codon, it tries to estimate alpha and beta separately
  • If a site has very little change, there is not enough information to estimate both rates well. Thus we get alpha = 0, beta = 0, or unstable ratios like very large beta / alpha
  • That makes many sites effectively uninformative, thus I removed them when estimating RSME:

• Similarly, I drop sites in SLAC when dS == 0, or dN==0 and dS==0

• FUBAR does not seem to have this problem because it rarely produced zero or missing denominators (I guess this might be one of the great advantages of estimating site values from the discrete grid)

---

Examples of problematic trees

• Simulation sim_0.02_0_.cleaned.nwk, dropped 990/1000 in FEL (rsme=43.37), and 992/1000 in SLAC (rsme = 0.575)

• Simulation sim_0.04_6_.cleaned.nwk, dropped 982 sites and got rsme = 0.93 for both SLAC and FEL

• Simulation sim_5.0_9_.cleaned.nwk, dropped 145 sites in SLAC (rsme = 0.82) and 148 sites in FEL (rsme = 0.71)

• Simulation sim_0.04_4_.cleaned.nwk, 966 sites dropped in FEL (rsme = 211.8), and 969 sites dropped in SLAC (rsme = 0.61)