diff --git a/Doc/Tutorial/chapter_pairwise.rst b/Doc/Tutorial/chapter_pairwise.rst index ce631887d..e4a6fa18b 100644 --- a/Doc/Tutorial/chapter_pairwise.rst +++ b/Doc/Tutorial/chapter_pairwise.rst @@ -2432,7 +2432,7 @@ The ``CodonAligner`` class defines additional attributes that refer to a number | | ``frameshift_plus_two_score`` | +----------------------------+---------------------------------+ -These parameters can be used as keyword arguments when creating the ``CodonAligner`` object; you can also use the following mnemonics shown in table :ref:`table:align-mnemonics`. The parameter values can also be set afterwards as attributes to the ``CodonAligner`` object. +These parameters can be used as keyword arguments when creating the ``CodonAligner`` object; you can also use the following mnemonics shown in table :ref:`table:codonalign-mnemonics`. The parameter values can also be set afterwards as attributes to the ``CodonAligner`` object. .. table:: Mnemonics to set codon alignment parameters :name: table:codonalign-mnemonics @@ -2457,7 +2457,7 @@ These parameters can be used as keyword arguments when creating the ``CodonAlign Using a ``CodonAligner`` to align a nucleotide sequence to an amino acid sequence ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -As an example, let's consider this section of the coding sequence of the Gag and Gag-Pol proteina of HIV-1: +As an example, let's consider this section of the coding sequence of the Gag and Gag-Pol proteins of HIV-1: .. cont-doctest @@ -2511,7 +2511,7 @@ well as on the ``alignment`` object: >>> alignment.score 10.0 -which is equal to the number of alignment amino acids (10) multiplied by the match score (1.0). +which is equal to the number of aligned amino acids (10) multiplied by the match score (1.0). For the Gag-Pol polyprotein, we find a -1 frameshift event: @@ -2553,7 +2553,7 @@ The alignment score for the Gag-Pol alignment is >>> alignment.score 7.0 -which is equal to the number of alignment amino acids (10) multiplied by the match score (1.0), minus the penalty for a -1 frameshift (-3.0). +which is equal to the number of aligned amino acids (10) multiplied by the match score (1.0), minus the penalty for a -1 frameshift (-3.0). By extracting rows of the alignment, you can obtain the aligned amino acid sequence and nucleotide sequence: @@ -2576,7 +2576,6 @@ As an example of a positive frameshift, here we align the amino acid sequence ar >>> from Bio.Seq import Seq >>> from Bio.SeqRecord import SeqRecord - >>> nuc_ecoli = Seq("TCCGACGTTCTTAGGGGGTATCTTTGACTACGACGCCAAGAAAGAGCGT") >>> nuc_ecoli = SeqRecord( ... Seq("TCCGACGTTCTTAGGGGGTATCTTTGACTACGACGCCAAGAAAGAGCGT"), id="Ecoli DNA" ... ) @@ -2649,20 +2648,20 @@ To calculate the alignment score without calculating the alignment itself, the 13.0 As before, you can obtain the aligned amino acid sequence and nucleotide sequence by -extracting rows from the alignment. Note that the A nucleotide read twice due to the --1 frameshift is included in the aligned nucleotide sequence: +extracting rows from the alignment. Note that the nucleotide skipped by the +1 +frameshift is included in the aligned nucleotide sequence, with a gap at the +corresponding position in the aligned amino acid sequence: .. cont-doctest .. code:: pycon - - >>> alignment[0] - 'QANFLREDLA' - >>> alignment[1] - 'CAGGCUAAUUUUUUAAGGGAAGAUCUGGCC' - >>> len(alignment[1]) - 30 + >>> alignment_ecoli[0] + 'SDVLRGYL-DYDAKKER' + >>> alignment_ecoli[1] + 'TCCGACGTTCTTAGGGGGTATCTTTGACTACGACGCCAAGAAAGAGCGT' + >>> len(alignment_ecoli[1]) + 49 The orthologous sequence in *Shigella dysenteriae* does not undergo a frame shift: @@ -2712,7 +2711,7 @@ to align the protein sequences to each other: >>> from Bio.Align import PairwiseAligner >>> protein_aligner = PairwiseAligner("blastp") >>> protein_alignments = protein_aligner.align(pro_ecoli, pro_shigella) - >>> len(alignments) + >>> len(protein_alignments) 1 >>> protein_alignment = protein_alignments[0] >>> print(protein_alignment) @@ -2784,7 +2783,7 @@ If you want to align just a few sequences, you may want to use the ``codon_align`` convenience function instead. This function can be used in two ways: - 1. With an amino acid sequence and a nucleotide sequence as the first and second rguments, respectively, this function creates a ``CodonAligner`` object, aligns the two sequences to each other, and returns an iterator over the alignments found. + 1. With an amino acid sequence and a nucleotide sequence as the first and second arguments, respectively, this function creates a ``CodonAligner`` object, aligns the two sequences to each other, and returns an iterator over the alignments found. For example, you can use ``codon_align`` as follows to align the amino acid sequence of the PrfB protein region in *Bacillus subtilis* against its corresponding genome sequence: @@ -2853,7 +2852,7 @@ The most important application of a codon alignment is to estimate the number of nonsynonymous substitutions per site (dN) and synonymous substitutions per site (dS). These can be calculated by the ``calculate_dn_ds`` function in ``Bio.Align.analysis``. This function -takes a pairwise codon alignment and input, as well as the optional +takes a pairwise codon alignment as input, as well as the optional arguments ``method`` specifying the calculation method, ``codon_table`` (defaulting to the Standard Code), the ratio ``k`` of the transition and transversion rates, and ``cfreq`` to specify the equilibrium codon @@ -2884,7 +2883,7 @@ frequency. Biopython currently supports three counting based methods - ``F61``: count the frequency of codons from the provided codon sequences, with a pseudocount of 0.1. -The ``calculate_dN_dS`` method can be applied to a pairwise codon +The ``calculate_dn_ds`` function can be applied to a pairwise codon alignment. In general, the different calculation methods will result in slightly different estimates for dN and dS: @@ -2908,12 +2907,14 @@ slightly different estimates for dN and dS: .. code:: pycon - >>> dN, dS = analysis.calculate_dn_ds(pairwise_codon_alignment, method="YN00") - >>> print(dN, dS) # doctest: +ELLIPSIS - 0.084526... 0.182723... >>> dN, dS = analysis.calculate_dn_ds(pairwise_codon_alignment, method="ML") >>> print(dN, dS) # doctest: +ELLIPSIS - 0.083040... 0.281696... + 0.083... 0.281... + +The ``ML`` estimates come from a numerical optimization, so only their first +three decimal places are shown. The ``YN00`` method is left out of this example because, +on this alignment, ``calculate_dn_ds`` with ``method="YN00"`` currently returns +``None`` instead of a pair of estimates. You can also calculate a matrix of dN and dS values, which is useful if you have a multiple sequence alignment: @@ -2942,10 +2943,10 @@ multiple sequence alignment: The objects ``dN`` and ``dS`` returned by ``calculate_dn_ds_matrix`` are instances of the ``DistanceMatrix`` class in -``Bio.Phylo.TreeConstruction``. This function only takes ``codon_table`` -as an optional argument. +``Bio.Phylo.TreeConstruction``. This function takes ``method`` and +``codon_table`` as optional arguments. -From these two sequences, you can create a dN tree and a dS tree using +From these two matrices, you can create a dN tree and a dS tree using ``Bio.Phylo.TreeConstruction``: .. cont-doctest @@ -2990,7 +2991,7 @@ will use Adh gene from fruit fly. The data includes 11 individuals from the nucleotide sequences are available in the ``Tests/codonalign`` directory as the files ``adh.aln`` and ``drosophila.fasta``, respectively, in the Biopython distribution. The function ``mktest`` in -``Bio.Align.analysis`` implements the Mcdonald-Kreitman test. +``Bio.Align.analysis`` implements the McDonald-Kreitman test. .. doctest ../Tests/codonalign lib:numpy diff --git a/UPSTREAM.md b/UPSTREAM.md index 6949a1ced..aa4bb1bdc 100644 --- a/UPSTREAM.md +++ b/UPSTREAM.md @@ -113,6 +113,8 @@ Verified unreported in the GitHub tracker and still present on upstream | 38 | A Gotoh local alignment with score 0 raises `IndexError` on `[0]`, while `len()` and iteration find its alignments | `GOTOH_LOCAL_ALIGN` raises `maximum` only by more than epsilon, so it can stay 0 while cells scoring exactly epsilon are end points; then [`if (maximum == 0) M[0][0].path = DONE;`](https://github.com/biopython/biopython/blob/f68b896031088a8ad11ae49f5e903d49228fd61f/Bio/Align/_pairwisealigner.c#L5480) (`master` `f68b89603`) starts the generator exhausted. `len()` counts the end points, and `PathGenerator_reset`, which iteration calls, clears `DONE`. `PairwiseAligner(mode="local", mismatch_score=-1, open_gap_score=-2, extend_gap_score=-1, epsilon=1).align("GAT", "TA")`: `len()` 2, `list()` 2 alignments, `[0]` `IndexError`. Smith-Waterman sets `NONE` there and has no end points at score 0, so it is consistent. Waterman-Smith-Beyer local has the same line ([L5982](https://github.com/biopython/biopython/blob/f68b896031088a8ad11ae49f5e903d49228fd61f/Bio/Align/_pairwisealigner.c#L5982)) and the same defect: with gap score functions `-2 - 0.5 * n` and `epsilon=1`, `align("A", "A")` has `len()` 1 and `align("GAT", "TA")` 2, and `[0]` raises `IndexError` for both. Both reproduced on PyPI 1.88. Tracker checked 2026-10-20 ("PairwiseAligner", "local alignment", "PairwiseAligner IndexError", "PairwiseAligner epsilon", and for Waterman-Smith-Beyer "Waterman-Smith-Beyer", "Waterman", "gap score function"): nothing filed. **Fixed here in [#238](https://github.com/dbolser/BioPAIthon/pull/238)** for Gotoh local, and in [#242](https://github.com/dbolser/BioPAIthon/pull/242) for Waterman-Smith-Beyer local. The same `maximum` ([`SELECT_TRACE_GOTOH_LOCAL_ALIGN`](https://github.com/biopython/biopython/blob/f68b896031088a8ad11ae49f5e903d49228fd61f/Bio/Align/_pairwisealigner.c#L4283)) also makes the reported score lag `score()` by up to epsilon, and as the first rise unmarks nothing, end points found before it stay: with the scores above, `"ACGT"` against itself gives score 4 and one path, `A` with `A`; `"AACAA"` against `"AAGAA"` gives score 2 (`score()` 3) and six paths. Separately, `GOTOH_LOCAL_SCORE` takes its maximum from Ix and Iy too, and lets a gap scoring 0 go on ([`SELECT_SCORE_LOCAL3`](https://github.com/biopython/biopython/blob/f68b896031088a8ad11ae49f5e903d49228fd61f/Bio/Align/_pairwisealigner.c#L4161)), so with a positive gap score `score()` reports a score no local alignment has: `match_score=1, mismatch_score=-1, open_gap_score=-0.5, extend_gap_score=0.5`, `"AAA"` against `"CCC"`: `score()` 0.5, `align()` 0. Both reproduced on PyPI 1.88. Waterman-Smith-Beyer local has the same lag ([L4393](https://github.com/biopython/biopython/blob/f68b896031088a8ad11ae49f5e903d49228fd61f/Bio/Align/_pairwisealigner.c#L4393)). **Fixed here in [#244](https://github.com/dbolser/BioPAIthon/pull/244)** for Gotoh local, as #239 did for Smith-Waterman. In Waterman-Smith-Beyer local, with `mismatch_score=-1`, gap score functions `-2 - (n - 1)` and `epsilon=1`, `"GAT"` against `"TA"` scores 0 (`score()` 1), `"ACGT"` against itself gives one path, `A` with `A`, and `"AACAA"` against `"AAGAA"` scores 2 (`score()` 3) with six paths. Its `score()` agrees with brute force, as local mode now rejects positive gap scores. Reproduced on PyPI 1.88. **Fixed here in [#245](https://github.com/dbolser/BioPAIthon/pull/245)** for Waterman-Smith-Beyer local. | | 39 | `PairwiseAligner` counts past an overflow, and Smith-Waterman's score, count and paths disagree | Two defects in [`Bio/Align/_pairwisealigner.c`](https://github.com/biopython/biopython/blob/f68b896031088a8ad11ae49f5e903d49228fd61f/Bio/Align/_pairwisealigner.c) on `master` `f68b89603`, both run on PyPI 1.86 and 1.88. `SAFE_ADD` (line 32) checks only the sum for the overflow marker, -1, so an overflowed count added to one that had not gives a wrong count instead of `OverflowError`, in every algorithm: global, all scores 0, `"A" * 37` against `"C" * 37` gives `len()` 8001394856851978725, where the true count is about 2e27. Smith-Waterman (`SELECT_TRACE_SMITH_WATERMAN_HVD`, line 4213) takes its score from every cell, including those a gap ends in, but marks end points only where a diagonal move ends, and unmarks earlier end points only on a rise of more than epsilon. So with a positive gap score `align()` and `score()` report a score no alignment has, the stale end points stay, and `PathGenerator_next_smithwaterman` raises "Unexpected trace" at one left without its diagonal trace (`"AAA"` against `"CAACC"`, gap score 0.5: score 2.5, the first path scores 1, `[1]` raises). With `epsilon=1`, `"CC"` against itself yields three paths scoring 1 for a score of 2. Tracker checked 2026-10-20: nothing filed. Report as two issues; whether local alignments may start with a gap when gaps score positive is upstream's call (see the fork PR). This fork has since chosen to reject positive gap scores in local mode (`ValueError`, see `DEPRECATED.rst`); offer that as one option, not as the fix. **Fixed here in [#239](https://github.com/dbolser/BioPAIthon/pull/239).** | | 40 | The EMBL writer has the same stray-line-break defect as [#3885](https://github.com/biopython/biopython/issues/3885), and [#3911](https://github.com/biopython/biopython/pull/3911) does not touch it | `EmblWriter.write_record` passes each feature straight to `_write_feature` ([`Bio/SeqIO/InsdcIO.py:1522`](https://github.com/biopython/biopython/blob/fc0c156786e86862112a05a720b6f4f47023f36a/Bio/SeqIO/InsdcIO.py#L1522) on `master` `fc0c15678`), so a qualifier value containing a newline, such as the SnapGene parser passes through, is written across two lines with the continuation unindented, and `SeqIO.read` of the writer's own output raises `ValueError: Problem with 'CDS' feature`. A record with `qualifiers={"note": ["line one\nline two"]}` written as `"embl"` or `"imgt"` reproduces it. #3911 fixes `GenBankWriter.write_record` only; the fix belongs on the shared `_InsdcWriter` base. Tracker checked 2026-10-22 ("EMBL line break", "EMBL writer newline qualifier", "EmblWriter", "qualifier newline"): nothing filed beyond #3885 itself, which names GenBank only. Comment on #3885 rather than a new issue. **Fixed here in [#246](https://github.com/dbolser/BioPAIthon/pull/246)** by moving the cleaning to `_InsdcWriter`, so both writers share it. | +| 41 | The Tutorial's codon-alignment section shows the HIV rows under the *E. coli* heading, checks `len(alignments)` after building `protein_alignments`, cites the wrong mnemonics table, has two typos, and labels LWL85's dN/dS numbers as `YN00` | Documentation only, all from [#5376](https://github.com/biopython/biopython/pull/5376) (2026-10-16). The copied *E. coli* block has an extra blank line after `.. code:: pycon`, so `Tests/test_Tutorial.py` never extracts it and upstream's CI never ran it; the `len(alignments)` line passes only because the HIV iterator is still in scope. Still on `master` `f10e88c58` on 2026-10-28; tracker searched for `alignment_ecoli`, `codonalign-mnemonics` and a `None` from `calculate_dn_ds`, nothing. Detail in [#269](https://github.com/dbolser/BioPAIthon/issues/269), which mentions #5376 so the cross-reference shows on that PR's timeline. **Fixed here in [#270](https://github.com/dbolser/BioPAIthon/pull/270)**, together with five older slips in the same section (`calculate_dN_dS`, "and input", `calculate_dn_ds_matrix` "only takes `codon_table`", "these two sequences", "Mcdonald"). | +| 42 | `calculate_dn_ds(method="YN00")` returns `None` when its iteration does not converge, and two defects make it fail on the Tutorial's own alignment | On the 16-codon *E. coli*/*Shigella* alignment: (1) the fold-class normalisation `for i, j in zip(fold0_cnt, fold4_cnt)` stops at the shorter dict, so with no `A` among the 4-fold sites `fold0_cnt["A"]` stays a raw count and the initial kappa is -35.0 (0-fold estimate -42.4); (2) with that fixed (kappa 0.98) the three-difference pair TAT/ATC has every substitution path through a codon with zero observed count, `sum(path_prob)` is 0 at `analysis.py:764`, numpy warns `invalid value encountered in scalar divide`, and every pass is NaN; (3) the convergence test never passes and `_yn00` falls off the end of its 20-pass loop, returning `None` with no error or warning of its own. The same `path_prob` warning is quoted in [#1611](https://github.com/biopython/biopython/issues/1611); see also the two `YN00` lines in the October list below. One report once checked against PAML's `yn00.c`. Detail in [#269](https://github.com/dbolser/BioPAIthon/issues/269). Not fixed here. | Rows 24–29 come from the §1.1 re-scope in `IMPROVEMENTS.md`. They were checked on 2026-10-06 against upstream `master` `372c71069`, and again on 2026-10-08