2016-08-08 17 views
1

Ich schreibe ein Programm, das Sequenz-Allele analysiert. Ich habe Code geschrieben, der eine Datei liest und ein Header-Array und ein Sequenz-Array erstellt. Hier ist ein Beispiel für eine Datei:Erstellen Sie ein neues Array

>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 

Die Header sind (‚> DQB1‘, ‚> GENDX‘ und ‚> OMIXON‘) und die drei Sequenzen sind die anderen drei Saiten wie oben zu sehen.

Der nächste Teil meines Codes erkennt, ob die Allelsequenz vollständig oder unvollständig ist. Ein Allel wird als "unvollständig" bestimmt, wenn es mehr als 4 Unterbrechungen innerhalb der> DQB1-Sequenz gibt. (Eine Pause wird durch ein '-' bezeichnet). Zum Beispiel ist die obige Sequenz unterbrochen, weil es fünf Unterbrechungen gibt.

Ich versuche Code zu schreiben, wo, wenn ein unvollständiges Allel entdeckt wird, das Programm ein neues Array mit nur den> GENDX und den> OMIXON Headern und Sequenzen erstellt.

Wie kann ich ein Array erstellen, das> DQB1 nicht enthält?

Hier ist mein Code wie:

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 

Das Problem mit dem obigen Code ist, dass, wenn ich drucken header2 es enthält noch die DQB1.

Antwort

2

Warum möchten Sie re.finditer verwenden?

Was

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

Sie dies sehr leicht Biopython mit tun können, vorausgesetzt, Ihre Sequenzen gespeichert werden in einer FASTA-Datei.

from Bio import SeqIO 

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

und fertig. Wenn Sie den Sequenzteil vom Objekt parse() erhalten möchten, verwenden Sie einfach record.seq.

+0

Ursprünglich schrieb ich den Code mit biopython und es wurde eine Menge Probleme haben, so ich herausgefunden, wie es zu tun, ohne. –

+0

@GiaConstantina, welche Probleme hatten Sie? – MattDMo

+0

Ich glaube nicht, dass es richtig heruntergeladen wurde bc im Allgemeinen würde es mich nicht auf Dinge wie SeqIO und andere Bibliotheken zugreifen lassen. –

0

Die re.finditer Funktion tut nicht, was Sie denken, dass es tut. Ein Beispiel finden Sie unter here.

Ich würde empfehlen, diese stattdessen mit:

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