2017-10-11 8 views
-4

Je suis en train de modifier le programme ci-dessous comme suit:La modification du code afin qu'il puisse lire un fichier, et de produire des sorties correspondantes

  1. La première ligne contient le nom de la protéine et le nombre de les lignes suivantes de sortie pour cette protéine (disons N)

  2. Chacune des N lignes suivantes contient une information de correspondance: les emplacements des GBoxes et les correspondances réelles (rappelez-vous qu'il y a des perturbations et des X, c'est-à-dire des jokers, qui sont permis).

Script:

import csv 
import matplotlib.pyplot as plt 
import numpy as np 

# all G boxes 
def match(x,y): 
     mismatch = 0 
     for i in range(len(x)): 
       if (x[i] == 'X' or x[i] == y[i]): 
         pass 
       else: 
         mismatch += 1 
     if(mismatch <= 1): 
       return True 
     else: 
       return False 



def H(protein,x1,x2,x3,x4): 
     pL1=[] 
     pL2=[] 
     pL3=[] 
     pL4=[] 
     L1=[] 
     L2=[] 
     L3=[] 
     L4=[] 


     for i in range(len(protein)-len(x1)): 
       if(match(x1, protein[i:i+len(x1)]) == True): 
#      global L1 
         pL1=pL1 + [i] 
         L1 = L1+[protein[i:i+len(x1)]] 


     for j in range(len(protein)-len(x2)): 
       if (match(x2, protein[j:j+len(x2)]) == True): 
#      global L2 
         pL2=pL2+[j] 
         L2 = L2+[protein[j:j+len(x2)]] 

     for k in range(len(protein)-len(x3)): 
       if (match(x3, protein[k:k+len(x3)]) == True): 
#       global L3 
         pL3=pL3+[k] 
         L3 = L3+[protein[k:k+len(x3)]] 

     for l in range(len(protein)-len(x4)): 
       if (match(x4, protein[l:l+len(x4)]) == True): 
#      global L3 
         pL4=pL4+[l] 
         L4 = L4+[protein[l:l+len(x4)]] 
     candidates = [] 

     for i in range(len(pL1)): 
       for j in range(len(pL2)): 
         for k in range(len(pL3)): 
           for l in range(len(pL4)): 
             if 40 <=pL2[j]-pL1[i] <= 80 and 40 <=pL3[k]- pL2[j] <= 80 and 20 <=pL4[l]- pL3[k] <= 40: 
               a = L1[i],pL1[i] 
               b = L2[j],pL2[j] 
               c = L3[k],pL3[k] 
               d = L4[l],pL4[l] 
               print a,b,c,d 
               candidates.append((a,b,c,d)) 


     offset = 5 
     for i in np.arange((np.array(candidates).transpose()).shape[1]): 
       barchartdata = np.unique(np.array(candidates).transpose()[:,i]) 
       barchartdata = barchartdata.reshape(2, len(barchartdata)/2) 
       print barchartdata 
       x_pos = np.arange(barchartdata.size/2) 
       print x_pos 
       print barchartdata[0,:] 
       plt.bar(x_pos + 5 * i, barchartdata[0,:]) 

     plt.show() 
     plt.xticks(x_pos, ('g1','g2','g3','g4')) 
     plt.yticks('Count') 
     plt.show() 

x1 = 'GXXXXGK' 
x2 = 'DXXG' 
x3 = 'NKXD' 
x4 = 'EXSAX' 
#input sequence 
protein = 'MAKGEFIRTKPHVNVGTIGHVDHGKTTLTAALTYVAAAENPNVEVKDYGEIDKAPEERARGITINTAHVEYETAKRHYSHVDCPGHADYIKNMITGAAQMDGAILVVSAADGPMPQTREHILLARQVGVPYIVVFMNKVDMVDDPELLDLVEMEVRDLLNQYEFPGDEVPVIRGSALLALEQMHRNPKTRRGENEWVDKIWELLDAIDEYIPTPVRDVDKPFLMPVEDVFTITGRGTVATGRIERGKVKVGDEVEIVGLAPETRKTVVTGVEMHRKTLQEGIAGDNVGVLLRGVSREEVERGQVLAKPGSITPHTKFEASVYVLKKEEGGRHTGFFSGYRPQFYFRTTDVTGVVQLPPGVEMVMPGDNVTFTVELIKPVALEEGLRFAIREGGRTVGAGVVTKILE' 
H(protein,x1,x2,x3,x4) 

