LG_matrice wrong?

[MAdriaansens]

Hey,

In the paper about LG matrice generation it is implied that the diagonal values should be negative (Qxx = -sum(Qxy)). With the sum of all rows equal to 0.
But why isnt the diagonal of the matrice for LG in IQtree negative then?

Mick

[bqminh]

Internally I'm sure that the diagonal of Q is negative, for all models, not just LG. Where do you see it's not?

[evolbeginner]

Hey,

In the paper about LG matrice generation it is implied that the diagonal values should be negative (Qxx = -sum(Qxy)). With the sum of all rows equal to 0. But why isnt the diagonal of the matrice for LG in IQtree negative then?

Mick

Hi Mick,

Here is the LG taken from IQ-Tree. https://github.com/evolbeginner/phyloHessianWrapper/blob/master/substitution_model/regular/LG.dat

It's given as an exchangeability matrix (top) and equilibrium frequencies (the last row). For simplicity, the exchangeability matrix doesn't have a diagonal shown (19*19).

The following is the r code for converting it to a Q matrix.

## Recover full 20x20 Q from:
## 1) 190 lower-triangle exchangeabilities (no diagonal)
## 2) 20 equilibrium frequencies (pi)

# ---- paste your lower-triangle block as one string ----
lower_txt <- "
0.425093 
0.276818 0.751878 
0.395144 0.123954 5.076149 
2.489084 0.534551 0.528768 0.062556 
0.969894 2.807908 1.695752 0.523386 0.084808 
1.038545 0.363970 0.541712 5.243870 0.003499 4.128591 
2.066040 0.390192 1.437645 0.844926 0.569265 0.267959 0.348847 
0.358858 2.426601 4.509238 0.927114 0.640543 4.813505 0.423881 0.311484 
0.149830 0.126991 0.191503 0.010690 0.320627 0.072854 0.044265 0.008705 0.108882 
0.395337 0.301848 0.068427 0.015076 0.594007 0.582457 0.069673 0.044261 0.366317 4.145067  
0.536518 6.326067 2.145078 0.282959 0.013266 3.234294 1.807177 0.296636 0.697264 0.159069 0.137500 
1.124035 0.484133 0.371004 0.025548 0.893680 1.672569 0.173735 0.139538 0.442472 4.273607 6.312358 0.656604 
0.253701 0.052722 0.089525 0.017416 1.105251 0.035855 0.018811 0.089586 0.682139 1.112727 2.592692 0.023918 1.798853 
1.177651 0.332533 0.161787 0.394456 0.075382 0.624294 0.419409 0.196961 0.508851 0.078281 0.249060 0.390322 0.099849 0.094464 
4.727182 0.858151 4.008358 1.240275 2.784478 1.223828 0.611973 1.739990 0.990012 0.064105 0.182287 0.748683 0.346960 0.361819 1.338132 
2.139501 0.578987 2.000679 0.425860 1.143480 1.080136 0.604545 0.129836 0.584262 1.033739 0.302936 1.136863 2.020366 0.165001 0.571468 6.472279 
0.180717 0.593607 0.045376 0.029890 0.670128 0.236199 0.077852 0.268491 0.597054 0.111660 0.619632 0.049906 0.696175 2.457121 0.095131 0.248862 0.140825 
0.218959 0.314440 0.612025 0.135107 1.165532 0.257336 0.120037 0.054679 5.306834 0.232523 0.299648 0.131932 0.481306 7.803902 0.089613 0.400547 0.245841 3.151815 
2.547870 0.170887 0.083688 0.037967 1.959291 0.210332 0.245034 0.076701 0.119013 10.649107 1.702745 0.185202 1.898718 0.654683 0.296501 0.098369 2.188158 0.189510 0.249313
"

# ---- paste your 20 frequencies ----
pi_vec <- scan(text = "
0.079066 0.055941 0.041977 0.053052 0.012937 0.040767 0.071586 0.057337 0.022355 0.062157
0.099081 0.0646 0.022951 0.042302 0.04404 0.061197 0.053287 0.012066 0.034155 0.069146
", quiet = TRUE)

# Optional state labels (AA order if needed)
states <- c("A","R","N","D","C","Q","E","G","H","I","L","K","M","F","P","S","T","W","Y","V")

# ---- parse lower triangle ----
r_lower <- scan(text = lower_txt, quiet = TRUE)
n <- length(pi_vec)

stopifnot(length(r_lower) == n*(n-1)/2)

# Exchangeability matrix R (symmetric, diag 0)
R <- matrix(0, n, n)
R[lower.tri(R)] <- r_lower
R <- R + t(R)
dimnames(R) <- list(states, states)

# Rate matrix Q: q_ij = r_ij * pi_j (i != j), q_ii = -sum_{j!=i} q_ij
Q <- R
for (i in seq_len(n)) {
    for (j in seq_len(n)) {
        if (i != j) Q[i, j] <- R[i, j] * pi_vec[j]
    }
}
diag(Q) <- -rowSums(Q)

dimnames(Q) <- list(states, states)

# Optional: scale so expected rate = 1
scale_factor <- -sum(pi_vec * diag(Q))
Q_scaled <- Q / scale_factor

cat("Check row sums (should be ~0):\n")
print(round(rowSums(Q_scaled), 12))

cat("\nR (exchangeabilities):\n")
print(R)

You may with to try the following to tse the sum of the diag multiplied by the eq freq equal to -1.

sum(diag(Q_scaled) %*% pi_vec)

And, what you might be asking

diag(Q_scaled) + rowSums(Q_scaled - diag(diag(Q_scaled)))

best,
sishuo