Biopython入门

Biopython入门

1. 什么是Biopython？

Biopython工程是一个使用Python来开发计算分子生物学工具的国际团体。
Biopython官网(http://www.biopython.org)为使用和研究生物信息学的开发者提供了一个在线资源库，包括模块、脚本以及一些基于Python的软件的网站链接。Biopython致力于通过创造高质量的和可重复利用的模块及类，从而使得Python在生物信息学中的应用变得更加容易。Biopython的特点包括解析各种生物信息学格式的文件(BLAST， Clustalw， FASTA， Genbank…)，访问在线的服务器(NCBI，Expasy…)，常见和不那么常见程序的接口(Clustalw， DSSP，MSMS…)，标准的序列类，各种收集的模块，KD树数据结构等等，还有一些文档。

2. Seq序列

Biopython序列是指一系列字母，用于表示生物体的蛋白质、DNA或RNA，它由Seq类表示。
Seq类在Bio.Seq模块中定义。

引入Seq

from Bio.Seq import Seq

创建一个序列

seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC")

3. 序列基本操作

序列类似于python字符串。可按python字符串操作，例如切片，计数，串联，查找，拆分和剥离。

获取第一个值

seq[0]

切取序列

seq[2:9]

和Python字符串一样，可以通过设置起始位置、终止位置和 步幅 （间隔数，默认为1） 进行切片，
例如，分别获取下面DNA序列密码子第一、第二、第三位的碱基。

seq[0::3]
seq[1::3]
seq[2::3]

长度

len(seq)

计数

seq.count('A')
seq.count('TA')

计算GC含量

from Bio.SeqUtils import GC
GC(seq)

添加序列

seq1 = Seq("AGCT")
seq2 = Seq("TCGA")
seq3 = seq1+seq2

添加多个序列，将其存储在python列表中，然后使用for循环进行检索，最后将其添加在一起

seqs = [Seq("AGCT"),Seq("TCGA"),Seq("AAA")]
final_seq = Seq("")
for seq in seqs:
    final_seq = final_seq + seq

更改序列的大小写

seq.lower()
seq.upper()

序列匹配

"GTAC" in seq
"GTAC" in seq.upper()

查找给定序列内的单个字母或字母序列

seq1.find('C')
seq1.find('CT')

拆分序列

seq.split('C')

剥离操作

seq = Seq(" AGCT ")
print(seq.strip())

将序列对象转换成字符串

str(seq)

反向序列

seq[::-1]

互补序列

seq.complement()

反向互补序列

seq.reverse_complement()

转录

转录DNA成mRNA
将T → U转换
transcribe()

seq.transcribe()

如果想从模板链去做一个真正的生物学上的转录，需要两步，结果相同

seq.reverse_complement().transcribe()

逆转录
back_transcribe()

seq1 = seq.transcribe()
print(seq1.back_transcribe())

转录翻译DNA成蛋白序列
translate()
可以直接从DNA链进行翻译

seq.translate()

或者将DNA转录为mRNA，再翻译成相应的蛋白质序列

mrna_seq = seq.transcribe()
mrna_seq.translate()

终止密码子用*表示，要让translate()方法在第一个终止密码子处终止，可以在translate()中分配参数：to_stop = True，当使用 to_stop 参数时，终止密码子本身是不翻译的，终止的符号也不显示在蛋白质序列中

seq.translate(to_stop=True)

biopython的翻译表是基于NCBI，可以根据实际设定table，默认情况下，翻译使用的是标准遗传密码（NCBI上 table id 为1)
若需要翻译一个线粒体序列，我们就需要告诉翻译函数使用相关的遗传密码

seq.translate(table="Vertebrate Mitochondrial", to_stop=True)  #设置为线粒体模式，遇到终止密码子则停止

或者，可以利用NCBI上表格的标号来指定所使用的遗传密码

seq.translate(table=2, to_stop=True)

在translate()中分配参数：cds=True 说明输入的是完整的cds序列，不含终止密码子，如果出现终止密码子则会抛出异常

print(seq.translate(cds=True))

翻译表
在这里先关注：标准的翻译表和脊椎动物线粒体DNA的翻译表
引入CodonTable，打印翻译表

from Bio.Data import CodonTable
standard_table = CodonTable.unambiguous_dna_by_name["Standard"]
mito_table = CodonTable.unambiguous_dna_by_name["Vertebrate Mitochondrial"]

