Biopython详解

安装：

pip install biopython
# 或
conda install -c conda-forge biopython

序列处理（Seq对象）：

from Bio.Seq import Seq
from Bio.SeqUtils import GC, molecular_weight
# 创建序列对象
dna_seq = Seq("ATGCGCTAGCTA")
# 序列属性
print(f"长度: {len(dna_seq)}")
print(f"GC含量: {GC(dna_seq):.1f}%")
print(f"分子量: {molecular_weight(dna_seq):.1f}")
# 序列操作
print(dna_seq.complement())        # 互补序列
print(dna_seq.reverse_complement()) # 反向互补
print(dna_seq.transcribe())        # DNA转RNA
print(dna_seq.translate())         # 翻译为蛋白质
# 转录和翻译
mrna = dna_seq.transcribe()
protein = mrna.translate()
print(f"Protein: {protein}")

解析FASTA/FASTQ文件：

from Bio import SeqIO
# 解析FASTA文件
for record in SeqIO.parse("sequences.fasta", "fasta"):
    print(f"ID: {record.id}")
    print(f"Description: {record.description}")
    print(f"Length: {len(record.seq)}")
    print(f"Sequence: {record.seq[:50]}...")
    print()
# 解析FASTQ文件（测序数据）
for record in SeqIO.parse("reads.fastq", "fastq"):
    print(f"ID: {record.id}")
    print(f"Sequence: {record.seq}")
    print(f"Quality: {record.letter_annotations['phred_quality'][:10]}")
# 将记录写入文件
records = []
for record in SeqIO.parse("input.fasta", "fasta"):
    if len(record.seq) > 100:  # 只保留长序列
        records.append(record)
SeqIO.write(records, "filtered.fasta", "fasta")
# 转换格式（FASTQ to FASTA）
SeqIO.convert("reads.fastq", "fastq", "reads.fasta", "fasta")

访问NCBI数据库：

from Bio import Entrez, SeqIO
# 设置邮箱（NCBI要求）
Entrez.email = "your.email@example.com"
# 搜索PubMed
handle = Entrez.esearch(db="pubmed", term="CRISPR[Title] AND 2023[PDAT]", retmax=10)
record = Entrez.read(handle)
print(f"Found {record['Count']} articles")
print(f"IDs: {record['IdList']}")
# 获取序列
handle = Entrez.efetch(db="nucleotide", id="NM_001301717", rettype="fasta", retmode="text")
record = SeqIO.read(handle, "fasta")
print(f"Sequence: {record.seq[:100]}...")

序列比对：

from Bio import pairwise2
from Bio.pairwise2 import format_alignment
# 全局比对（Needleman-Wunsch）
alignments = pairwise2.align.globalxx("ATCG", "ATG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 局部比对（Smith-Waterman）
alignments = pairwise2.align.localxx("ATCGGCTA", "CGG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 带参数的比对（匹配+1，错配-1，空位开启-2，空位延伸-1）
alignments = pairwise2.align.globalms("ATCG", "ATG", 1, -1, -2, -1)

BLAST解析：

from Bio.Blast import NCBIXML
# 解析BLAST XML结果
with open("blast_result.xml") as result_handle:
    blast_record = NCBIXML.read(result_handle)
    for alignment in blast_record.alignments:
        print(f"Hit: {alignment.title}")
        for hsp in alignment.hsps:
            print(f"  E-value: {hsp.expect}")
            print(f"  Identity: {hsp.identities}/{hsp.align_length}")
            print(f"  Query: {hsp.query[:50]}...")
            print(f"  Match: {hsp.match[:50]}...")
            print(f"  Sbjct: {hsp.sbjct[:50]}...")