Kodomo

Пользователь

Учебная страница курса биоинформатики,
год поступления 2014

BioPython

Python активно используется биоинформатиками. Для работы с биоинформатическими данными создан пакет BioPhython. В нем есть много полезных модулей для работы с последовательностями, геномами, филогенитическими деревьями. Сегодня разберем некоторые его возможности для работы с последовательностями. Ссылка на станицу BioPython:

http://biopython.org/wiki/Main_Page

Seq

Модуль для работы с последовательностями. В Biopython последовательности представлены в виде объектов Seq, которые содержат строку последовательности и алфавит, из которого составлена строка: 

   1 >>> from Bio.Seq import Seq
   2 >>> my_seq = Seq("AGTACACTGGT")
   3 >>> my_seq
   4 Seq('AGTACACTGGT', Alphabet())
   5 >>> my_seq.alphabet
   6 Alphabet()

В данном примере мы не задали алфавит при создании последовательности, поэтому, в переменной alphabet лежит алфавит по-умолчанию. Можно при создании последовательности указать конкретный алфавит, например, ДНК, РНК или белковая последовательность: 

   1 >>> from Bio.Seq import Seq
   2 >>> from Bio.Alphabet import generic_dna, generic_protein
   3 >>> my_seq = Seq("AGTACACTGGT")
   4 >>> my_seq
   5 Seq('AGTACACTGGT', Alphabet())
   6 >>> my_dna = Seq("AGTACACTGGT", generic_dna)
   7 >>> my_dna
   8 Seq('AGTACACTGGT', DNAAlphabet())
   9 >>> my_protein = Seq("AGTACACTGGT", generic_protein)
  10 >>> my_protein
  11 Seq('AGTACACTGGT', ProteinAlphabet())

Зачем это нужно и почему это важно? При работе с биологическими последовательностями использование алфавита является дополнительным контролем, это позволяет отслеживать ошибки. Например, не позволит объединять в одну последовательности ДНК и белка: 

   1 >>> my_protein + my_dna
   2 Traceback (most recent call last):
   3 ...
   4 TypeError: Incompatable alphabets ProteinAlphabet() and DNAAlphabet()

Кроме того, например, Biopython не позволит применить методы, используемые только для последовательностей ДНК (такие как трансляция) к белковым последовательностям.

Основные методы

Объект Seq имеет набор методов, которые работают также как методы для строк в Python. Например, метод find: 

   1 >>> from Bio.Seq import Seq
   2 >>> from Bio.Alphabet import generic_dna
   3 >>> my_dna = Seq("AGTACACTGGT", generic_dna)
   4 >>> my_dna
   5 Seq('AGTACACTGGT', DNAAlphabet())
   6 >>> my_dna.find("ACT")
   7 5
   8 >>> my_dna.find("TAG")
   9 -1

Есть также метод для посчета количества вхождений: 

   1 >>> my_dna.count("A")
   2 3
   3 >>> my_dna.count("ACT")
   4 1

Однако, будьте осторожны, также как и в методе для строк в Python, count считает число неперекрывающихся строк! 

   1 >>> "AAAA".count("AA")
   2 2
   3 >>> Seq("AAAA", generic_dna).count("AA")
   4 2

Методы для нуклеотидных последовательностей

Для нуклеотидных последовательностей, например, можно построить комплементарную и обратную комплементарную последовательность: 

   1 >>> from Bio.Seq import Seq
   2 >>> from Bio.Alphabet import generic_dna
   3 >>> my_dna = Seq("AGTACACTGGT", generic_dna)
   4 >>> my_dna
   5 Seq('AGTACACTGGT', DNAAlphabet())
   6 >>> my_dna.complement()
   7 Seq('TCATGTGACCA', DNAAlphabet())
   8 >>> my_dna.reverse_complement()
   9 Seq('ACCAGTGTACT', DNAAlphabet())

Транскрибция и обратная транскрибция

Можно преобразовать последовательность ДНК в РНК: 

   1 >>> my_dna
   2 Seq('AGTACACTGGT', DNAAlphabet())
   3 >>> my_dna.transcribe()
   4 Seq('AGUACACUGGU', RNAAlphabet())

А можно наоборот, преобразовать РНК в ДНК: 

   1 >>> my_rna = my_dna.transcribe()
   2 >>> my_rna
   3 Seq('AGUACACUGGU', RNAAlphabet())
   4 >>> my_rna.back_transcribe()
   5 Seq('AGTACACTGGT', DNAAlphabet())
   6 >>> my_rna.back_transcribe().reverse_complement()
   7 Seq('ACCAGTGTACT', DNAAlphabet())

Трансляция

Можно транслировать RNA: 

   1 >>> from Bio.Seq import Seq
   2 >>> from Bio.Alphabet import generic_rna
   3 >>> messenger_rna = Seq("AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG", generic_rna)
   4 >>> messenger_rna.translate()
   5 Seq('MAIVMGR*KGAR*', HasStopCodon(ExtendedIUPACProtein(), '*'))

Или ДНК - которая предполагается на кодирующей цепи: 

   1 >>> from Bio.Seq import Seq
   2 >>> from Bio.Alphabet import generic_dna
   3 >>> coding_dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG", generic_dna)
   4 >>> coding_dna.translate()
   5 Seq('MAIVMGR*KGAR*', HasStopCodon(ExtendedIUPACProtein(), '*'))

