基于对齐修剪序列答案

【问题标题】：Trim sequences based on alignment基于对齐修剪序列
【发布时间】：2019-01-08 16:14:18
【问题描述】：

我正在尝试使用 BioPython 编辑 ClustalW 生成的 MSA（多序列比对）文件，以修剪一致序列之前的序列。 xxx 指的是这里不相关的其他碱基

这是 I/O 示例：

输入

ITS_primer_fw               --------------------------------CGCGTCCACTMTCCAGTT
RBL67ITS_full_sequence      CCACCCCAACAAGGGCGGCCACGCGGTCCGCTCGCGTCCACTCTCCAGTTxxxxxxxxxxxxxxxxxxxxxxx
PRL2010                     ACACCCCCGAAAGGGCGTCC------CCTGCTCGCGTCCACTATCCAGTTxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
BBF32_3                     ACACACCCACAAGGGCGAGCAGGCG----GCTCGCGTCCACTATCCAGTTxxxxxxxxxxxxxx
BBFCG32                     CAACACCACACCGGGCGAGCGGG-------CTCGCGTCCACTGTCGAGTTxxxxxxxxxxxxxxxxxxxxxx

预期输出

ITS_primer_fw               CGCGTCCACTMTCCAGTT
RBL67ITS_full_sequence      CGCGTCCACTCTCCAGTTxxxxxxxxxxxxxxxxxxxx
PRL2010                     CGCGTCCACTATCCAGTTxxxxxxxxxxxxxxxxxxxxx
BBF32_3                     CGCGTCCACTATCCAGTTxxxxxxxxxxxxxxxxxxx
BBFCG32                     CGCGTCCACTGTCGAGTTxxxxxxxxxxxxxxxxxxxx

AlignIO 的文档化代码仅描述了一种通过将比对视为数组来提取序列的方法。在这个例子中

align = AlignIO.read(input_file, "clustal")
sub_alignment = align[:,20:]

我能够提取从第 20 个核苷酸开始的所有序列 (:) 的子比对。我正在寻找一种方法将示例中的20 替换为共有序列的第一个核苷酸的位置。

【问题讨论】：

添加了正则表达式标签，可能作为 biopython 的替代解决方案。
输入文件总是3行吗？
输入文件有多行，这里以三行为例。这怎么能用正则表达式来完成？我正在学习python，所以我目前不太流利。
我不知道python，认为这可以使用正则表达式来完成。如果您使用 R 添加更多示例行，我可以尝试一下？
我在 I/O 中添加了另外两行。使用 R 不是首选方法，因为我在不使用 R 的 python 3.x 控制台中运行它。这就是为什么我不使用像 awk 这样的 bash 文本操作可能性来执行此操作的原因。

标签： python regex numpy bioinformatics biopython

【解决方案1】：

感谢 Biostars 用户找到了答案。

抽搐正在查看列以找到起点，这将在最后一个“-”之后如预期的那样。默认情况下，我对齐的第一行是最短的，并且在对齐之前以“-”开头。

代码如下：

aln = AlignIO.read(input_file, "clustal")
for col in range(aln.get_alignment_length()):  # search into column
    res = list(aln[:,col])
    if not '-' in res:
        position = col                         # find start point
        print('First full column is {}'.format(col))
        break
print(aln[:,position::])                       # print the whole alignment starting from the position variable found

【讨论】：

我发现您的回答非常有用！这是反转对齐的代码，因此您也可以修剪另一端。 def reverse_ma(aln): for i in range(len(aln)): aln._records[i].seq = aln._records[i].seq[::-1] return aln