Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
53 changes: 27 additions & 26 deletions Doc/Tutorial/chapter_pairwise.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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:

Expand Down Expand Up @@ -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:

Expand All @@ -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"
... )
Expand Down Expand Up @@ -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:

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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:

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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:

Expand All @@ -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:
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down
2 changes: 2 additions & 0 deletions UPSTREAM.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading