HTSeq项目教程：深入解析htseq-count的工作原理与实现细节

引言

在RNA-Seq数据分析中，基因表达定量是一个基础而关键的步骤。HTSeq项目提供的htseq-count工具因其准确性和灵活性而广受欢迎。本文将深入解析这个"黑箱"工具的内部工作机制，帮助用户更好地理解其原理并掌握定制化使用的方法。

准备工作

要运行htseq-count，我们需要两个核心输入文件：

1. 基因组注释文件：通常为GTF格式，包含基因、外显子等特征的位置信息
2. 比对结果文件：SAM/BAM/CRAM格式，包含测序reads与基因组的比对信息

在本教程示例中，我们使用酵母(Saccharomyces cerevisiae)的数据：

• 注释文件：Saccharomyces_cerevisiae.SGD1.01.56.gtf.gz
• 测序数据：yeast_RNASeq_excerpt.sam

第一步：加载基因组注释GTF文件

1.1 读取GTF文件

import HTSeq
gtffile = HTSeq.GFF_Reader("Saccharomyces_cerevisiae.SGD1.01.56.gtf.gz")

1.2 构建基因组特征数据结构

feature_scan = HTSeq.make_feature_genomicarrayofsets(
    gtffile,
    id_attribute='gene_id',
    feature_type='exon',
)

这里创建了一个GenomicArrayOfSets数据结构，它高效地存储了基因组上各个区间的特征信息。关键参数说明：

• id_attribute='gene_id'：使用GTF文件中的gene_id属性作为特征标识
• feature_type='exon'：只处理GTF中的外显子特征

1.3 数据结构解析

feature_scan将基因组划分为多个区间(GenomicInterval)，每个区间关联一个或多个外显子。例如：

• 如果区间1566-1589仅被一个基因的外显子覆盖，则该区间关联该gene_id
• 如果区间1580-1589被两个基因的外显子重叠覆盖，则该区间关联这两个gene_id

这种设计使得后续的reads分类既高效又准确。

第二步：处理比对文件中的reads

2.1 读取比对文件

bamfile = HTSeq.BAM_Reader("yeast_RNASeq_excerpt.sam")

2.2 初始化计数数据结构

attributes = feature_scan['attributes']
feature_attr = sorted(attributes.keys())
counts = {key: 0 for key in feature_attr}

# 特殊分类计数
counts['notaligned'] = 0    # 未比对reads
counts['no_feature'] = 0    # 比对但不在任何特征上
counts['ambiguous'] = 0     # 比对到多个特征

2.3 遍历比对reads

for read in bamfile:
    if not read.aligned:
        counts['notaligned'] += 1
        continue

第三步：reads与基因特征的比对

3.1 获取reads的比对区间

aligned_codes = ('M', '=', 'X')
iv_read = (co.ref_iv for co in read.cigar if co.type in aligned_codes)

这里通过CIGAR字符串解析reads的实际比对区间，跳过插入和删除等比对特征。

3.2 查找重叠基因

gene_ids_read = None
for iv in iv_read:
    for _, gene_ids in feature_scan['features'][iv].steps():
        if gene_ids_read is None:
            gene_ids_read = gene_ids.copy()
        else:
            gene_ids_read.intersection(gene_ids)

这段代码实现了"intersection-strict"模式，即只有当read完全位于某个基因内时才计数。其他模式如"union"会有不同的处理逻辑。

第四步：reads分类计数

4.1 无特征重叠

if gene_ids_read is None or len(gene_ids_read) == 0:
    counts['no_feature'] += 1
    continue

4.2 多特征重叠(ambiguous)

if len(gene_ids_read) > 1:
    counts['ambiguous'] += 1
    continue

4.3 单特征重叠(唯一比对)

gene_id = list(gene_ids_read)[0]
counts[gene_id] += 1

第五步：资源清理

bamfile.close()
gtffile.close()

高级主题

外显子水平计数

通过修改make_feature_genomicarrayofsets的参数，可以实现外显子水平的计数：

feature_scan = HTSeq.make_feature_genomicarrayofsets(
    gtffile,
    id_attribute=['gene_id', 'exon_number'],
    feature_type='exon',
)

双端测序数据的特殊处理

双端测序数据的分析更为复杂，需要注意：

1. 比对文件可能是按名称排序或位置排序
2. 需要缓冲机制等待配对的reads
3. 建议使用未排序或按名称排序的BAM文件以提高效率

总结

本文详细剖析了htseq-count的内部工作机制，包括：

1. 基因组注释的加载与处理
2. 比对reads的分类策略
3. 不同计数模式的区别
4. 特殊情况的处理方法

理解这些底层原理不仅有助于正确使用工具，也为用户定制自己的分析流程奠定了基础。对于有特殊需求的用户，可以参考这些原理开发自己的计数脚本。