Modifier

sortie précédente (My Script) - Correct:

('GAGGVGK', 9) ('DILD', 53) ('NKCD', 115) ('ETSAK', 142) 
('GAGGVGK', 9) ('DTAG', 56) ('NKCD', 115) ('ETSAK', 142) 
('GAGGVGK', 9) ('DQYM', 68) ('NKCD', 115) ('ETSAK', 142) 
('GAGGVGK', 9) ('MRTG', 71) ('NKCD', 115) ('ETSAK', 142) 
('GAGGVGK', 9) ('TGEG', 73) ('NKCD', 115) ('ETSAK', 142) 

Obtention d'une sortie dans votre script:

((17, 'GAGGVGK'), (61, 'DILD'), (123, 'NKCD'), (150, 'ETSAK')) 
((17, 'GAGGVGK'), (64, 'DTAG'), (123, 'NKCD'), (150, 'ETSAK')) 
((17, 'GAGGVGK'), (76, 'DQYM'), (123, 'NKCD'), (150, 'ETSAK')) 
((17, 'GAGGVGK'), (79, 'MRTG'), (123, 'NKCD'), (150, 'ETSAK')) 
((17, 'GAGGVGK'), (81, 'TGEG'), (123, 'NKCD'), (150, 'ETSAK')) 

La longueur est incorrecte. J'ai besoin de courir pour plusieurs séquences, mais il ne fonctionne que pour une séquence. S'il vous plaît me guider

J'essaye aussi de tracer un graphique mais il n'est pas capable d'obtenir la sortie attendue.

graphique attendu:

Dans cette image est le graphique attendu - il faut calculer sage la colonne en pourcentage - s'il vous plaît vérifier « Sortie précédente (My Script) - Correct: s'il vous plaît vérifier ci-dessous l'image par exemple.

Expected graph image

Modifier 1

fichier d'entrée est un fichier CSV a le format suivant (plusieurs lignes):

PDB ID Macromolecule Name Sequence Source 
121P H-RAS P21 PROTEIN MTEYKLVVVGAGGVGKSALTIQLIQNHFVDEYDPTIEDSYRKQVVIDGETCLLDILDTAGQEEYSAMRDQYMRTGEGFLCVFAINNTKSFEDIHQYREQIKRVKDSDDVPMVLVGNKCDLAARTVESRQAQDLARSYGIPYIETSAKTRQGVEDAFYTLVREIRQH Homo sapiens 
1A12 REGULATOR OF CHROMOSOME CONDENSATION 1 RRSPPADAIPKSKKVKVSHRSHSTEPGLVLTLGQGDVGQLGLGENVMERKKPALVSIPEDVVQAEAGGMHTVCLSKSGQVYSFGCNDEGALGRDTSVEGSEMVPGKVELQEKVVQVSAGDSHTAALTDDGRVFLWGSFRDNNGVIGLLEPMKKSMVPVQVQLDVPVVKVASGNDHLVMLTADGDLYTLGCGEQGQLGRVPELFANRGGRQGLERLLVPKCVMLKSRGSRGHVRFQDAFCGAYFTFAISHEGHVYGFGLSNYHQLGTPGTESCFIPQNLTSFKNSTKSWVGFSGGQHHTVCMDSEGKAYSLGRAEYGRLGLGEGAEEKSIPTLISRLPAVSSVACGASVGYAVTKDGRVFAWGMGTNYQLGTGQDEDAWSPVEMMGKQLENRVVLSVSSGGQHTVLLVKDKEQS Homo sapiens 
1A2B TRANSFORMING PROTEIN RHOA SMAAIRKKLVIVGDVACGKTCLLIVFSKDQFPEVYVPTVFENYVADIEVDGKQVELALWDTAGQEDYDRLRPLSYPDTDVILMCFSIDSPDSLENIPEKWTPEVKHFCPNVPIILVGNKKDLRNDEHTRRELAKMKQEPVKPEEGRDMANRIGAFGYMECSAKTKDGVREVFEMATRAALQA Homo sapiens 
1A2K NUCLEAR TRANSPORT FACTOR 2 MGDKPIWEQIGSSFIQHYYQLFDNDRTQLGAIYIDASCLTWEGQQFQGKAAIVEKLSSLPFQKIQHSITAQDHQPTPDSCIISMVVGQLKADEDPIMGFHQMFLLKNINDAWVCTNDMFRLALHNFG Rattus norvegicus 

