pysam

Pysam - 基因组比对

用于高通量测序工作流。它支持高效访问比对到参考基因组的数十亿个 DNA 片段。

适用场景

• 处理新一代测序（NGS）数据。
• 分析基因组变异（SNP、indel）。
• 提取特定基因组区域的读段。
• 构建自定义生物信息学工作流。
• 对测序数据进行质量控制。

核心原则

索引访问

BAM 文件必须建立索引（.bai），才能高效地随机访问基因组区域。

坐标系统

基因组位置坐标从 0 开始（Python 风格），但在某些情况下，范围坐标从 1 开始。

读段属性

每个读段都包含序列、质量分数、比对位置和标志。

快速参考

标准导入

import pysam

基本模式

# 1. Open BAM file
samfile = pysam.AlignmentFile("aligned_reads.bam", "rb")

# 2. Iterate over reads in a specific genomic region
for read in samfile.fetch("chr1", 10000, 10100):
    print(f"Read: {read.query_name}, Quality: {read.mapping_quality}")
    print(f"Sequence: {read.query_sequence}")
    print(f"Position: {read.reference_start}")

# 3. Variant analysis (VCF)
vcf = pysam.VariantFile("mutations.vcf")
for rec in vcf.fetch("chr1", 10000, 10100):
    print(f"Pos: {rec.pos}, Ref: {rec.ref}, Alt: {rec.alts}")
    print(f"Genotype: {rec.samples['sample1']['GT']}")

# 4. Writing aligned reads
outfile = pysam.AlignmentFile("output.bam", "wb", template=samfile)
for read in samfile:
    if read.mapping_quality > 30:
        outfile.write(read)
outfile.close()

关键规则

✅ 应该做

• 始终使用已建立索引的文件 - 使用 pysam.index("file.bam") 创建索引，以便快速访问。
• 检查读段标志 - 使用 read.is_paired、read.is_unmapped 筛选读段。
• 处理未比对的读段 - 未比对读段的 reference_start = -1。
• 显式关闭文件 - 使用上下文管理器或 .close()，避免资源泄漏。

❌ 不要做

• 不要遍历整个 BAM 文件 - 使用带区域参数的 fetch() 提高效率。
• 不要忽略质量分数 - 低质量碱基可能导致误报变异。
• 不要混用坐标系统 - 统一使用从 0 开始或从 1 开始的索引方式。

高级模式

按基因统计读段

# Using a gene annotation file
genes = {}  # gene_name -> (chr, start, end)
for read in samfile.fetch():
    # Check if read overlaps any gene
    for gene, (chr, start, end) in genes.items():
        if read.reference_name == chr and start <= read.reference_start < end:
            genes[gene]['count'] += 1

变异筛选

# Filter variants by quality and depth
for rec in vcf.fetch():
    depth = rec.samples['sample1']['DP']
    qual = rec.qual
    if depth > 10 and qual > 20:
        # Process high-quality variant
        pass

Pysam 提供了基因组数据处理所需的底层访问能力，让研究人员能够直接处理生命本身的原始数据。