Extraction de la séquence et de l'en-tête du fichier fasta par des séquences connues

Question

J'essaie de comparer deux fichiers et d'extraire les séquences qui ont le sous-ensemble des autres. Et, je veux extraire les identifiants aussi. Cependant, ce que je peux faire, c'est pouvoir extraire les séquences, y compris les sous-ensembles. Les exemples de fichiers sont:

text.fa 
>header1 
ETTTHAASCISATTVQEQ*TLFRLLP 
>header2 
SKSPCSDSDY**AAA 
>header3 
SSGAVAAAPTTA 

et,

textref.fa 
>textref.fa 
CISA 
AAAP 
AATP 

Quand je lance le code, j'ai cette sortie:

ETTTHAASCISATTVQEQ*TLFRLLP 
SSGAVAAAPTTA 

Cependant, mon résultat attendu est avec les en-têtes:

>header1 
ETTTHAASCISATTVQEQ*TLFRLLP 
>header3 
SSGAVAAAPTTA 

Mon code est en deux parties, tout d'abord je crée le fichier avec ces séquences, puis j'essaie de les extraire avec leurs têtes à partir du fichier FASTA d'origine:

def get_nucl(filename): 
    with open(filename,'r') as fd: 
     nucl = [] 
     for line in fd: 
      if line[0]!='>': 
       nucl.append(line.strip()) 
     return nucl 
def finding(filename,reffile): 
     nucl = get_nucl(filename) 
     with open(reffile,'r') as reffile2: 
      for line in reffile2: 
       for element in nucl: 
        if line.strip() in element: 
          yield(element) 



    with open('sequencesmatched.txt','w') as output: 
      results = finding('text.fa','textref.fa',) 
      for res in results: 
       print(res) 
       output.write(res + '\n') 

Alors, je suis dans ce sequencesmatched.txt, ayant les séquences de text.fa qui ont sous-chaînes de textref.fa. comme:

ETTTHAASCISATTVQEQ*TLFRLLP 
SSGAVAAAPTTA 

Ainsi, dans l'autre partie, pour extraire les en-têtes respectifs, et ces séquences:

def finding(filename,seqfile): 
     with open(filename,'r') as fastafile: 
       with open(seqfile,'r') as sequf: 
         alls=[] 
         for line in fastafile: 
           alls.append(line.strip()) 
         print(alls) 
         sequfs = [] 
         for line2 in sequf: 
           sequfs.append(line2.strip()) 
           if str(line.strip()) == str(line2.strip()): 
             num = alls.index(line.strip()) 
             print(alls[num-1] + line) 


print(finding('text.fa','sequencesmatched.txt')) 

Cependant, en tant que sortie, je ne peux extraire une séquence, qui est la première partie:

>header1 
ETTTHAASCISATTVQEQ*TLFRLLP 

Peut-être que je pouvais le faire sans le second fichier, mais je ne pouvais pas faire le droit des boucles pour obtenir les séquences et leurs têtes respectives. Par conséquent, je suis allé au long chemin ..

Je serais heureux si vous pouvez aider!

Answer 1

Vous pouvez faire quelque chose beaucoup plus facile si votre fichier est toujours la même structure:

def get_nucl(filename): 
    with open(filename, 'r') as fd: 
     headers = {} 
     key = '' 
     for line in fd.readlines():  
      if '>' in line: 
       key = line.strip()[1:] # to remove the '>' 
      else: 
       headers[key] = line.strip() 

    return headers 

Ici, je suis en supposant que votre fichier commence par « > headern » que ce soit, si vous avez de ne pas ajouter un peu de test. Maintenant, vous avez un dictionnaire comme headers['header1'] = 'ETTTHAASCISATTVQEQ*TLFRLLP'.

Alors maintenant, pour trouver les matches que vous utilisez juste que dict:

def finding(filename, reffile): 
    headers = get_nucl(filename) 
    with open(reffile, 'r') as f: 
     matches = {} 
     for line in f.readlines(): 
      for key, value in headers.items(): 
       if line.stip() in value and key not in matches: 
        matches[key] = value 

    return matches 

Alors que vous avez un dict avec en-têtes correspondant à leurs valeurs, il vous suffit de vérifier dans le dict si vous avez une chaîne de sous, et vous avez déjà la valeur de l'en-tête en tant que clé.

Juste vu que vous avez fait print(finding(....), votre fonction est déjà imprimée, il suffit donc de l'appeler.