Скрипт на языке Python


сайт ФББ

сайт МГУ

Cкрипт получает на вход fasta-файл ххх, выводит в консоль информацию о последовательности в fasta-формате, создает файл yyy.txt с таблицей, в которой подсчитана доля в процентах встречаемости аминокислот в последовательности белка/оснований в посдедовательности нуклейновой кислоты.

Чтобы запустить скрипт, нужно скачать его, скачать какой-либо fasta-файл в ту же директорию, что и скрипт, в 7-ой строке скрипта изменить "xxx.fasta" на имя скачанного fasta-файла. Далее зайти в директорию, где лежат скипт и файл fasta, через Far Manager, написать в командной строке: python Gorbacheva_script.py

        # MSU FBB. Bioinformatics. Term 1. Block 2. Practice 10. Task 2.     
        # AAFrequency&FastaReader v0.1                                             
        # Author: Gorbacheva D. V.                                      
        # Last modification date: 24.12.2013                            


        filename = "xxx.fasta"                                                        
        frequency = open("yyy.txt", "w")                                              
        fasta_file = open(filename)                                                   
                                                                                      
        class Protein:                                                                
            ide = ""                                                                  
            description = ""                                                          
            sequence = ""                                                             
            def __init__(self):                                                       
               0                                                                      
            def print_fasta(self, length):                                            
                if length < len(self.sequence):                                       
                   start = 0                                                          
                   print self.sequence[0: length]                                     
                   while start <= len(self.sequence):                                 
                         start = start + length                                       
                         s = self.sequence[start: start + length]                     
                         print s                                                      
                else:                                                                 
                   print self.sequence                                                
                                                                                      
        protein_list = []                                                             
        mydict = {}                                                                   
                                                                                      
        for line in fasta_file:                                                       
            if line[0] == ">":                                                        
               otherprotein = Protein()                                               
               otherprotein.ide = line.split(None, 1)[0].strip(">")                   
               otherprotein.description = line.split(None, 1)[1].strip()              
               otherprotein.sequence = ""                                             
               protein_list.append(otherprotein)                                      
            if len(line)==0:                                                          
               continue                                                               
            if line[0] != ">":                                                        
               protein_list[-1].sequence = protein_list[-1].sequence + line.strip()   
               for symbol in line:                                                    
                   if symbol in mydict:                                               
                      mydict[symbol] = mydict[symbol] + 1                             
                   elif symbol == "\t" or symbol=="\n":                               
                      continue                                                        
                   else:                                                              
                      mydict[symbol] = 1                                              
                                                                                      
        x = 0                                                                         
        proteomseq = protein_list[x].sequence                                         
        while x < len(protein_list)-1:                                                
              x = x + 1                                                               
              proteomseq = proteomseq + protein_list[x].sequence                      
                                                                                      
        def percent(symbol, mydict, proteomseq):                                      
            aapercent = str((float(mydict[symbol]))  / len(proteomseq) * 100)         
            frequency.write(symbol + "\t" + aapercent + "%\n")                        
                                                                                      
        for sbl in mydict.keys():                                                     
            percent(sbl, mydict, proteomseq)                                          
                                                                                      
        length = int(raw_input("String length: "))                                    
        protnumber = range(0 , len(protein_list))                                     
        for n in protnumber:                                                          
            firstline = ">" +  protein_list[n].ide + " " + protein_list[n].description
            print firstline                                                           
            protein_list[n].print_fasta(length)                                       
                                                                                      
        frequency.close()                                                             
        fasta_file.close()

Скачать скрипт

Пример входного fasta-файла

© Дарья Горбачева

изменено 8.08.2014