Можно установить опцию - транслировать до стоп кодона: 

   1 >>> coding_dna.translate(to_stop=True)
   2 Seq('MAIVMGR', ExtendedIUPACProtein())

Отметим, что если применить методы для нуклеотидных последовательностей к белковым, то это вызовет исключение: 

   1 >>> my_protein.complement()
   2 Traceback (most recent call last):
   3 ...
   4 ValueError: Proteins do not have complements!

SeqIO

Данный модуль в Bio.SeqIO позволяет осуществлять ввод и вывод последовательностей в различных форматах (включая множественные выравнивания). Последовательности преобразуются в объекты SeqRecord.

Форматы файлов:

опция

формат файла

clustal

выравнивания в формате Clustal X и Clustal W

fasta

Fasta

stockholm

PFAM

swiss

Swiss-Prot

И многие другие (см. http://biopython.org/wiki/SeqIO )

Ввод последовательности

Для загрузки последовательности используется функция Bio.SeqIO.parse() которая берет на вход файл и имя формата, и возвращает итератор SeqRecord (штука, позволяющая перебирать последовательности одну за другой). Например: 

   1 from Bio import SeqIO
   2 handle = open("example.fasta")
   3 for record in SeqIO.parse(handle, "fasta") :
   4     print record.id
   5 handle.close()

Итератор удобен, когда последовательности нужно перебрать одну за другой в порядке, в котором они идут в файле.Для некоторых задач бывает полезно получить все записи в виде списка: 

   1 from Bio import SeqIO
   2 handle = open("example.fasta")
   3 records = list(SeqIO.parse(handle, "fasta"))
   4 handle.close()
   5 print records[0].id  #first record
   6 print records[-1].id #last record

Другая распространенная задача - извлечь последовательность по id. Для работы с маленькими файлами используется функция Bio.SeqIO.to_dict(), которая преобразует список в словарь: 

   1 from Bio import SeqIO
   2 handle = open("example.fasta")
   3 record_dict = SeqIO.to_dict(SeqIO.parse(handle, "fasta"))
   4 handle.close()
   5 print record_dict["gi:12345678"] #use any record ID

Для файлов большого размера может не хватить памяти, в этом случае используется Bio.SeqIO.index(). 

   1 from Bio import SeqIO
   2 record_dict = SeqIO.index("example.fasta", "fasta")
   3 print record_dict["gi:12345678"] #use any record ID

Для записи последовательностей в файл используется функция Bio.SeqIO.write(), которая берет на вход SeqRecord итератор (или список), выходной файл и формат: 

   1 from Bio import SeqIO
   2 sequences = ... # add code here
   3 output_handle = open("example.fasta", "w")
   4 SeqIO.write(sequences, output_handle, "fasta")
   5 output_handle.close()

Преобразование форматов

Предположим, у нас есть файл в формате genbank, мы хотим преобразовать его в fasta формат: 

   1 from Bio import SeqIO
   2  
   3 input_handle = open("cor6_6.gb")
   4 output_handle = open("cor6_6.fasta", "w")
   5  
   6 sequences = SeqIO.parse(input_handle, "genbank")
   7 count = SeqIO.write(sequences, output_handle, "fasta")
   8  
   9 output_handle.close()
  10 input_handle.close()
  11  
  12 print "Converted {0} records".format(count)

Или еще более коротко, с использованием функции Bio.SeqIO.convert(): 

   1 from Bio import SeqIO
   2 count = SeqIO.convert("cor6_6.gb", "genbank", "cor6_6.fasta", "fasta")
   3 print "Converted {0} records".format(count)

Бывает необходимо произвести какие-нибудь манипуляции с последовательностями, отфильтровать по какому-то признаку (например, по длине): 

   1 from Bio import SeqIO
   2  
   3 short_sequences = [] # Setup an empty list
   4 for record in SeqIO.parse(open("cor6_6.gb"), "genbank") :
   5     if len(record.seq) < 300 :
   6         # Add this record to our list
   7         short_sequences.append(record)
   8  
   9 print "Found {0} short sequences".format(len(short_sequences))

Или еще короче (применяем сокращенную запись из предыдущего занятия): 

   1 from Bio import SeqIO
   2  
   3 input_seq_iterator = SeqIO.parse(open("cor6_6.gb"), "genbank")
   4  
   5 #Build a list of short sequences:
   6 short_sequences = [record for record in input_seq_iterator \
   7                    if len(record.seq) < 300]
   8  
   9 print "Found {0} short sequences".format(len(short_sequences))

Можно вырезать часть последовательности и сделать из нее новую, чтобы затем записать в файл: 

   1 >>>from Bio.SeqRecord import SeqRecord
   2 >>> handle=open("copy.fasta")
   3 >>> records = list(SeqIO.parse(handle, "fasta"))
   4 >>> seq=records[0].seq[0:20]
   5 >>> records[0]=SeqRecord(seq, id=records[0].id, name=records[0].name, description=records[0].description)
   6 >>> print records[0]
   7 ID: 2a06_C
   8 Name: 2a06_C
   9 Description: 2a06_C
  10 Number of features: 0
  11 Seq('MTNIRKSHPLMKIVNNAFID', SingleLetterAlphabet())