выравнивание последовательностей

Этот код, надеюсь, то, что вы ищете:

class test():

    def get_seq(self, pdb,fasta_file): # Get sequences
        from Bio.PDB.PDBParser import PDBParser
        from Bio import SeqIO
        aa = {'ARG':'R','HIS':'H','LYS':'K','ASP':'D','GLU':'E','SER':'S','THR':'T','ASN':'N','GLN':'Q','CYS':'C','SEC':'U','GLY':'G','PRO':'P','ALA':'A','ILE':'I','LEU':'L','MET':'M','PHE':'F','TRP':'W','TYR':'Y','VAL':'V'}
        p=PDBParser(PERMISSIVE=1)
        structure_id="%s" % pdb[:-4]
        structure=p.get_structure(structure_id, pdb)
        residues = structure.get_residues()
        seq_pdb = ''
        for res in residues:
            res = res.get_resname() 
            if res in aa:
                seq_pdb = seq_pdb+aa[res]           

        handle = open(fasta_file, "rU")
        for record in SeqIO.parse(handle, "fasta") :
            seq_fasta = record.seq
        handle.close()
        self.seq_aln(seq_pdb,seq_fasta)

    def seq_aln(self,seq1,seq2): # Align the sequences
        from Bio import pairwise2
        from Bio.SubsMat import MatrixInfo as matlist

        matrix = matlist.blosum62
        gap_open = -10
        gap_extend = -0.5

        alns = pairwise2.align.globalds(seq1, seq2, matrix, gap_open, gap_extend)
        top_aln = alns[0]
        aln_seq1, aln_seq2, score, begin, end = top_aln
        with open('aln.fasta', 'w') as outfile:
            outfile.write('> PDB_seq\n'+str(aln_seq1)+'\n> Uniprot_seq\n'+str(aln_seq2))
        print aln_seq1+'\n'+aln_seq2
        self.seq_id('aln.fasta')

    def seq_id(self,aln_fasta): # Get sequence ID
        import string
        from Bio import AlignIO

        input_handle = open("aln.fasta", "rU")
        alignment = AlignIO.read(input_handle, "fasta")
        j=0 # counts positions in first sequence
        i=0 # counts identity hits
        for record in alignment:
            #print record
            for amino_acid in record.seq:
                if amino_acid == '-':
                    pass
                else:
                    if amino_acid == alignment[0].seq[j]:
                        i += 1
                j += 1
            j = 0
            seq = str(record.seq)
            gap_strip = seq.replace('-', '')

            percent = 100*i/len(gap_strip)
            print record.id+' '+str(percent)
            i=0


a = test()
a.get_seq('1DBH.pdb','Q07889.fasta')

Это выводит:

-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------EQTYYDLVKAF-AEIRQYIRELNLIIKVFREPFVSNSKLFSANDVENIFSRIVDIHELSVKLLGHIEDTVE-TDEGSPHPLVGSCFEDLAEELAFDPYESYARDILRPGFHDRFLSQLSKPGAALYLQSIGEGFKEAVQYVLPRLLLAPVYHCLHYFELLKQLEEKSEDQEDKECLKQAITALLNVQSG-EKICSKSLAKRRLSESA-------------AIKK-NEIQKNIDGWEGKDIGQCCNEFI-EGTLTRVGAKHERHIFLFDGL-ICCKSNHGQPRLPGASNAEYRLKEKFF-RKVQINDKDDTNEYKHAFEIILKDENSVIFSAKSAEEKNNW-AALISLQYRSTL---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
MQAQQLPYEFFSEENAPKWRGLLVPALKKVQGQVHPTLESNDDALQYVEELILQLLNMLCQAQPRSASDVEERVQKSFPHPIDKWAIADAQSAIEKRKRRNPLSLPVEKIHPLLKEVLGYKIDHQVSVYIVAVLEYISADILKLVGNYVRNIRHYEITKQDIKVAMCADKVLMDMFHQDVEDINILSLTDEEPSTSGEQTYYDLVKAFMAEIRQYIRELNLIIKVFREPFVSNSKLFSANDVENIFSRIVDIHELSVKLLGHIEDTVEMTDEGSPHPLVGSCFEDLAEELAFDPYESYARDILRPGFHDRFLSQLSKPGAALYLQSIGEGFKEAVQYVLPRLLLAPVYHCLHYFELLKQLEEKSEDQEDKECLKQAITALLNVQSGMEKICSKSLAKRRLSESACRFYSQQMKGKQLAIKKMNEIQKNIDGWEGKDIGQCCNEFIMEGTLTRVGAKHERHIFLFDGLMICCKSNHGQPRLPGASNAEYRLKEKFFMRKVQINDKDDTNEYKHAFEIILKDENSVIFSAKSAEEKNNWMAALISLQYRSTLERMLDVTMLQEEKEEQMRLPSADVYRFAEPDSEENIIFEENMQPKAGIPIIKAGTVIKLIERLTYHMYADPNFVRTFLTTYRSFCKPQELLSLIIERFEIPEPEPTEADRIAIENGDQPLSAELKRFRKEYIQPVQLRVLNVCRHWVEHHFYDFERDAYLLQRMEEFIGTVRGKAMKKWVESITKIIQRKKIARDNGPGHNITFQSSPPTVEWHISRPGHIETFDLLTLHPIEIARQLTLLESDLYRAVQPSELVGSVWTKEDKEINSPNLLKMIRHTTNLTLWFEKCIVETENLEERVAVVSRIIEILQVFQELNNFNGVLEVVSAMNSSPVYRLDHTFEQIPSRQKKILEEAHELSEDHYKKYLAKLRSINPPCVPFFGIYLTNILKTEEGNPEVLKRHGKELINFSKRRKVAEITGEIQQYQNQPYCLRVESDIKRFFENLNPMGNSMEKEFTDYLFNKSLEIEPRNPKPLPRFPKKYSYPLKSPGVRPSNPRPGTMRHPTPLQQEPRKISYSRIPESETESTASAPNSPRTPLTPPPASGASSTTDVCSVFDSDHSSPFHSSNDTVFIQVTLPHGPRSASVSSISLTKGTDEVPVPPPVPPRRRPESAPAESSPSKIMSKHLDSPPAIPPRQPTSKAYSPRYSISDRTSISDPPESPPLLPPREPVRTPDVFSSSPLHLQPPPLGKKSDHGNAFFPNSPSPFTPPPPQTPSPHGTRRHLPSPPLTQEVDLHSIAGPPVPPRQSTSQHIPKLPPKTYKREHTHPSMHRDGPPLLENAHSS
PDB_seq 100 # pdb to itself would obviously have 100% identity
Uniprot_seq 24 # pdb sequence has 24% identity to the uniprot sequence

Чтобы это работало с вашим входным файлом, вам нужно поместить мой a.get_seq() в цикл for с входными данными из вашего текстового файла.

ИЗМЕНИТЬ:

Замените функцию seq_id на эту:

def seq_id(self,aln_fasta):
    import string
    from Bio import AlignIO
    from Bio import SeqIO

    record_iterator = SeqIO.parse(aln_fasta, "fasta")
    first_record = record_iterator.next()
    print '%s has a length of %d' % (first_record.id, len(str(first_record.seq).replace('-','')))
    second_record = record_iterator.next()
    print '%s has a length of %d' % (second_record.id, len(str(second_record.seq).replace('-','')))

    lengths = [len(str(first_record.seq).replace('-','')), len(str(second_record.seq).replace('-',''))]
    if lengths.index(min(lengths)) == 0: # If both sequences have the same length the PDB sequence will be taken as the shortest
        print 'PDB sequence has the shortest length'
    else:
        print 'Uniport sequence has the shortes length'

    idenities = 0   
    for i,v in enumerate(first_record.seq):
        if v == '-':
            pass
            #print i,v, second_record.seq[i]
        if v == second_record.seq[i]:
            idenities +=1
            #print i,v, second_record.seq[i], idenities

    print 'Sequence Idenity = %.2f percent' % (100.0*(idenities/min(lengths)))

чтобы передать аргументы классу, используйте:

with open('input_file.txt', 'r') as infile:
    next(infile)
    next(infile) # Going by your input file
    for line in infile:
        line = line.split()
        a.get_seq(segs[0]+'.pdb',segs[1]+'.fasta')

Это может быть что-то вроде этого; повторяемый пример (например, с короткими файлами, размещенными в Интернете) поможет...

library(Biostrings)
pdb = readAAStringSet("pdb.fasta")
uniprot = readAAStringSet("uniprot.fasta")

для ввода всех последовательностей в два объекта. pairwiseAlignment принимает вектор в качестве первого аргумента (запроса), поэтому, если вы хотите выровнять все pdb со всеми uniprot, предварительно выделите матрицу результатов

pids = matrix(numeric(), length(uniprot), length(pdb),
              dimnames=list(names(uniprot), names(pdb)))

а потом делать расчеты

for (i in seq_along(uniprot)) {
    globalAlignment = pairwiseAlignment(pdb, uniprot[i])
    pids[i,] = pid(globalAlignment)
}