做一个新的数组

Question

我正在写一个分析序列等位基因的程序。我已经编写了读取文件并创建头数组和序列数组的代码。下面是一个文件的示例：

>DQB1*04:02:01
------------------------------------------------------------
--ATGTCTTGGAAGAAGGCTTTGCGGAT-------CCCTGGAGGCCTTCGGGTAGCAACT
GTGACCTT----GATGCTGGCGATGCTGAGCACCCCGGTGGCTGAGGGCAGAGACTCTCC
CGAGGATTTCGTGTTCCAGTTTAAGGGCATGTGCTACTTCACCAACGGGACCGAGCGCGT
GTTGGAGCTCCGCACGACCTTGCAGCGGCGA-----------------------------
---GTGGAGCCCACAGTGACCATCTCCCCATCCAGGACAGAGGCCCTCAACCACCACAAC
CTGCTGGTCTGCTCAGTGACAG----CATTGGAGGCTTCGTGCTGGGGCTGATCTTCCTC
GGGCTGGGCCTTATTATC--------------CATCACAGGAGTCAGAAAGGGCTCCTGC
ACTGA-------------------------------------------------------
>OMIXON_CONSENSUS_M_155_09_4890_DQB1*04:02:01
-------------------ATCAGGTCCAAGCTGTGTTGACTACCACTACTTTTCCCTTC
GTCTCAATTATGTCTTGGAAGAAGGCTTTGCGGATCCCTGGAGGCCTTCGGGTAGCAACT
GTGACCTTGATGCTGGCGATGCTGAGCACCCCGGTGGCTGAGGGCAGAGACTCTCCCGGT
AAGTGCAGGGCCACTGCTCTCCAGAGCCGCCACTCTGGGAACAGGCTCTCCTTGGGCTGG
GGTAGGGGGATGGTGATCTCCATGATCTCGGACACAATCTTTCATCAACATTTCCTCTCT
TTGGGGAAAGAGAACGATGTTGCATTCCCATTTATCTTT---------------------
>GENDX_CONSENSUS_M_155_09_4890_DQB1*04:02:01
TGCCAGGTACATCAGATCCATCAGGTCCAAGCTGTGTTGACTACCACTACTTTTCCCTTC
GTCTCAATTATGTCTTGGAAGAAGGCTTTGCGGATCCCTGGAGGCCTTCGGGTAGCAACT
GTGACCTTGATGCTGGCGATGCTGAGCACCCCGGTGGCTGAGGGCAGAGACTCTCCCGGT
AAGTGCAGGGCCACTGCTCTCCAGAGCCGCCACTCTGGGAACAGGCTCTCCTTGGGCTGG
GGTAGGGGGATGGTGATCTCCATGATCTCGGACACAATCTTTCATCAACATTTCCTCTCT

标题是(‘>DQB1 1’，'>GENDX'，和'>OMIXON')，这三个序列是另外三个字符串，如上文所示.

代码的下一部分将检测等位基因序列是完整的还是不完整的。如果>DQB1 1序列中有4个以上的断裂，则一个等位基因被确定为“不完全”。(休息以“-”表示)。例如，上面的序列被破坏了，因为有五个中断。

我正在编写代码，如果检测到一个不完整的等位基因，程序将创建一个新的数组，其中只有>GENDX和>OMIXON报头和序列。

如何制作不包括>DQB1 1的数组？

下面是我的代码：

import sys, re

max_num_breaks=4
filename=sys.argv[1]
f=open(filename,"r")
header=[]
header2=[]
sequence=[]
sequence2=[]
string=""
for line in f:
    if ">" in line and string=="":
        header.append(line[:-1])
    elif ">" in line and string!="":
        sequence.append(string)
        header.append(line[:-1])
        string=""
    else:
        string=string+line[:-1]
sequence.append(string)
s1=sequence[0]
breaks=sum(1 for m in re.finditer("-+",''.join(s1.splitlines())))
if breaks>max_num_breaks:
    print "Incomplete Reference Allele Detected"
    for m in range(len(header)):
        if re.finditer(header[m], 'OMIXON') or re.finditer(header[m], 'GENDX'):
            header2.append(header[m])
            sequence2.append(sequence[m])
    print header2

上面的代码的问题是，每当我打印header2时，它仍然包括DQB1。

Answer 1

你为什么要使用re.finditer

关于

if header[m].find('OMIXON') > -1 or header[m].find('GENDX') > -1:

Answer 2

假设您的序列保存在FASTA文件中，那么使用Biopython可以很容易地做到这一点。

from Bio import SeqIO

headers = [record.id for record in SeqIO.parse("myfile.fasta", "fasta")][1:]

你就完蛋了。

如果要从parse()对象获取序列部分，只需使用record.seq即可。

Answer 3

re.finditer函数不执行您认为它所做的事情。有关示例，请参见这里。

我建议用这个代替：

if header[m][1:7] == 'OMIXON' or header[m][1:6]=='GENDX':