Kodomo

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

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

BioPython

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

http://biopython.org/wiki/Main_Page

Seq

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

>>> from Bio.Seq import Seq
>>> my_seq = Seq("AGTACACTGGT")
>>> my_seq
Seq('AGTACACTGGT', Alphabet())
>>> my_seq.alphabet
Alphabet()

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

>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import generic_dna, generic_protein
>>> my_seq = Seq("AGTACACTGGT")
>>> my_seq
Seq('AGTACACTGGT', Alphabet())
>>> my_dna = Seq("AGTACACTGGT", generic_dna)
>>> my_dna
Seq('AGTACACTGGT', DNAAlphabet())
>>> my_protein = Seq("AGTACACTGGT", generic_protein)
>>> my_protein
Seq('AGTACACTGGT', ProteinAlphabet())

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

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

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

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

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

>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import generic_dna
>>> my_dna = Seq("AGTACACTGGT", generic_dna)
>>> my_dna
Seq('AGTACACTGGT', DNAAlphabet())
>>> my_dna.find("ACT")
5
>>> my_dna.find("TAG")
-1

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

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

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

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

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

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

>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import generic_dna
>>> my_dna = Seq("AGTACACTGGT", generic_dna)
>>> my_dna
Seq('AGTACACTGGT', DNAAlphabet())
>>> my_dna.complement()
Seq('TCATGTGACCA', DNAAlphabet())
>>> my_dna.reverse_complement()
Seq('ACCAGTGTACT', DNAAlphabet())

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

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

>>> my_dna
Seq('AGTACACTGGT', DNAAlphabet())
>>> my_dna.transcribe()
Seq('AGUACACUGGU', RNAAlphabet())

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

>>> my_rna = my_dna.transcribe()
>>> my_rna
Seq('AGUACACUGGU', RNAAlphabet())
>>> my_rna.back_transcribe()
Seq('AGTACACTGGT', DNAAlphabet())
>>> my_rna.back_transcribe().reverse_complement()
Seq('ACCAGTGTACT', DNAAlphabet())

Трансляция

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

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

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

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

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

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

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

>>> my_protein.complement()
Traceback (most recent call last):
...
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 (штука, позволяющая перебирать последовательности одну за другой). Например:

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

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

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

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

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

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

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

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

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

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

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

from Bio import SeqIO
 
input_handle = open("cor6_6.gb")
output_handle = open("cor6_6.fasta", "w")
 
sequences = SeqIO.parse(input_handle, "genbank")
count = SeqIO.write(sequences, output_handle, "fasta")
 
output_handle.close()
input_handle.close()
 
print "Converted {0} records".format(count)

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

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

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

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

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

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

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

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

AlignIO

Модуль Bio.SeqIO работает с одной или несколькими записями как с объектами SeqRecord. Аналогично устроен модуль Bio.AlignIO, который предназначен для работы с выравниваниями последовательностей как с объектами Alignment. Bio.AlignIO использует такие же функции, как Bio.SeqIO, и такие же имена для поддерживаемых форматов.

Ввод

Есть две функции для ввода последовательностей. Это Bio.AlignIO.read() для тех случаев, когда файл содержит только одно выравнивание, и более общая Bio.AlignIO.parse(), когда файл содержит несколько отдельных выравниваний. Обе этих функции принимают на вход файл и формат. Например, у нас есть выравнивание в формате PFAM или Stockholm. Файл содержит одно выравнивание, можем использовать метод Bio.AlignIO.read():

from Bio import AlignIO
alignment = AlignIO.read(open("PF09395_seed.sth"), "stockholm")
print "Alignment length {0}".format(alignment.get_alignment_length())
for record in alignment :
    print record.seq, record.id

Также как и в Bio.SeqIO, для вывода существует единственная функция Bio.AlignIO.write().

Метод format() возвращает выравнивание в виде строки в указанном формате:

AlignIO.read(open("PF09395_seed.sth"), "stockholm")
print alignment.format("fasta")

Точно также, как в Bio.SeqIO, можно преобразовать один формат в другой:

from Bio import AlignIO
 
input_handle = open("example.phy")
output_handle = open("example.sth", "w")
 
alignment = AlignIO.read(input_handle, "phylip")
AlignIO.write(alignment, output_handle, "stockholm")
 
output_handle.close()
input_handle.close()

Срезы

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

alignment[0].seq[2]

или

alignment[0,2]

а можно выбрать диапозон последовательностей и диапозон колонок (кусок выравнивания):

alignment[0:10,2:10]

столбец:

alignment[:,2]

2012/4/Python/11/Record (последним исправлял пользователь stavrovskaya 2014-05-12 17:53:00)