Restricting match output for multiple sequence alignment using mafft?

Question

So I aligned roughly 5k sequences and I got my output using mafft. However, I want to restrict the output to only present conserved segments of a specific length (let's say 25bp). Does mafft have parameters to set this that I seem to be missing?

Answer 1

Q - Is there some parameter for this?

No, there is no parameter for this. Getting the conserved regions is something that you will have to perform post-alignment.

How to get the conserved regions?

You will need to define what exactly you mean by conservation and calculate a conservation score for each position of the alignment.

For instance you can look at each column of the alignment and calculate the percent of queries that had the same base/AA at that position. If using AAs you can define conservation more loosely depending on what properties are important to you, or using a scoring matrix.

Once you have a co-occurrence matrix you can define the regions of interest based on a conservation threshold that you define.

Is there software to find conserved columns?

There is a program MstatX which outputs conservation scores of individual columns.

How MstatX works

The example of MstatX provides this example alignment. Each column is progressively less conserved than the last. The last column is conserved, but has some gaps -- this lowers the conservation score.

   col          1
 seq1  1234567890
 seq2  DDDDDDIPDL
 seq3  DDDDDDIPVL
 seq4  DDDDDDIPYL
 seq5  DDDDDDIPA-
 seq6  DDDDDDLWT-
 seq7  DDEDEELWK-
 seq8  DDEDEELWP-
 seq9  DDEDEFVSR-
 seq10 DEEFFFVSH-

The conservation scores for each column were:

1       2.71272e-08
2       0.243433
3       0.312789
4       0.243433
5       0.483364
6       0.480391
7       0.498757
8       0.498757
9       0.973514
10      0.267319

How to get conserved regions

The python program below will output all windows of conserved regions when using a fasta alignment and the default output from MstatX.

• Positional args
  • 1. MstatX_alignment_example.mali, is the fasta alignment file.
  • 2. output.txt is the default output of mstatx when called with mstatx -i MstatX_alignment_example.mali.
  • 3. 0.4 is the minimum conservation score
  • 4. 3 is the minimum number of bases that must have that conservation score.

Calling this program with the following command on the above example:

python print_conserved_regions.py MstatX_alignment_example.mali output.txt 0.4 3 prints:

{1,4}

Meaning that there was only one window at least 3 bases long in which all positions had a conservation score below 0.4. This window was positions 1 through 4 in the alignment.

#call this file cons_windows.py
import sys
from Bio import AlignIO

align_filename = sys.argv[1]
cons_filename = sys.argv[2]
cutoff = float(sys.argv[3])
minchar = int(sys.argv[4])

def parse_conservation(filename):
    cons = {}
    with open(filename, "r") as f:
        for line in f:
            if line.strip():
                temp = line.split()
                cons[int(temp[0])] = float(temp[1])
    return cons

align = AlignIO.read(align_filename, "fasta")
cons = parse_conservation(sys.argv[2])

i = 1
while i < max(cons):
    if cons[i] <= cutoff:
        j = i
        done = False
        while j < max(cons) and not done:
            if (cons[j] > cutoff) and (j-i >= minchar):
                print("{{{},{}}}".format(i,j-1))
                done = True
                i = j
            j+=1
    i+=1