MSA

The MSA module of MIToS has utilities for working with Multiple Sequence Alignments of protein Sequences (MSA).

Features

• Read and write MSAs in Stockholm, FASTA or Raw format

• Handle MSA annotations

• Edit the MSA, e.g. delete columns or sequences, change sequence order, shuffling...

• Keep track of positions and annotations after modifications on the MSA

• Describe a MSA, e.g. mean percent identity, sequence coverage, gap percentage...

using MIToS.MSA

MIToS MSA and aligned sequences (aligned objects) are subtypes of AbstractMatrix{Residue}, because MSAs and sequences are stored as Matrix of Residues.

A MIToS aligned sequence is an AbstractMatrix{Residue} with only 1 row/sequence.

MSAs are stored as Matrix{Residue}. It's possible to use a NamedResidueMatrix as the most simple MSA with sequence identifiers and column names.

An AlignedSequence wraps a NamedResidueMatrix with only 1 row/sequence. The NamedArray stores the sequence name and original column numbers as Strings.

This type represent an aligned sequence, similar to AlignedSequence, but It also stores its Annotations.

This type represent an MSA, similar to MultipleSequenceAlignment, but It also stores Annotations. This annotations are used to store residue coordinates (i.e. mapping to UniProt residue numbers).

The Annotations type is basically a container for Dicts with the annotations of a multiple sequence alignment. Annotations was designed for storage of annotations of the Stockholm format.

MIToS also uses MSA annotations to keep track of:

• Modifications of the MSA (MIToS_...) as deletion of sequences or columns.

• Positions numbers in the original MSA file (column mapping: ColMap)

• Position of the residues in the sequence (sequence mapping: SeqMap)

Data structure to represent sequence clusters. The sequence data itself is not included.

This type defines the usual alphabet of the 20 natural residues and a gap character.

julia> GappedAlphabet()
MIToS.MSA.GappedAlphabet of length 21. Residues : res"ARNDCQEGHILKMFPSTWYV-"

This MSA type include a NamedArray wrapping a Matrix of Residues. The use of NamedArray allows to store sequence names and original column numbers as Strings, and fast indexing using them.

Use NoClustering() to avoid the use of clustering where a Clusters type is needed.

ReducedAlphabet allows the construction of reduced residue alphabets, where residues inside parenthesis belong to the same group.

julia> ab = ReducedAlphabet("(AILMV)(RHK)(NQST)(DE)(FWY)CGP")
MIToS.MSA.ReducedAlphabet of length 8 : "(AILMV)(RHK)(NQST)(DE)(FWY)CGP"

julia> ab[Residue('K')]
2

Most of the MIToS design is created around the Residue bitstype. It has representations for the 20 natural amino acids, a value representing insertions and deletions (GAP, '-') and one representing unknown, ambiguous and non standard residues (XAA, 'X'). Each Residue is encoded as an integer number, with the same bit representation and size than a Int. This allows fast indexing operation of probability or frequency matrices.

Residue creation and conversion

Creation and conversion of Residues should be treated carefully. Residue is encoded as a 32 or 64 bits type similar to Int, to get fast indexing using Int(x::Residue). Int simply calls reinterpret without checking if the residue is valid. Valid residues have integer values in the closed interval [1,22]. convert from Int and Char always returns valid residues, however it's possible to find invalid residues (they are shown using the character '�') after the creation of uninitialized Residue arrays (i.e. using Array). You can use zeros, ones or rand to get initialized Residue arrays with valid residues. Conversions to and from Chars changes the bit representation and allows the use of the usual character representation of residues and amino acids. This conversions are used in IO operations and always return valid residues. In conversions from Char, lowercase letters, '*', '-' and '.' are translated to GAP, letters representing the 20 natural amino (ARNDCQEGHILKMFPSTWYV) acids are translated to their corresponding Residue and any other character is translated to XAA. Since lowercase letters and dots are translated to gaps, Pfam MSA insert columns are converted to columns full of gaps.

julia> alanine = Residue('A')
A

julia> Char(alanine)
'A'

julia> for residue in res"ARNDCQEGHILKMFPSTWYV-X"
           println(residue, " ", Int(residue))
       end
A 1
R 2
N 3
D 4
C 5
Q 6
E 7
G 8
H 9
I 10
L 11
K 12
M 13
F 14
P 15
S 16
T 17
W 18
Y 19
V 20
- 21
X 22

Abstract type to define residue alphabet types.

This type defines the usual alphabet of the 20 natural residues, without the gap character.

