python awk 的替代品？

Question

我有两个fasta文件，我想搜索序列ID，只将ID对应的序列赋值给Python中的一个字符串。

我目前有：

import os
#use awk on the command line to search reference file and cut the reference sequence
os.system("awk '/LOC_OS05G45410.1/{getline;print}' Ref_seqs.fasta > sangerRef") 

#use awk on the command line to cut the aligned sequence
os.system("awk '/seq1/{getline;print}' Sanger_seq_1.fasta > sangerAlign")

Ref_seq = open('sangerRef', 'r').read()
Sanger_seq = open('sangerAlign', 'r').read()

当我打印这些变量时，一切看起来都很好：

TGGTGAGGCTTTTGACAGGGTTGAGCTGAGCCTGGTCTCCCTGGAGAAACTCTTCCAGAGAGCAAATGATGCTTGCACAGCTGCTGAAGAAATGTACTCCCATGGTCATGGTGGTACTGAACCCAG

CTGCTGCCCAAGTACTTCAAGCACAACAACTTCTCCAGCTTCATCAGGCAGCTCAACGCCTACGGTTTCCGAAAAATCGATCCTGAGAGATGGGAGTTCGCAAACGAGGATTTCATAAGAGGGCACACGCACCTT

但是，当我尝试将这些变量读入另一个函数时，它不起作用:

from Bio import pairwise2
from Bio.Align import substitution_matrices

#load sequences
s1=Ref_seq
s2=Sanger_seq
matrix = substitution_matrices.load("NUC.4.4")
gap_open = -10
gap_extend = -0.5

align = pairwise2.align.globalds(s1, s2, matrix, gap_open, gap_extend)

align

我认为用 Python 命令替换 awk 命令可能会更好？

Answer 1

我认为这是因为你没有解析序列。不过，我不知道我是否用对了 'Parse' 这个词。

我认为这应该可行

from Bio import SeqIO

s1 = SeqIO.read('filepath/filename.fasta','fasta')
s2 = SeqIO.read('filepath/file.fasta','fasta')

matrix = substitution_matrices.load("NUC.4.4")

gap_open = -10
gap_extend = -0.5

align = pairwise2.align.globalds(s1.seq, s2.seq, matrix, gap_open, gap_extend)
align

Answer 2

眼前的问题是 read() returns 所有行的末尾都有一个换行符。

但实际上，您的 Awk 命令应该很容易用原生命令替换 Python。

def getseq(filename, search):
    with open(filename) as reffile:
        for line in reffile:
            if search in line:
                return seqfile.__next__().rstrip('\n')

s1 = getseq("Ref_seqs.fasta", "LOC_OS05G45410.1")
s2 = getseq("Sanger_seq_1.fasta", "seq1")

可能 BioPython 已经包含一个更好的函数来执行此操作。特别是，您的 Awk 脚本（以及因此盲目重新实现）假定每个序列仅占用文件中的一行。

python awk 的替代品？

python alternative for awk?

python

awk

bioinformatics