Handle feature references another sequence

[peterjc]

This pull request addresses issue #808, including the work of Adam Sjøgren via the mailing list.

• I hereby agree to dual licence this and any previous contributions under both
  the Biopython License Agreement AND the BSD 3-Clause License.

• I have read the CONTRIBUTING.rst file, have run flake8 locally, and
  understand that AppVeyor and TravisCI will be used to confirm the Biopython unit
  tests and style checks pass with these changes.

• I have added my name to the alphabetical contributors listings in the files
  NEWS.rst and CONTRIB.rst as part of this pull request, am listed
  already, or do not wish to be listed. (This acknowledgement is optional.)

[codecov]

Codecov Report

Merging #2334 into master will increase coverage by 0.00%.
The diff coverage is 87.50%.

@@           Coverage Diff           @@
##           master    #2334   +/-   ##
=======================================
  Coverage   83.88%   83.89%           
=======================================
  Files         317      317           
  Lines       51548    51556    +8     
=======================================
+ Hits        43242    43251    +9     
+ Misses       8306     8305    -1     

Impacted Files	Coverage Δ
Bio/SeqFeature.py	87.12% <87.50%> (+0.38%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update db43482...585fc9a. Read the comment docs.

[peterjc]

Rebasing and dealing with conflicts and new release...

[peterjc]

I haven't forgotten about this, but we've been distracted with the Python 2/3 cleanup in the meantime.

[peterjc]

Rebasing again (and leaving out the dual licensing changes)...

[peterjc]

Code seems to be working fine, but have not found a good example yet - for instance none of the features here give the expected translation:

from Bio import SeqIO

files = ("FP885871.gbk", "FP885876.gbk")

references = SeqIO.to_dict(SeqIO.read(_, "gb") for _ in files)

for acc, ref in references.items():
    print(f"{acc} length {len(ref)} has {len(ref.features)} features")
    for f in ref.features:
        if f.type == "CDS" and len(f.location.parts) > 1:
            print(f"{f.type} feature at {f.location}")
            #print(f)
            expected = f.qualifiers["translation"][0]
            try:
                translation = f.extract(ref.seq, references).translate(cds=True)
            except:
                translation = f.extract(ref.seq, references).translate(to_stop=True)
            print(expected + "<-- expected")
            print(translation + "<-- computed", expected == translation)
            if "exception" in f.qualifiers and f.qualifiers["exception"] == ["RNA editing"]:
                continue
            assert expected == translation, f"{expected} vs {translation}"

Some (all?) of them have human readable comments explaining why (in more detail than just "RNA editing") at least.

[peterjc]

Sadly this example from ftp://ftp.ncbi.nlm.nih.gov/genomes/archive/old_refseq/Arabidopsis_thaliana/ used in the Tutorial from drawing tRNA on chromosomes does not appear to have any clear trans-splicing annotation:

filenames = [
    "CHR_I/NC_003070.gbk",
    "CHR_II/NC_003071.gbk",
    "CHR_III/NC_003074.gbk",
    "CHR_IV/NC_003075.gbk",
    "CHR_V/NC_003076.gbk",
]

references = SeqIO.index_db("Arabidopsis.idx", filenames, "genbank")

Perhaps a newer version might... or maybe this is just not an appropriate example.

[peterjc]

Zebrafish to the rescue - a well behaved example from Adam Sjøgren, BX465837.11 which has a couple of CDS with trans-splicing with BX537337.9

https://www.ncbi.nlm.nih.gov/nuccore/BX465837.11
https://www.ncbi.nlm.nih.gov/nuccore/BX537337.9

Should work in EMBL format too:

https://www.ebi.ac.uk/ena/browser/api/embl/BX465837.11
https://www.ebi.ac.uk/ena/browser/api/embl/BX537337.9

[peterjc]

@MarkusPiotrowski - would you be a good person to review this?

I'm hoping to turn the Zebrafish example into something for the Tutorial and/or docstrings.

[peterjc]

This is going to need some merge conflicts resolved now...

[peterjc]

@mdehoon do you want to review?

[mdehoon]

I guess these exceptions:

            raise ValueError("Feature references another sequence.")

in _shift and _flip can now be removed?

[peterjc]

Yes, as per my comment #808 (comment) in those cases _shift and _flip are relative to the local sequence, and thus are a no-operation relative to an external sequence.

[peterjc]

Thank you.