Répondre

1

Votre code est très difficile à manier. Je l'ai réécrit pour le nettoyer un peu. J'ai enlevé toutes les informations de base sur ce que j'ai choisi de changer et pourquoi, car il y en avait simplement trop. Cependant, voici quelques-unes des choses que je fixe:

  • Vous avez tendance à écrire des boucles unpythonic dans le style

    for i in range(len(x)): 
        do_something_with(x[i]) 
    

    En Python, il est préférable d'écrire

    for element in x: 
        do_something_with(element) 
    
  • Vous avez eu beaucoup de parenthèses superflues autour de if -conditions

  • Vous avez réécrit essentiellement le même gros for -loop quatre fois plutôt que d'écrire une fonction contenant la boucle et de la réutiliser.

Malheureusement, le 20 qui brise le modèle établi par les 40 s dans la for -loop quadruplement-imbriquée, fait le remplacer par une seule boucle sur itertools.product un peu plus de mal que ça valait la peine.

Quoi qu'il en soit, voici la version nettoyée du code.

import csv 
import matplotlib.pyplot as plt 
import numpy as np 

def match(X,Y): 
     mismatch = 0 
     for x,y in zip(X,Y): 
       if not (x == 'X' or x == y): 
         mismatch += 1 
         if mismatch > 1: 
          return False 
     return True 

