BioPython（一）

BioPython使用手册https://biopython-cn.readthedocs.io/zh_CN/latest/index.html#

1.1 Seq 对象

Biopython处理序列的机制为Seq 对象，Seq对象支持常规字符串的很多方法，Seq 对象具有一个重要的属性alphabet ，这一对象用于描述由单个字母构成的序列字符串的mean（意义），以及如何解释这一字符串。

>>>from Bio.Seq import Seq
>>>from Bio.Alphabet import IUPAC
>>>my_seq = Seq("GATCG", IUPAC.unambiguous_dna)
>>>my_seq
Seq('GATCG', IUPACUnambiguousDNA())

1.2 计算GC含量

使用Bio.SeqUtils.GC() 函数时会自动处理序列和可代表G或者C的歧意核苷酸 字母S混合的情况：

>>>from Bio.SeqUtils import GC
>>>my_seq = Seq('GATCGATGGGCCTATATAGGATCGAAAATCGC', IUPAC.unambiguous_dna)
>>>GC(my_seq)
46.875

1.3 核苷酸序列和（反向）互补序列：

.reverse_complement()方法:

>>>my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC", IUPAC.unambiguous_dna)
>>>my_seq
Seq('GATCGATGGGCCTATATAGGATCGAAAATCGC', IUPACUnambiguousDNA())
>>>my_seq.complement() #互补
Seq('CTAGCTACCCGGATATATCCTAGCTTTTAGCG', IUPACUnambiguousDNA())
>>>my_seq.reverse_complement() #反向互补
Seq('GCGATTTTCGATCCTATATAGGCCCATCGATC', IUPACUnambiguousDNA())

1.4 转录和翻译

.transcribe()方法，.translate()方法：

>>>from Bio.Seq import Seq
>>>from Bio.Alphabet import IUPAC
>>>coding_dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG", IUPAC.unambiguous_dna)
>>>coding_dna
Seq('ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG', IUPACUnambiguousDNA())
>>>template_dna = coding_dna.reverse_complement()
>>>template_dna
Seq('CTATCGGGCACCCTTTCAGCGGCCCATTACAATGGCCAT', IUPACUnambiguousDNA())
>>>messenger_rna = coding_dna.transcribe()#编码链直接转录
messenger_rna
Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG', IUPACUnambiguousRNA())
>>>template_dna.reverse_complement().transcribe() #模板链转录
Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG', IUPACUnambiguousRNA())
>>messenger_rna.back_transcribe() #逆转录
Seq('ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG', IUPACUnambiguousDNA())
>>>messenger_rna.translate() #翻译
Seq('MAIVMGR*KGAR*', HasStopCodon(IUPACProtein(), '*'))
>>>coding_dna.translate() #从编码链直接翻译
Seq('MAIVMGR*KGAR*', HasStopCodon(IUPACProtein(), '*'))
>>>coding_dna.translate(table="Vertebrate Mitochondrial")#通过指定密码子表进行特殊翻译如“线粒体密码子表”
Seq('MAIVMGRWKGAR*', HasStopCodon(IUPACProtein(), '*'))

1.5 SeqRecord (Sequence Record) 类

SeqRecord (Sequence Record) 类包含在Bio.SeqRecord 模块中，可以把identifiers 和featuresidentifiers等高级属性与序列关联起来。SeqRecord类中的format() 能将字符串转换成被 Bio.SeqIO 支持的格式。使用复杂，详见：https://biopython-cn.readthedocs.io/zh_CN/latest/cn/chr04.html

>>>from Bio.Seq import Seq
>>>from Bio.SeqRecord import SeqRecord
>>>from Bio.Alphabet import generic_protein
>>>record = SeqRecord(Seq("MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD" \
...                       +"GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK" \
...                       +"NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM" \
...                       +"SSAC", generic_protein),
...                    id="gi|14150838|gb|AAK54648.1|AF376133_1",
...                    description="chalcone synthase [Cucumis sativus]")
>>>record.format("fasta")
'>gi|14150838|gb|AAK54648.1|AF376133_1 chalcone synthase [Cucumis sativus]
MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD
GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK
NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM
SSAC'

1.6 序列输入和输出

1.6.1 解析/读取序列
Bio.SeqIO.parse(‘file’，‘format’) 用于读取序列文件生成 SeqRecord 对象，需要传入文件和格式两个参数，它返回一个 SeqRecord 对象迭代器（ iterator ），迭代器通常用在循环中。需要处理只包含一个序列条目的文件，可使用函数Bio.SeqIO.read()。

>>>from Bio import SeqIO
>>>for seq_record in SeqIO.parse("ls.fasta", "fasta"):
    print seq_record.id
    print repr(seq_record.seq)
    print len(seq_record)

