Example

In this simple demonstration, you will see how to calculate ZBLMIp (Z score of the corrected MIp using BLOSUM62 pseudo frequencies) for a Pfam MSA from the Julia REPL or using a MIToS script in the system command line.

MIToS in the Julia REPL

If you load the Pfam module from MIToS, you will get access to a set of functions that work with Pfam MSAs. In this case, we are going to use it for download a Stockholm MSA from the Pfam website and read it into Julia.

using MIToS.Pfam
pfam_file = downloadpfam("PF10660")
msa = read(pfam_file, Stockholm, generatemapping=true, useidcoordinates=true)
AnnotatedMultipleSequenceAlignment with 1256 annotations : 614×64 Named Matrix{MIToS.MSA.Residue}
               Seq ╲ Col │  32   35   36   37   38  …  107  108  109  110  111
─────────────────────────┼────────────────────────────────────────────────────
A0A2I2Y8P5_GORGO/9-41    │   -    -    -    -    -  …    K    D    H    R    N
A0A670IAA5_PODMU/1-34    │   -    -    -    -    -       K    D    R    C    -
A0A158NWR3_ATTCE/11-73   │   M    E    P    I    A       R    -    -    -    -
V9KQ82_CALMI/5-38        │   -    -    -    -    -       K    D    K    -    -
A0A1J1J3N7_9DIPT/1-64    │   M    E    L    I    S       A    S    -    -    -
A0A091D4B4_FUKDA/1-66    │   M    E    S    V    A       K    K    K    Q    Q
A0A0A0AJQ4_CHAVO/1-31    │   -    -    -    -    -       K    K    K    Q    Q
A0A5F8GCI8_MONDO/1-59    │   -    E    S    M    A       -    -    -    -    -
A0A3M6TGR7_9CNID/1-64    │   M    D    A    V    S       S    K    P    -    -
⋮                            ⋮    ⋮    ⋮    ⋮    ⋮  ⋱    ⋮    ⋮    ⋮    ⋮    ⋮
A0A091RJP3_9GRUI/1-31    │   -    -    -    -    -       K    K    K    Q    Q
G1MRN1_MELGA/16-53       │   -    -    -    -    -       K    D    K    C    -
A0A2K6DQQ1_MACNE/2-41    │   -    -    -    -    -       K    D    N    R    -
G3H4L8_CRIGR/2-34        │   -    -    -    -    -       K    E    N    R    -
A0A556V466_BAGYA/1-66    │   -    E    T    I    S       Q    R    K    K    Q
G3SVU7_LOXAF/8-41        │   -    -    -    -    -       K    D    H    R    -
A0A5J5MRP2_MUNRE/1-30    │   -    -    -    -    -       K    D    H    R    -
A0A6H5GFP8_9HEMI/1-64    │   M    E    P    V    H       K    L    -    -    -
A0A194RSG1_PAPMA/1-63    │   M    Y    F    V    S  …    -    -    -    -    -
Note

Generation of sequence and column mappings The keyword argument generatemapping of read allows to generate sequence and column mappings for the MSA. Column mapping is the map between of each column on the MSA object and the column number in the file. Sequence mappings will use the start and end coordinates in the sequence ids for enumerate each residue in the sequence if useidcoordinates is true.

You can plot this MSA and other MIToS’ objects using the Plots package. The installation of Plots is described in the Installation section of this site:

using Plots
gr()
plot(msa)
QStandardPaths: XDG_RUNTIME_DIR not set, defaulting to '/tmp/runtime-juliateam'
QPainter::begin: Paint device returned engine == 0, type: 3
QPainter::setCompositionMode: Painter not active
QWidget::paintEngine: Should no longer be called
QPainter::begin: Paint device returned engine == 0, type: 1

The Information module of MIToS has functions to calculate measures from the Information Theory, such as Entropy and Mutual Information (MI), on a MSA. In this example, we will estimate covariation between columns of the MSA with a corrected MI that use the BLOSUM62 matrix for calculate pseudo frequencies (BLMI).

using MIToS.Information
ZBLMIp, BLMIp = BLMI(msa)
ZBLMIp # shows ZBLMIp scores
63×63 Named PairwiseListMatrices.PairwiseListMatrix{Float64, false, Vector{Float64}}
Col1 ╲ Col2 │           35            36  …           110           111
────────────┼──────────────────────────────────────────────────────────
35          │          NaN    0.00556154  …    -0.0477258     -0.261566
36          │   0.00556154           NaN       -0.0621987     -0.413523
37          │  -0.00114558     -0.031544        -0.222168     -0.305511
38          │    0.0579368    -0.0709433        -0.276483      -0.41669
39          │   -0.0685672     0.0421717       -0.0385996     -0.524873
40          │  0.000717448     0.0591496       -0.0422035     -0.325192
41          │    0.0178404    -0.0515125         0.021301      -0.23127
42          │     0.106721    -0.0967409       -0.0229173     -0.339495
43          │      -0.0342   -0.00752816       -0.0248318     -0.438674
⋮                        ⋮             ⋮  ⋱             ⋮             ⋮
103         │   -0.0366782     0.0111343        0.0572733     0.0991304
104         │     0.176451      0.114633        -0.133585     -0.329759
105         │    0.0437504     0.0835655        -0.111321      0.242734
106         │  -0.00526741      -0.04554        0.0445883      0.121984
107         │     0.150476      0.087743        0.0137499     -0.227666
108         │   -0.0711401    0.00974643        0.0508096    -0.0942058
109         │     0.266546      0.203712        -0.137783      0.157585
110         │   -0.0477258    -0.0621987              NaN      0.188203
111         │    -0.261566     -0.413523  …      0.188203           NaN

Once the Plots package is installed and loaded, you can use its capabilities to visualize this results:

heatmap(ZBLMIp, yflip=true, c=:grays)
┌ Warning: Attribute alias `ratio` detected in the user recipe defined for the signature (::NamedArrays.NamedMatrix{Float64, PairwiseListMatrices.PairwiseListMatrix{Float64, false, Vector{Float64}}, Tuple{OrderedCollections.OrderedDict{String, Int64}, OrderedCollections.OrderedDict{String, Int64}}}). To ensure expected behavior it is recommended to use the default attribute `aspect_ratio`.
└ @ Plots ~/.julia/packages/Plots/SVksJ/src/pipeline.jl:26
QStandardPaths: XDG_RUNTIME_DIR not set, defaulting to '/tmp/runtime-juliateam'
QPainter::begin: Paint device returned engine == 0, type: 3
QPainter::setCompositionMode: Painter not active
QWidget::paintEngine: Should no longer be called
QPainter::begin: Paint device returned engine == 0, type: 1

MIToS in system command line

Calculate ZBLMIp on the system shell is easy using the MIToS script called BLMI.jl. This script reads a MSA file, and writes a file with the same base name of the input but with the .BLMI.csv extension.

BLMI.jl PF14972.stockholm.gz