Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

error during SNaQ score calculation: "strange qnet", "Bad Diamond I" and "qnet should have only one of the gammaz" #158

Open
cecileane opened this issue Mar 15, 2021 · 0 comments

Comments

@cecileane
Copy link
Member

This is from a network obtained with SNaQ (thanks to Damien Esquerré). The error comes when trying to re-calculate the expected quartet CFs for a goodness-of-fit test. Here is a reproducible example, using the attached file tableCFs.txt:

qCF = DataFrame!(CSV.File("tableCFs.txt"));
netstring = "((a,b):0.741,((((((c,(d)#H29:::0.893):0.039,(e,#H29:::0.107):0.039):1.682,((f,(g,#H25:0.0::0.331):0.519):0.591,((h,i):1.058,j):0.128):0.185):0.105,k):0.126,(((((l,(m)#H1:::0.897):0.714,(#H1:::0.103,((n,o):0.502,p):0.076):0.165):1.97,((q,((r,((s,t):0.499,(u)#H28:::0.765):0.199):0.134,(v,#H28:::0.235):0.437):1.889):0.0)#H30:0.308::0.976):0.247,#H30:0.0::0.024):0.367)#H25:0.003::0.669):1.142,(w,#H27:::0.273):0.25):0.24,(x)#H27:::0.727);"
net = readTopology(netstring)
topologyQPseudolik!(net,readTableCF(qCF)) # recalculate SNaQ score and all expected CFs

and here is the output with error:

between 49.0 and 838.0 gene trees per 4-taxon set
ERROR: strange qnet when net has node 5 Bad Diamond I: qnet should have only one of the gammaz if it does not have node, but it has two
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] updateHasEdge!(::PhyloNetworks.QuartetNetwork, ::HybridNetwork) at /Users/.../src/pseudolik.jl:445
 [3] extractQuartet!(::HybridNetwork, ::Quartet) at /Users/.../src/pseudolik.jl:492
  ...

The error is strange, in that it goes away easily. I have a hard time providing a simpler example. Here is a function to remove some taxon set from both the network and the data, to see if the error persists.

function rmtax(taxa; unroot=true)
  qCF = DataFrame!(CSV.File("tableCFs.txt"));
  net = readTopology(netstring)
  for taxon in taxa deleteleaf!(net, taxon; unroot=unroot); end
  filter!(row -> !any(row[i] in taxa for i in [:t1,:t2,:t3,:t4]), qCF)
  return topologyQPseudolik!(net, readTableCF(qCF))
end
rmtax([]) # same error as above: original full taxon set
rmtax(["x"]; unroot=false) # error after removing taxon x, without unrooting the network
rmtax(["x"]) # no error!!
rmtax(["h"]) # error
rmtax(["h","i"]) # no error
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant