我一直在尝试使用 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) 我想得到一个替代矩阵,它解释了三个可能的阅读框架
非常感谢任何帮助。