1.6.2如果需要从文件中提取序列ID列表:

>>> from Bio import SeqIO
>>> identifiers = [seq_record.id for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank")]
>>> identifiers
['Z78533.1', 'Z78532.1', 'Z78531.1', 'Z78530.1', 'Z78529.1', 'Z78527.1', ..., 'Z78439.1']

1.6.3遍历序列
除了使用for循环，还可以使用迭代器的.next()方法遍历序列条目，使用 .next() 方法，当没有序列条目时，将抛出StopIteration 异常。序列文件包含多个序列条目，而只需要第一个条目时，此方法更简洁；Bio.SeqIO.read() 函数同样用于读取第一条目录。

>>>from Bio import SeqIO
>>>record_iterator = SeqIO.parse("ls_orchid.fasta", "fasta")
>>>first_record = record_iterator.next()
>>>second_record = record_iterator.next()

1.7序列文件作为字典

（1）Bio.SeqIO.to_dict()最灵活但内存占用最大，每个条目以SeqRecord对象形式存储在内存中，允许修改这些条目。genbank格式下默认会使用每条序列条目的ID（i.e. .id 属性）作为键，fasta格式下默认使用序列名作为键，使用可选参数key_function指定键，然后通过键读取操作相关序列。

>>>SM667 = SeqIO.to_dict(SeqIO.parse("SM667.faa", "fasta"))
>>>SM667.keys()
dict_keys(['fig|1548889.4.peg.1', 'fig|1548889.4.peg.2' ...'fig|1548889.4.peg.3']
>>>len(SM667.keys())  #获取序列数
4658

（2）对于更大的文件，应该考虑使用 Bio.SeqIO.index()，它并不将所有的信息存储在内存中，它仅仅记录每条序列条目在文件中的位置 ，当需要读取某条特定序列条目时，它才进行解析，它不支持比对文件格式，如PHYLIP或Clustal。
（3）Bio.SeqIO.index_db() 。由于它将序列信息以文件方式存储在硬盘上（使用SQLite3数据库）而不是内存中，因此它可以处理超大文件。同时，你可以同时对多个文件建立索引（前提是所有序列条目的ID是唯一的）。它需要三个参数：
1）索引文件名，建议使用以.idx 结尾的字符，改索引文件实质上是SQLite3数据库
2）要建立索引的文件列表（或者单个文件名）
3）文件格式

>>>gb_vrl = SeqIO.index_db("gbvrl.idx", [files_list], "genbank")

对于较大的文件，第一次运行会消耗较长时间，但是建立索引后再次使用便很快了。三种方式各有优缺点，视情况选择使用。
（4）要建立索引的文件非常大时，可以使用samtools的命令行工具 bgzip创建BGZF格式压缩文件，Bio.SeqIO.index()和 Bio.SeqIO.index_db() 函数均可以用于BGZF压缩文件。

$ bgzip -c ls.gbk > ls.gbk.bgz

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.index("ls.gbk.bgz", "genbank")

（5）获取序列条目原始数据
因为Biopython的GenBank和EMBL格式输出并不会保留每一点注释信息，所以当需要大文件中提取出一个序列子集或者需要完整地保留源文件时，使用 get_raw() 方法，它仅需要一个参数序列ID，然后返回一个字符串（提取自文件的未处理数据）。

根据序列号提取子集：

>>>from Bio import SeqIO
>>> uniprot = SeqIO.index("uniprot_sprot.dat", "swiss")
>>> handle = open("selected.dat", "w")
>>> for acc in ["P33487", "P19801", "P13689", "Q8JZQ5", "Q9TRC7"]:
...     handle.write(uniprot.get_raw(acc))
>>> handle.close()

1.8 写入序列文件

使用 Bio.SeqIO.parse()输入序列（读取文件），使用Bio.SeqIO.write()输出序列（写入文件）。该函数需要三个参数：某些 SeqRecord 对象，要写入的句柄或文件名，序列格式：

>>>from Bio import SeqIO
>>>SeqIO.write(my_records, "my_example.faa", "fasta")

1.9格式转换：

使用SeqRecord 对象列表作为Bio.SeqIO.write() 函数的输入，但是它也接受如来自于Bio.SeqIO.parse() 的 SeqRecord 迭代器 ，因此通过结合使用这两个函数实现文件转换。

>>>from Bio import SeqIO
>>>gbk = SeqIO.parse("file1.gbk", "genbank")
>>>fasta = SeqIO.write(gbk, "file2.fasta", "fasta")

或者使用.convert()方法

fasta = SeqIO.convert("file1.gbk", "genbank", "file2.fasta", "fasta")