1

我一直在尝试使用 Biopython 建立一个经验密码子替换矩阵,给出一个 fasta 格式的多序列比对。

当比对的序列具有相同的长度时,使用AlignInfo模块的单核苷酸替换矩阵似乎相对简单。这是我使用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()

使用具有相同长度的序列 (aln.fa) 的 fasta 格式的多序列比对文件,输出是一个半矩阵,对应于比对中观察到的核苷酸替换的数量(请注意,允许有空位 (-)):

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

我的目标是计算类似的经验替换矩阵,但对于多序列比对中存在的所有核苷酸三联体(密码子)。

我试图调整 AlignInfo 模块的_pair_replacement函数,以便通过更改来接受核苷酸三联体:

第 305 至 308 行

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

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

在这个阶段,它可以从比对中检索密码子,但抱怨字母表(模块只接受单字符字母表?)。

注意

(i) 我想得到一个替代矩阵,它解释了三个可能的阅读框架

非常感谢任何帮助。

4

0 回答 0