1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349
|
r"""
Clustal format (:mod:`skbio.io.format.clustal`)
===============================================
.. currentmodule:: skbio.io.format.clustal
Clustal format (``clustal``) stores multiple sequence alignments. This format
was originally introduced in the Clustal package [1]_.
Format Support
--------------
**Has Sniffer: Yes**
+------+------+---------------------------------------------------------------+
|Reader|Writer| Object Class |
+======+======+===============================================================+
|Yes |Yes |:mod:`skbio.alignment.TabularMSA` |
+------+------+---------------------------------------------------------------+
Format Specification
--------------------
A clustal-formatted file is a plain text format. It can optionally have a
header, which states the clustal version number. This is followed by the
multiple sequence alignment, and optional information about the degree of
conservation at each position in the alignment [2]_.
Alignment Section
^^^^^^^^^^^^^^^^^
Each sequence in the alignment is divided into subsequences each at most
60 characters long. The sequence identifier for each sequence precedes each
subsequence. Each subsequence can optionally be followed by the cumulative
number of non-gap characters up to that point in the full sequence (not
included in the examples below). A line containing conservation information
about each position in the alignment can optionally follow all of the
subsequences (not included in the examples below).
.. note:: scikit-bio ignores conservation information when reading and does not
support writing conservation information.
.. note:: When reading a clustal-formatted file into an
``skbio.alignment.TabularMSA`` object, sequence identifiers/labels are
stored as ``TabularMSA`` index labels (``index`` property).
When writing an ``skbio.alignment.TabularMSA`` object as a clustal-formatted
file, ``TabularMSA`` index labels will be converted to strings and written
as sequence identifiers/labels.
Format Parameters
-----------------
The only supported format parameter is ``constructor``, which specifies the
type of in-memory sequence object to read each aligned sequence into. This must
be a subclass of ``GrammaredSequence`` (e.g., ``DNA``, ``RNA``, ``Protein``)
and is a required format parameter. For example, if you know that the clustal
file you're reading contains DNA sequences, you would pass ``constructor=DNA``
to the reader call.
Examples
--------
Assume we have a clustal-formatted file of RNA sequences:
.. code-block:: none
CLUSTAL W (1.82) multiple sequence alignment
abc GCAUGCAUCUGCAUACGUACGUACGCAUGCAUCA
def ----------------------------------
xyz ----------------------------------
abc GUCGAUACAUACGUACGUCGUACGUACGU-CGAC
def ---------------CGCGAUGCAUGCAU-CGAU
xyz -----------CAUGCAUCGUACGUACGCAUGAC
We can use the following code to read the clustal file into a ``TabularMSA``:
>>> from skbio import TabularMSA, RNA
>>> clustal_f = ['CLUSTAL W (1.82) multiple sequence alignment\n',
... '\n',
... 'abc GCAUGCAUCUGCAUACGUACGUACGCAUGCA\n',
... 'def -------------------------------\n',
... 'xyz -------------------------------\n',
... '\n',
... 'abc GUCGAUACAUACGUACGUCGGUACGU-CGAC\n',
... 'def ---------------CGUGCAUGCAU-CGAU\n',
... 'xyz -----------CAUUCGUACGUACGCAUGAC\n']
>>> msa = TabularMSA.read(clustal_f, constructor=RNA)
>>> msa
TabularMSA[RNA]
--------------------------------------------------------------
Stats:
sequence count: 3
position count: 62
--------------------------------------------------------------
GCAUGCAUCUGCAUACGUACGUACGCAUGCAGUCGAUACAUACGUACGUCGGUACGU-CGAC
----------------------------------------------CGUGCAUGCAU-CGAU
------------------------------------------CAUUCGUACGUACGCAUGAC
>>> msa.index
Index(['abc', 'def', 'xyz'], dtype='object')
We can use the following code to write a ``TabularMSA`` to a clustal-formatted
file:
>>> from io import StringIO
>>> from skbio import DNA
>>> seqs = [DNA('ACCGTTGTA-GTAGCT', metadata={'id': 'seq1'}),
... DNA('A--GTCGAA-GTACCT', metadata={'id': 'sequence-2'}),
... DNA('AGAGTTGAAGGTATCT', metadata={'id': '3'})]
>>> msa = TabularMSA(seqs, minter='id')
>>> msa
TabularMSA[DNA]
----------------------
Stats:
sequence count: 3
position count: 16
----------------------
ACCGTTGTA-GTAGCT
A--GTCGAA-GTACCT
AGAGTTGAAGGTATCT
>>> msa.index
Index(['seq1', 'sequence-2', '3'], dtype='object')
>>> fh = StringIO()
>>> _ = msa.write(fh, format='clustal')
>>> print(fh.getvalue()) # doctest: +NORMALIZE_WHITESPACE
CLUSTAL
<BLANKLINE>
<BLANKLINE>
seq1 ACCGTTGTA-GTAGCT
sequence-2 A--GTCGAA-GTACCT
3 AGAGTTGAAGGTATCT
<BLANKLINE>
<BLANKLINE>
References
----------
.. [1] http://www.sciencedirect.com/science/article/pii/0378111988903307
.. [2] http://web.mit.edu/meme_v4.9.0/doc/clustalw-format.html
"""
# ----------------------------------------------------------------------------
# Copyright (c) 2013--, scikit-bio development team.
#
# Distributed under the terms of the Modified BSD License.
#
# The full license is in the file COPYING.txt, distributed with this software.
# ----------------------------------------------------------------------------
from skbio.io import create_format, ClustalFormatError
from skbio.alignment import TabularMSA
clustal = create_format('clustal')
def _label_line_parser(record):
"""Returns dict mapping list of data to labels, plus list with field order.
Field order contains labels in order encountered in file.
NOTE: doesn't care if lines are out of order in different blocks. This
should never happen anyway, but it's possible that this behavior should
be changed to tighten up validation.
"""
labels = []
result = {}
for line in record:
split_line = line.strip().rsplit(None, 1)
if len(split_line) == 2:
key, val = split_line
else:
raise ClustalFormatError(
"Failed to parse sequence identifier and subsequence from "
"the following line: %r" % line)
if key in result:
result[key].append(val)
else:
result[key] = [val]
labels.append(key)
return result, labels
def _is_clustal_seq_line(line):
"""Returns True if line starts with a non-blank character but not 'CLUSTAL'
Useful for filtering other lines out of the file.
"""
return line and (not line[0].isspace()) and\
(not line.startswith('CLUSTAL')) and (not line.startswith('MUSCLE'))
def _delete_trailing_number(line):
"""Deletes trailing number from a line.
WARNING: does not preserve internal whitespace when a number is removed!
(converts each whitespace run to a single space). Returns the original
line if it didn't end in a number.
"""
pieces = line.split()
try:
int(pieces[-1])
return ' '.join(pieces[:-1])
except ValueError: # no trailing numbers
return line
def _check_length(data, labels, num_seqs_check=None):
"""
Checks the lengths of the clustal sequences to make
sure that they are lining up right
num_seqs_check: The number of sequences to check
Return True if all of the subsequence lengths are equal
or if data is empty
Return False if one of the subsequence lengths differs
"""
if len(labels) == 0:
return True
num_subseqs = len(data[labels[0]])
if num_seqs_check is None:
num_seqs_check = num_subseqs
else:
if num_seqs_check > num_subseqs:
num_seqs_check = num_subseqs
subseq_length = len(data[labels[0]][0])
end_lengths = set() # subsequence lengths at end of file
for i in range(num_seqs_check):
for label in labels:
seq = data[label][i]
if len(seq) > subseq_length:
return False
elif i+1 == num_subseqs: # Last subsequence
end_lengths.add(len(seq))
elif len(seq) < subseq_length:
return False
# All trailing subsequences must be the same
if len(end_lengths) > 1:
return False
return True
@clustal.sniffer()
def _clustal_sniffer(fh):
# Strategy
# The following conditions preclude a file from being clustal
# * It is an empty file
# * The whole sequences have differing lengths
# * The sub-sequences have differing lengths
# * One of the sequence ids is not immediately
# followed by a subsequence
empty = True
if fh.read(7) != 'CLUSTAL':
return False, {}
fh.seek(0)
try:
records = map(_delete_trailing_number,
filter(_is_clustal_seq_line, fh))
data, labels = _label_line_parser(records)
if len(data) > 0:
empty = False
# Only check first 50 sequences
aligned_correctly = _check_length(data, labels, 50)
if not aligned_correctly:
raise ClustalFormatError("Sequences not aligned properly")
except ClustalFormatError:
return False, {}
return not empty, {}
@clustal.writer(TabularMSA)
def _tabular_msa_to_clustal(obj, fh):
if not obj.index.is_unique:
raise ClustalFormatError(
"TabularMSA's index labels must be unique.")
clen = 60 # Max length of clustal lines
seqs = [str(s) for s in obj]
names = [str(label) for label in obj.index]
nameLen = max(map(len, names))
seqLen = max(map(len, seqs))
fh.write('CLUSTAL\n\n\n')
for i in range(0, seqLen, clen):
for label, seq in zip(names, seqs):
name = ('{:<%d}' % (nameLen)).format(label)
fh.write("%s\t%s\n" % (name, seq[i:i+clen]))
fh.write("\n")
@clustal.reader(TabularMSA)
def _clustal_to_tabular_msa(fh, constructor=None):
r"""yields labels and sequences from msa (multiple sequence alignment)
Parameters
----------
fh : open file object
An open Clustal file.
Returns
-------
skbio.TabularMSA
MSA containing aligned sequences.
Raises
------
skbio.util.exception.ClustalFormatError
If the sequences in `fh` don't have the same sequence length
or if the sequence ids don't properly match with the subsequences
Notes
-----
Skips any line that starts with a blank.
``_clustal_to_tabular_msa`` preserves the order of the sequences from the
original file. However, it does use a dict as an intermediate, so
two sequences can't have the same label. This is probably OK since
Clustal will refuse to run on a FASTA file in which two sequences have
the same label, but could potentially cause trouble with manually
edited files (all the segments of the conflicting sequences would
be interleaved, possibly in an unpredictable way).
If the lines have trailing numbers (i.e. Clustal was run with
`-LINENOS=ON`), silently deletes them. Does not check that the numbers
actually correspond to the number of chars in the sequence printed so far.
References
----------
.. [1] Thompson JD, Higgins DG, Gibson TJ, "CLUSTAL W: improving the
sensitivity of progressive multiple sequence alignment through sequence
weighting, position-specific gap penalties and weight matrix choice.
Thompson", Nucleic Acids Res. 1994 Nov 11;22(22):4673-80.
"""
if constructor is None:
raise ValueError("Must provide `constructor`.")
records = map(_delete_trailing_number,
filter(_is_clustal_seq_line, fh))
data, labels = _label_line_parser(records)
aligned_correctly = _check_length(data, labels)
if not aligned_correctly:
raise ClustalFormatError("Sequences not aligned properly")
seqs = []
for label in labels:
seqs.append(constructor(''.join(data[label])))
return TabularMSA(seqs, index=labels)
|