-
Notifications
You must be signed in to change notification settings - Fork 18
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
Biopython update #157
base: master
Are you sure you want to change the base?
Biopython update #157
Changes from all commits
9493f88
3472a94
7ca2b37
ad552f4
fb35568
00260a1
00a497a
cb3a39f
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,6 +1,7 @@ | ||
""" | ||
Command Line Interface for opencadd | ||
""" | ||
|
||
import argparse | ||
import logging | ||
|
||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,6 +1,7 @@ | ||
""" | ||
Base class for all the superposition engines. | ||
""" | ||
|
||
import logging | ||
|
||
_logger = logging.getLogger(__name__) | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,6 +1,7 @@ | ||
""" | ||
Utilities for sequence alignment | ||
""" | ||
|
||
import logging | ||
|
||
import numpy as np | ||
|
@@ -77,6 +78,27 @@ def sequence_alignment(seq1: str, seq2: str, matrix: str, gap: int, local: bool | |
return alignment[0] | ||
|
||
|
||
def find_gap_character(seq): | ||
""" | ||
Finds the gap character of an alignment sequence | ||
|
||
Args: | ||
seq: str | ||
aligned sequence | ||
|
||
Returns: | ||
gap character if found: str | ||
None else | ||
""" | ||
if "-" in seq: | ||
return "-" | ||
elif "=" in seq: | ||
return "=" | ||
else: | ||
return "" | ||
# raise ValueError("No standard gap character found!") | ||
|
||
|
||
def fasta2select( | ||
fastafilename, | ||
ref_resids=None, | ||
|
@@ -112,9 +134,8 @@ def fasta2select( | |
dictionary with 'reference' and 'mobile' selection string | ||
|
||
""" | ||
protein_gapped = Bio.Alphabet.Gapped(Bio.Alphabet.IUPAC.protein) | ||
with open(fastafilename) as fasta: | ||
alignment = Bio.AlignIO.read(fasta, "fasta", alphabet=protein_gapped) | ||
alignment = Bio.AlignIO.read(fasta, "fasta") | ||
|
||
nseq = len(alignment) | ||
if nseq != 2: | ||
|
@@ -127,7 +148,7 @@ def fasta2select( | |
if orig_resids[iseq] is None: | ||
# build default: assume consecutive numbering of all | ||
# residues in the alignment | ||
GAP = a.seq.alphabet.gap_char | ||
GAP = find_gap_character(a.seq) | ||
length = len(a.seq) - a.seq.count(GAP) | ||
orig_resids[iseq] = np.arange(1, length + 1) | ||
else: | ||
|
@@ -136,7 +157,7 @@ def fasta2select( | |
if orig_segids[iseq] is None: | ||
# build default: assume consecutive numbering of all | ||
# residues in the alignment | ||
GAP = a.seq.alphabet.gap_char | ||
GAP = find_gap_character(a.seq) | ||
length = len(a.seq) - a.seq.count(GAP) | ||
orig_segids[iseq] = np.full(length, "") | ||
else: | ||
|
@@ -175,7 +196,8 @@ def resid_factory(alignment, seq2resids): | |
t = np.zeros((nseq, alignment.get_alignment_length()), dtype=int) | ||
s = np.zeros((nseq, alignment.get_alignment_length()), dtype=object) | ||
for iseq, a in enumerate(alignment): | ||
GAP = a.seq.alphabet.gap_char | ||
print(a.seq) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Leftover from debugging? |
||
GAP = find_gap_character(a.seq) | ||
indices = np.cumsum(np.where(np.array(list(a.seq)) == GAP, 0, 1)) - 1 | ||
t[iseq, :] = seq2resids[iseq][0][indices] | ||
s[iseq, :] = seq2resids[iseq][1][indices] | ||
|
@@ -190,8 +212,8 @@ def resid(nseq, ipos, t=t, s=s): | |
res_list = [] # collect individual selection string | ||
|
||
# should be the same for both seqs | ||
GAP = alignment[0].seq.alphabet.gap_char | ||
if GAP != alignment[1].seq.alphabet.gap_char: | ||
GAP = find_gap_character(a.seq) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Shouldn't this (Line 214) be GAP = alignment[0].seq.alphabet.gap_char
if GAP != alignment[1].seq.alphabet.gap_char: and the second line has changed to if GAP != find_gap_character(alignment[1].seq): so it looks like the first line must also be: GAP = find_gap_character(alignment[0].seq) instead of: GAP = find_gap_character(a.seq) |
||
if GAP != find_gap_character(alignment[1].seq): | ||
raise ValueError("Different gap characters in sequence 'target' and 'mobile'.") | ||
for ipos in range(alignment.get_alignment_length()): | ||
aligned = list(alignment[:, ipos]) | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Shouldn't the gap characters be removed explicitly now that BioPython doesn't recognize them?
From BioPython Documentation:
Another case where the alphabet was used was in declaring the gap character, by default
-
in the various Biopython sequence and alignment parsers. If you are using a different character, you will need to pass this to theSeq
object.replace()
method explicitly now: