How to build an empirical codon substitution matrix from a multiple sequence alignment

Question

I have been trying to build an empirical codon substitution matrix given a multiple sequence alignment in fasta format using Biopython.

It appears to be relatively straigh-forward for single nucleotide substitution matrices using the AlignInfo module when the aligned sequences have the same length. Here is what I managed to do using python2.7:

#!/usr/bin/env python 
import os
import argparse
from Bio import AlignIO
from Bio.Align import AlignInfo
from Bio import SubsMat
import sys

version = "0.0.1 (23.04.20)"
name = "Aln2SubMatrix.py"
parser=argparse.ArgumentParser(description="Outputs a codon substitution matrix given a multi-alignment in FastaFormat. Will raise error if alignments contain dots (\".\"), so replace those with dashes (\"-\") beforehand (e.g. using sed)")
parser.add_argument('-i','--input', action = "store", dest = "input", required = True, help = "(aligned) input fasta")
parser.add_argument('-o','--output', action = "store", dest = "output", help = "Output filename (default = <Input-file>.codonSubmatrix")
args=parser.parse_args()
if not args.output:
    args.output = args.input + ".codonSubmatrix"  #if no outputname was specified set outputname based on inputname

def main():
    infile = open(args.input, "r")
    outfile = open(args.output, "w")
    align = AlignIO.read(infile, "fasta")
    summary_align = AlignInfo.SummaryInfo(align)
    replace_info = summary_align.replacement_dictionary()
    mat = SubsMat.SeqMat(replace_info)
    print >> outfile, mat
    infile.close()
    outfile.close()
    sys.stderr.write("\nfinished\n")

main()

Using a multiple sequence alignment file in fasta format with sequences of same length (aln.fa), the output is a half-matrix corresponding to the number of nucleotide substitutions oberved in the alignment (Note that gaps (-) are allowed):

python Aln2SubMatrix.py -i aln.fa

-   0
a 860 232
c 596  75 129
g 571 186  75 173
t 892  58 146  59 141
   -   a   c   g   t

What I am aiming to do is to compute similar empirical substitution matrix but for all nucleotide triplets (codons) present in a multiple sequence alignment.

I have tried to tweak the _pair_replacement function of the AlignInfo module in order to accept nucleotide triplets by changing:

line 305 to 308

for residue_num in range(len(seq1)):
            residue1 = seq1[residue_num]
            try:
                residue2 = seq2[residue_num]

to

for residue_num in range(0, len(seq1), 3):
            residue1 = seq1[residue_num:residue_num+3]
            try:
                residue2 = seq2[residue_num:residue_num+3]

At this stage it can retrieve the codons from the alignment but complains about the alphabet (the module only accepts single character alphabet?).

Note that

(i) I would like to get a substitution matrix that accounts for the three possible reading frames

Any help is highly appreciated.

How to build an empirical codon substitution matrix from a multiple sequence alignment

0 Answers0