iqbal-lab-org / make_prg

Code to create a PRG from a Multiple Sequence Alignment file
Other
21 stars 7 forks source link

Properly handling Ns in make_prg #60

Closed leoisl closed 1 year ago

leoisl commented 1 year ago

Our handling of N bases in make_prg has not been appropriate I guess. An extreme example comes from @Danderson123 . He has an MSA with 789 alleles, and just a single allele has two runs of 10 Ns. At some point during the recursive cluster and collapse algorithm, we get this part of the allele isolated, and we skip this gene:

2023-07-11 10:09:41.175 | WARNING  | make_prg.subcommands.from_msa:process_MSA:148 - Skipping building PRG for alsB. Error: All sequences in this slice contained N. Redo sequence curation.
Sequences: ['GAANNNNNNNNNN']

I think this might be worst than erroring out, specially in cases where we build PRGs from thousands of genes, this just gets buried in the log files, even though it is a WARNING log message...

The rationale to refuse to build the PRG from MSAs where we have full columns with N might be good, but the recursive nature of this problem makes it hard to generalise (e.g. we should refuse building a graph from a 10-allele MSA if some columns are all Ns, but if this was a 1000-allele MSA and just these had 10 had the Ns, it would be ok...). Anyway, we chatted and it seemed to us the best solution would be to replace N with the consensus of each column. If the whole column is N and/or gaps, we replace it by a random base, with the seed initialised with the alignment sequences, i.e. this replacement is reproducible if the same alignment is given. This is done upfront when loading the MSA so that we take the whole MSA as the context to find the consensus for each column.

This is a small PR, the huge number of additions comes from the integration tests I added (I used the 3 genes @Danderson123 pointed out in his minimum working example).

This should fix https://github.com/rmcolq/pandora/issues/333

codecov[bot] commented 1 year ago

Codecov Report

Merging #60 (7454a89) into dev (80aca37) will increase coverage by 0.00%. The diff coverage is 100.00%.

@@           Coverage Diff           @@
##              dev      #60   +/-   ##
=======================================
  Coverage   99.73%   99.73%           
=======================================
  Files          18       18           
  Lines        1494     1525   +31     
=======================================
+ Hits         1490     1521   +31     
  Misses          4        4           
Impacted Files Coverage Δ
make_prg/utils/io_utils.py 98.41% <100.00%> (+0.23%) :arrow_up:
make_prg/utils/seq_utils.py 100.00% <100.00%> (ø)
leoisl commented 1 year ago

Hey @mbhall88 , thanks for the comments. I addressed both your comments, added unit tests covering the two changes, and further refactored the get_majority_consensus_from_MSA() function.