julia> UngappedAlphabet()
MIToS.MSA.UngappedAlphabet of length 20. Residues : res"ARNDCQEGHILKMFPSTWYV"

GAP is the Residue representation on MIToS for gaps ('-', insertions and deletions). Lowercase residue characters, dots and '*' are encoded as GAP in conversion from Strings and Chars. This Residue constant is encoded as Residue(21).

XAA is the Residue representation for unknown, ambiguous and non standard residues. This Residue constant is encoded as Residue(22).

The MIToS macro @res_str takes a string and returns a Vector of Residues (sequence).

julia> res"MIToS"
5-element Array{MIToS.MSA.Residue,1}:
 M
 I
 T
 -
 S

It chooses from the 20 natural residues (it doesn't generate gaps).

julia> srand(1); # Reseed the random number generator.

julia> rand(Residue)
Y

julia> rand(Residue, 4, 4)
4×4 Array{MIToS.MSA.Residue,2}:
 E  L  K  P
 V  N  A  M
 I  D  A  F
 F  Y  C  P

It's like Base.shuffle. When a Matrix{Residue} is used, you can indicate if the gaps should remain their positions using the last boolean argument. The previous argument should be the dimension to shuffle, 1 for shuffling residues in a sequence (row) or 2 for shuffling residues in a column.

julia> msa = hcat(res"RRE",res"DDK", res"G--")
3×3 Array{MIToS.MSA.Residue,2}:
 R  D  G
 R  D  -
 E  K  -

julia> srand(42);

julia> shuffle(msa, 1, true)
3×3 Array{MIToS.MSA.Residue,2}:
 G  D  R
 D  R  -
 E  K  -

julia> srand(42);

julia> shuffle(msa, 1, false)
3×3 Array{MIToS.MSA.Residue,2}:
 G  D  R
 D  -  R
 -  E  K

It's like shuffle but in-place. When a Matrix{Residue} or a AbstractAlignedObject (sequence or MSA) is used, you can indicate if the gaps should remain their positions using the last boolean argument.

isvalid(res::Residue)

It returns true if the encoded integer is in the closed interval [1,22].

It returns the name of each group. The name is a string with the one letter code of each residue that belong to the group.

julia> ab = ReducedAlphabet("(AILMV)(RHK)(NQST)(DE)(FWY)CGP")
MIToS.MSA.ReducedAlphabet of length 8 : "(AILMV)(RHK)(NQST)(DE)(FWY)CGP"

julia> names(ab)
8-element Array{String,1}:
 "AILMV"
 "RHK"
 "NQST"
 "DE"
 "FWY"
 "C"
 "G"
 "P"

parse(io, format[, output; generatemapping, useidcoordinates, deletefullgaps])

The keyword argument generatemapping (false by default) indicates if the mapping of the sequences ("SeqMap") and columns ("ColMap") and the number of columns in the original MSA ("NCol") should be generated and saved in the annotations. If useidcoordinates is true (default: false) the sequence IDs of the form "ID/start-end" are parsed and used for determining the start and end positions when the mappings are generated. deletefullgaps (true by default) indicates if columns 100% gaps (generally inserts from a HMM) must be removed from the MSA.

Get a vector of assignments, where the i value is the index/number of the cluster to which the i-th sequence is assigned.

Get the number of clusters in a Clusters object.

Creates a new matrix of residues. This function deletes positions/columns of the MSA with gaps in the reference (first) sequence.

It removes positions/columns of the MSA with gaps in the reference (first) sequence.

Annotates on file annotations the modifications realized by MIToS on the MSA. It always returns true, so It can be used in a boolean context.

annotations returns the Annotations of an MSA or aligned sequence.

Fraction of gaps per column/position on the MSA

columnnames(msa)

It returns a Vector{String} with the sequence names/identifiers. If the msa is a Matrix{Residue} this function returns the actual column numbers as strings. Otherwise it returns the column number of the original MSA through the wrapped NamedArray column names.

Initialize an empty PairwiseListMatrix for a pairwise measure in sequence pairs. It uses the sequence names if they are available, otherwise it uses the actual sequence numbers. You can use the positional argument to indicate the number Type (default: Float64), if the PairwiseListMatrix should store the diagonal values on the list (default: false) and a default value for the diagonal (default: NaN).

Coverage of the sequences with respect of the number of positions on the MSA

Deletes all the MIToS annotated modifications

Deletes columns with 100% gaps, this columns are generated by inserts.

filtercolumns!(msa, mask[, annotate::Bool=true])

It allows to filter MSA or aligned sequence columns/positions using a AbstractVector{Bool} mask. Annotations are updated if annotate is true (default).

filtercolumns!(data::Annotations, mask)

It is useful for deleting column annotations (creating a subset in place).

It's similar to filtercolumns! but for an AbstractMatrix{Residue}

filtersequences!(msa, mask[, annotate::Bool=true])

It allows to filter msa sequences using a AbstractVector{Bool} mask (It removes sequnces with false values). AnnotatedMultipleSequenceAlignment annotations are updated if annotate is true (default).

filtersequences!(data::Annotations, ids::Vector{String}, mask::AbstractArray{Bool,1})

It is useful for deleting sequence annotations. ids should be a list of the sequence names and mask should be a logical vector.

It's similar to filtersequences! but for an AbstractMatrix{Residue}

It calculates the fraction of gaps on the Array (alignment, sequence, column, etc.). This function can take an extra dimension argument for calculation of the gap fraction over the given dimension.

This functions deletes/filters sequences and columns/positions on the MSA on the following order:

1. Removes all the columns/position on the MSA with gaps on the reference (first) sequence.

2. Removes all the sequences with a coverage with respect to the number of

columns/positions on the MSA less than a coveragelimit (default to 0.75: sequences with 25% of gaps).

1. Removes all the columns/position on the MSA with more than a gaplimit

(default to 0.5: 50% of gaps).

Creates a new matrix of Residues (MSA) with deleted sequences and columns/positions. The MSA is edited in the following way:

1. Removes all the columns/position on the MSA with gaps on the reference (first) sequence

2. Removes all the sequences with a coverage with respect to the number of

columns/positions on the MSA less than a coveragelimit (default to 0.75: sequences with 25% of gaps)

1. Removes all the columns/position on the MSA with more than a gaplimit

(default to 0.5: 50% of gaps)

getannotcolumn(ann[, feature[,default]])

It returns per column annotation for feature

getannotfile(ann[, feature[,default]])

It returns per file annotation for feature

getannotresidue(ann[, seqname, feature[,default]])

It returns per residue annotation for (seqname, feature)

getannotsequence(ann[, seqname, feature[,default]])

It returns per sequence annotation for (seqname, feature)

It returns a Vector{Int} with the original column number of each column on the actual MSA. The mapping is annotated in the "ColMap" file annotation of an AnnotatedMultipleSequenceAlignment or in the column names of an NamedArray or MultipleSequenceAlignment.

It takes a ResidueAlphabet and returns a dictionary from group name to group position.

julia> ab = ReducedAlphabet("(AILMV)(RHK)(NQST)(DE)(FWY)CGP")
MIToS.MSA.ReducedAlphabet of length 8 : "(AILMV)(RHK)(NQST)(DE)(FWY)CGP"

julia> getnamedict(ab)
DataStructures.OrderedDict{String,Int64} with 8 entries:
  "AILMV" => 1
  "RHK"   => 2
  "NQST"  => 3
  "DE"    => 4
  "FWY"   => 5
  "C"     => 6
  "G"     => 7
  "P"     => 8

getresidues allows you to access the residues stored inside an MSA or aligned sequence as a Matrix{Residue} without annotations nor column/row names.

getresiduesequences returns a Vector{Vector{Residue}} with all the MSA sequences without annotations nor column/sequence names.

getsequence takes an MSA and a sequence number or identifier and returns an aligned sequence object. If the MSA is an AnnotatedMultipleSequenceAlignment, it returns an AnnotatedAlignedSequence with the sequence annotations. From a MultipleSequenceAlignment, It returns an AlignedSequence object. If an Annotations object and a sequence identifier are used, this function returns the annotations related to the sequence.

It returns the sequence coordinates as a Vector{Int} for an MSA sequence. That vector has one element for each MSA column. If the number if 0 in the mapping, there is a gap in that column for that sequence.

getweight(c[, i::Int])

This function returns the weight of the sequence number i. getweight should be defined for any type used for count!/count in order to use his weigths. If i isn't used, this function returns a vector with the weight of each sequence.

Sequence clustering using the Hobohm I method from Hobohm et. al. 1992.

Hobohm, Uwe, et al. "Selection of representative protein data sets." Protein Science 1.3 (1992): 409-417.

Returns the mean of the percent identity between the sequences of a MSA. If the MSA has 300 sequences or less, the mean is exact. If the MSA has more sequences and the exact keyword is false (defualt), 44850 random pairs of sequences are used for the estimation. The number of samples can be changed using the second argument. Use exact=true to perform all the pairwise comparison (the calculation could be slow).

namedmatrix returns the NamedResidueMatrix stored in an MSA or aligned sequence.

ncolumns returns the number of MSA columns or positions.

ncolumns(ann::Annotations) returns the number of columns/residues with annotations. This function returns -1 if there is not annotations per column/residue.

nsequences returns the number of sequences on the MSA.

percentidentity(seq1, seq2, threshold)

Computes quickly if two aligned sequences have a identity value greater than a given threshold value. Returns a boolean value. Positions with gaps in both sequences doesn't count to the length of the sequences. Positions with a XAA in at least one sequence aren't counted.

percentidentity(seq1, seq2)

Calculates the fraction of identities between two aligned sequences. The identity value is calculated as the number of identical characters in the i-th position of both sequences divided by the length of both sequences. Positions with gaps in both sequences doesn't count to the length of the sequences. Positions with a XAA in at least one sequence aren't counted.

percentidentity(msa[, out::Type=Float64])

Calculates the identity between all the sequences on a MSA. You can indicate the output element type with the last optional parameter (Float64 by default). For a MSA with a lot of sequences, you can use Float32 or Flot16 in order to avoid the OutOfMemoryError().

Calculates the similarity percent between two aligned sequences. The 100% is the length of the aligned sequences minus the number of columns with gaps in both sequences and the number of columns with at least one residue outside the alphabet. So, columns with residues outside the alphabet (other than the specially treated GAP) aren't counted to the protein length. Two residues are considered similar if they below to the same group in a ReducedAlphabet. The alphabet (third positional argument) by default is:

ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP")

The first group is composed of the non polar residues (AILMV), the second group is composed of polar residues, the third group are positive residues, the fourth group are negative residues, the fifth group is composed by the aromatic residues (FWY). C, G and P are considered unique residues.

Other residue groups/alphabets:

SMS (Sequence Manipulation Suite) Ident and Sim:

ReducedAlphabet("(GAVLI)(FYW)(ST)(KRH)(DENQ)P(CM)")

Stothard P (2000) The Sequence Manipulation Suite: JavaScript programs for analyzing and formatting protein and DNA sequences. Biotechniques 28:1102-1104.

Bio3D 2.2 seqidentity:

ReducedAlphabet("(GA)(MVLI)(FYW)(ST)(KRH)(DE)(NQ)PC")

Grant, B.J. et al. (2006) Bioinformatics 22, 2695–2696.

Calculates the similarity percent between all the sequences on a MSA. You can indicate the output element type with the out keyword argument (Float64 by default). For an MSA with a lot of sequences, you can use out=Float32 or out=Flot16 in order to avoid the OutOfMemoryError().

Prints MIToS annotated modifications

This function returns the three letter name of the Residue.

julia> residue2three(Residue('G'))
"GLY"

It calculates the fraction of residues (no gaps) on the Array (alignment, sequence, column, etc.). This function can take an extra dimension argument for calculation of the residue fraction over the given dimension

sequencenames(msa)

It returns a Vector{String} with the sequence names/identifiers.

Initialize an empty PairwiseListMatrix for a pairwise measure in column pairs. It uses the column mapping (column number in the input MSA file) if it’s available, otherwise it uses the actual column numbers. You can use the positional argument to indicate the number Type (default: Float64), if the PairwiseListMatrix should store the diagonal values on the list (default: false) and a default value for the diagonal (default: NaN).

setannotcolumn!(ann, feature, annotation)

It stores per column annotation (1 char per column) for feature

setannotfile!(ann, feature, annotation)

It stores per file annotation for feature

setannotresidue!(ann, seqname, feature, annotation)

It stores per residue annotation (1 char per residue) for (seqname, feature)

setannotsequence!(ann, seqname, feature, annotation)

It stores per sequence annotation for (seqname, feature)

It puts the sequence i (name or position) as reference (first sequence) of the MSA. This function swaps the sequences 1 and i.

stringsequence(seq)
stringsequence(msa, i::Int)
stringsequence(msa, id::String)

It returns the selected sequence as a String.

It swaps the sequences on the positions i and j of an MSA. Also it's possible to swap sequences using their sequence names/identifiers when the MSA object as names.

It takes a three letter residue name and returns the corresponding Residue. If the name isn't in the MIToS dictionary, a XAA is returned.

julia> three2residue("ALA")
A

Get sample counts of clusters as a Vector. Each k value is the number of samples assigned to the k-th cluster.