J'écris un programme qui analyse les allèles de séquence. J'ai écrit du code qui lit un fichier et crée un tableau d'en-têtes et un tableau de séquences. Voici un exemple de fichier:

>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

Les en-têtes sont ('> DQB1', '> GENDX' et '> OMIXON') et les trois séquences sont les trois autres chaînes comme vu ci-dessus.

La partie suivante de mon code détecte si la séquence d'allèles est complète ou incomplète. Un allèle est déterminé comme "incomplet" s'il y a plus de 4 ruptures dans la séquence> DQB1. (Une pause est signifiée par un «-»). Par exemple, la séquence ci-dessus est interrompue car il y a cinq ruptures.

J'essaie d'écrire du code où s'il y a un allèle incomplet détecté, le programme crée un nouveau tableau avec uniquement les en-têtes et les séquences> GENDX et> OMIXON.

Comment puis-je créer un tableau qui ne comprend pas> DQB1?

Voici mon code tel quel:

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

Le problème avec le code ci-dessus est que chaque fois que j'imprime header2, il inclut toujours le DQB1.

1
Gia Constantina 8 août 2016 à 19:29

3 réponses

Meilleure réponse

Pourquoi voulez-vous utiliser re.finditer?

Qu'en est-il de

if header[m].find('OMIXON') > -1 or header[m].find('GENDX') > -1:
2
Michael Westwort 8 août 2016 à 16:56

Vous pouvez le faire très facilement en utilisant Biopython, en supposant que vos séquences sont enregistrées dans un fichier FASTA.

from Bio import SeqIO

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

Et tu as fini.

Si vous souhaitez obtenir la partie séquence de l'objet parse(), utilisez simplement record.seq.

0
MattDMo 8 août 2016 à 17:09

La fonction re.finditer ne fait pas ce que vous pensez qu'elle fait. Pour un exemple, voir ici.

Je recommanderais d'utiliser ceci à la place:

if header[m][1:7] == 'OMIXON' or header[m][1:6]=='GENDX':
0
Noam 8 août 2016 à 17:12