或者以标号1和2来标识

standard_table = CodonTable.unambiguous_dna_by_id[1]
mito_table = CodonTable.unambiguous_dna_by_id[2]

打印standard_table

查找终止密码子

standard_table.stop_codons

查找起始密码子

standard_table.start_codons

查找对应的氨基酸

standard_table.forward_table["ACG"]

4. MutableSeq对象

Seq 对象是 “只读的” ，不可变的。
可以使用 tomutable对象将它转换成可变的序列

mutable_seq = seq.tomutable()

或者直接建立一个 MutableSeq 对象

from Bio.Seq import MutableSeq
mutable_seq = MutableSeq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA")

对序列进行变换

mutable_seq[5] = "C"
mutable_seq.remove("T")
mutable_seq.reverse()

与Seq对象不同的是， MutableSeq 对象的各种方法都是实时呈现的，比如reverse_complement() 和 reverse() 方法。差别就是 MutableSeq 对象不可以作为字典的键值 ，但是Seq对象就可以。

5. SeqRecord对象

SeqRecord (Sequence Record) 类包含在 Bio.SeqRecord 模块中。该类是 Bio.SeqIO 序列输入/输出交互的基本数据类型。

SeqRecord 类包括下列属性:

.seq：序列自身（即 Seq 对象）。
.id： 序列主ID（-字符串类型）。通常类同于accession number。
.name：序列名/id （-字符串类型）。 可以是accession number, 也可是clone名（类似GenBank record中的LOCUS id）。
.description：序列描述（-字符串类型）。
.letter_annotations：对照序列的每个字母逐字注释（per-letter-annotations），以信息名为键（keys），信息内容为值（value）所构成的字典。值与序列等长，用Python列表、元组或字符串表示。.letter_annotations可用于质量分数或二级结构信息 (如 Stockholm/PFAM 比对文件)等数据的存储。
.annotations：用于储存附加信息的字典。信息名为键（keys），信息内容为值（value）。用于保存序列的零散信息（如unstructured information）。
.features：SeqFeature 对象列表，储存序列的结构化信息（structured information），如：基因位置, 蛋白结构域。
.dbxrefs：储存数据库交叉引用信息（cross-references）的字符串列表。

使用 SeqRecord 对象非常简单，因为所有的信息都存储在该类的属性中；通常不必手动新建，用 Bio.SeqIO 从序列文件读取即可。

6. 序列输入和输出

Biopython提供了Bio.SeqIO模块来分别从文件读取序列和向文件写入序列。它支持生物信息学中几乎所有可用的文件格式，通过SeqRecord对象向用户显示已解析的序列数据。

通过ncbi下载LOC9271075基因序列

引入SeqIO

from Bio import SeqIO

读取单条序列

seq_record = SeqIO.read("filename.fna","fasta")
print(seq_record.seq)

**读取多条SeqRecord **

seq_records = SeqIO.parse("filename.fna","fasta")

遍历多条序列

for seq_record in seq_records:
    print(repr(seq_record.seq))

除了使用for循环，还可以使用迭代器的 .next() 方法遍历序列条目

seq_records = SeqIO.parse("filename.fasta", "fasta")

first_record = seq_records.next()
print(first_record.id)
print(first_record.description)

second_record = seq_records.next()
print(second_record.id)
print(second_record.description)

写入序列文件
Bio.SeqIO.write()

创建两个新的 SeqRecord 对象，得到一个 SeqRecord 对象列表，将它写入一个FASTA格式文件

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio import SeqIO

rec1 = SeqRecord(Seq("MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD" \
                    +"GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK" \
                    +"NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM" \
                    +"SSAC"),
                 id="gi|14150838|gb|AAK54648.1|AF376133_1",
                 description="chalcone synthase [Cucumis sativus]")
rec2 = SeqRecord(Seq("YPDYYFRITNREHKAELKEKFQRMCDKSMIKKRYMYLTEEILKENPSMCEYMAPSLDARQ" \
                    +"DMVVVEIPKLGKEAAVKAIKEWGQ"),
                 id="gi|13919613|gb|AAK33142.1|",
                 description="chalcone synthase [Fragaria vesca subsp. bracteata]")
my_records = [rec1, rec2]
SeqIO.write(my_records, "my_example.fna", "fasta")