def H(protein,x1,x2,x3,x4): 

     def find_matches(x): 
      match_positions = [] 
      matches   = [] 
      for i in range(len(protein) - len(x)): 
       candidate = protein[i : i + len(x)] 
       if match(x, candidate): 
        match_positions.append(i) 
        matches  .append(candidate) 
      return matches, match_positions 

     L1, pL1 = find_matches(x1) 
     L2, pL2 = find_matches(x2) 
     L3, pL3 = find_matches(x3) 
     L4, pL4 = find_matches(x4) 

     candidates = [] 
     for a in zip(pL1, L1): 
      for b in zip(pL2, L2): 
       for c in zip(pL3, L3): 
        for d in zip(pL4, L4): 
         if (40 <= b[0] - a[0] <= 80 and 
          40 <= c[0] - b[0] <= 80 and 
          20 <= d[0] - c[0] <= 80 ): 
          print(a,b,c,d) 
          candidates.append((a,b,c,d)) 

     for i in np.arange((np.array(candidates).transpose()).shape[1]): 
       barchartdata = np.unique(np.array(candidates).transpose()[:,i]) 
       barchartdata = barchartdata.reshape(2, len(barchartdata)//2) 
       print (barchartdata) 
       x_pos = np.arange(barchartdata.size/2) 
       print (x_pos) 
       print (barchartdata[0,:]) 
       plt.bar(x_pos + 5 * i, barchartdata[0,:]) 

     plt.show() 
     plt.xticks(x_pos, ('g1','g2','g3','g4')) 
     plt.yticks('Count') 
     plt.show() 

x1 = 'GXXXXGK' 
x2 = 'DXXG' 
x3 = 'NKXD' 
x4 = 'EXSAX' 

protein = 'MAKGEFIRTKPHVNVGTIGHVDHGKTTLTAALTYVAAAENPNVEVKDYGEIDKAPEERARGITINTAHVEYETAKRHYSHVDCPGHADYIKNMITGAAQMDGAILVVSAADGPMPQTREHILLARQVGVPYIVVFMNKVDMVDDPELLDLVEMEVRDLLNQYEFPGDEVPVIRGSALLALEQMHRNPKTRRGENEWVDKIWELLDAIDEYIPTPVRDVDKPFLMPVEDVFTITGRGTVATGRIERGKVKVGDEVEIVGLAPETRKTVVTGVEMHRKTLQEGIAGDNVGVLLRGVSREEVERGQVLAKPGSITPHTKFEASVYVLKKEEGGRHTGFFSGYRPQFYFRTTDVTGVVQLPPGVEMVMPGDNVTFTVELIKPVALEEGLRFAIREGGRTVGAGVVTKILE' 
H(protein,x1,x2,x3,x4) 

Lorsque j'exécute ce code ci-dessus je reçois la sortie

((3, 'GEFIRTK'), (57, 'RARG'), (136, 'NKVD'), (172, 'RGSAL')) 
((3, 'GEFIRTK'), (81, 'DCPG'), (136, 'NKVD'), (172, 'RGSAL')) 
((18, 'GHVDHGK'), (81, 'DCPG'), (136, 'NKVD'), (172, 'RGSAL')) 
((18, 'GHVDHGK'), (87, 'DYIK'), (136, 'NKVD'), (172, 'RGSAL')) 
((18, 'GHVDHGK'), (92, 'MITG'), (136, 'NKVD'), (172, 'RGSAL')) 

Lorsque j'exécute le code qui est visible dans votre question au moment de l'écriture, je reçois la sortie

('GEFIRTK', 3) ('RARG', 57) ('NKVD', 136) ('RGSAL', 172) 
('GEFIRTK', 3) ('DCPG', 81) ('NKVD', 136) ('RGSAL', 172) 
('GHVDHGK', 18) ('DCPG', 81) ('NKVD', 136) ('RGSAL', 172) 
('GHVDHGK', 18) ('DYIK', 87) ('NKVD', 136) ('RGSAL', 172) 
('GHVDHGK', 18) ('MITG', 92) ('NKVD', 136) ('RGSAL', 172) 

Etes-vous absolument sûr que vous avez exécuté les deux codes sur la sa moi les données? Notez que votre code utilise le protein qui est câblé dans le code, alors que mon code était en utilisant les données que vous avez dans le contenu du fichier CSV exemple. Le premier commence par 'MAKG', le second commence par 'MTEY'. Ce ne sont pas les mêmes données, donc vous vous attendez à ce qu'ils donnent des résultats différents!

Notez que la version modifiée du code que je l'ai montré ci-dessus est revenu à l'aide du protein dur câblé il le code lui-même, et donne des résultats identiques à la vôtre

Voici un exemple de code qui montre comment pour lire vos données de protéines à partir du fichier CSV:

with open('xxx.csv') as infile: 
    lines = csv.reader(infile, delimiter=' ', skipinitialspace=True) 
    next(lines) # skip header 
    for line in lines: 
     protein = line[4] 
     # Do whatever you want with protein here ... 
+0

Merci, Pouvez-vous s'il vous plaît me suggérer comment ajouter un fichier CSV en tant que séquence d'entrée? Dans le fichier d'entrée, il n'y a pas de séquence de protéines. Je veux courir pour toute la séquence – vishnu

+0

@vishnu J'ai édité la réponse, pour inclure une réponse à ce que j'imagine que votre question pourrait signifier. Bonne chance. – jacg

+0

Merci beaucoup, mais pas de résultat. J'ai édité ma question s'il vous plaît vérifier une fois – vishnu