Biopython 序列处理¶

pip install biopython                                # 安装
python -c "import Bio; print(Bio.__version__)"       # 验证版本

from Bio import SeqIO                                # 导入序列 IO
from Bio.Seq import Seq                              # 导入序列类

# 创建序列对象
dna = Seq("ATGCGATCGATCG")                           # 创建 DNA 序列
print(dna.complement())                              # 互补链
print(dna.reverse_complement())                      # 反向互补
print(dna.transcribe())                              # 转录为 RNA
print(dna.translate())                               # 翻译为蛋白质

# 读取 FASTA 文件
for record in SeqIO.parse("sequences.fasta", "fasta"):  # 逐条读取
    print(f"ID: {record.id}")                        # 序列 ID
    print(f"长度: {len(record.seq)}")                # 序列长度
    print(f"序列: {record.seq[:50]}...")              # 前 50 个碱基

# 读取单条序列
record = SeqIO.read("single.fasta", "fasta")         # 读取单条

# 读取 FASTQ 文件
for record in SeqIO.parse("reads.fastq", "fastq"):  # 读取 FASTQ
    qualities = record.letter_annotations["phred_quality"]  # 质量值
    avg_qual = sum(qualities) / len(qualities)       # 平均质量
    if avg_qual >= 30:                               # 过滤低质量
        print(f"{record.id}: avg_qual={avg_qual:.1f}")

# 写入文件
records = list(SeqIO.parse("input.fasta", "fasta"))  # 读取所有记录
SeqIO.write(records, "output.fasta", "fasta")        # 写入新文件

# 格式转换
SeqIO.convert("input.fastq", "fastq",               # 输入
              "output.fasta", "fasta")               # 输出（FASTQ→FASTA）

from Bio import Entrez                               # 导入 Entrez

Entrez.email = "your@email.com"                      # 必须设置邮箱

# 搜索 PubMed
handle = Entrez.esearch(
    db="pubmed",                                     # 数据库
    term="gut microbiome type 2 diabetes",           # 搜索词
    retmax=5                                         # 最多返回5条
)
results = Entrez.read(handle)
print(f"找到 {results['Count']} 篇文献")             # 结果数
print(results['IdList'])                              # 文献 ID 列表

# 获取 GenBank 序列
handle = Entrez.efetch(
    db="nucleotide",                                 # 核酸数据库
    id="NM_001301717",                               # 序列 ID
    rettype="gb",                                    # GenBank 格式
    retmode="text"
)
record = SeqIO.read(handle, "genbank")               # 解析
print(f"基因: {record.description}")                 # 描述
for feature in record.features:                      # 遍历特征
    if feature.type == "CDS":                        # 编码区
        print(f"CDS: {feature.location}")

from Bio.Blast import NCBIXML                        # 导入 BLAST 解析

# 解析 BLAST XML 结果
with open("blast_results.xml") as f:
    blast_records = NCBIXML.parse(f)                 # 解析 XML
    for record in blast_records:
        for alignment in record.alignments:          # 每个比对
            for hsp in alignment.hsps:               # 每个 HSP
                if hsp.expect < 1e-10:               # E-value 过滤
                    print(f"匹配: {alignment.title[:60]}")
                    print(f"E-value: {hsp.expect}")
                    print(f"Identity: {hsp.identities}/{hsp.align_length}")

from Bio import Phylo                                # 导入系统发育模块

# 读取树文件
tree = Phylo.read("tree.nwk", "newick")              # 读取 Newick 格式
Phylo.draw_ascii(tree)                               # ASCII 画树
tree.count_terminals()                               # 叶节点数

# 可视化
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(10, 8))
Phylo.draw(tree, axes=ax)                            # 绘制系统发育树
plt.savefig("phylo_tree.png", dpi=300)

from Bio import AlignIO                              # 导入比对 IO

# 读取比对文件
alignment = AlignIO.read("aligned.fasta", "fasta")   # 读取
print(f"序列数: {len(alignment)}")                   # 序列数量
print(f"比对长度: {alignment.get_alignment_length()}") # 比对长度

# 计算一致性
from Bio.Align import AlignInfo
summary = AlignInfo.SummaryInfo(alignment)
consensus = summary.dumb_consensus()                 # 一致序列
print(f"Consensus: {consensus}")

# === 序列操作 ===
seq.complement()                  # 互补链
seq.reverse_complement()          # 反向互补
seq.transcribe()                  # DNA → RNA
seq.translate()                   # DNA → 蛋白质
len(seq)                          # 序列长度

# === 文件读写 ===
SeqIO.parse(file, "fasta")       # 读取多条（迭代器）
SeqIO.read(file, "fasta")        # 读取单条
SeqIO.write(records, file, fmt)  # 写入
SeqIO.convert(in, fmt, out, fmt) # 格式转换

# === NCBI 查询 ===
Entrez.esearch(db, term)         # 搜索
Entrez.efetch(db, id)            # 获取记录

# === 支持格式 ===
# fasta, fastq, genbank, embl, clustal, nexus, phylip, stockholm