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?
1 Answer
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 ofmstatx
when called withmstatx -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.
- 1.
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