使用HTSeq计数转录组比对至细菌基因组CDS时无特征返回求助

• 参考序列ID完全不匹配：你用BWA比对的是单个CDS组成的参考序列库，SAM文件里reads的比对目标是每个CDS的ID（比如fig|203124.6.peg.1），但给HTSeq的GFF是基于全基因组染色体（accn|NC_008312）的注释，两者的参考序列ID没有任何交集，HTSeq自然找不到reads对应的特征区间，所有CDS计数为0。
• -i参数选得不对：你指定了-i product，但product是基因产物名称，不仅不唯一（比如示例里有两个CDS的product完全相同），而且和SAM里的参考ID完全无关，HTSeq无法通过这个字段关联reads和注释。
• 比对逻辑和注释不匹配：HTSeq的核心是统计reads在基因组区间的重叠，如果你直接比对到CDS序列，本质上每个CDS就是一个独立的“参考序列”，根本不需要用全基因组GFF来计数。

方案1：直接统计CDS参考的reads数（最省事） #

既然已经比对到CDS参考序列，没必要用HTSeq，直接用工具统计每个参考序列的reads数：

• 用samtools快速统计：

  samtools idxstats Test.sam > cds_counts.txt
  

  输出结果每一行的第二列就是对应CDS的总reads数，直接用就行。
• 用featureCounts（更灵活，支持处理配对端、过滤重复等）：

  featureCounts -a your_cds_ref.fasta -o cds_counts.txt Test.sam
  

  这里your_cds_ref.fasta就是你用来建BWA索引的CDS参考序列文件，featureCounts会自动对应每个CDS的reads数。

方案2：修改GFF适配当前SAM（非要用HTSeq的话） #

把全基因组GFF改成基于CDS序列的注释文件，让序列ID和SAM里的匹配：

1. 从你的CDS参考fasta文件中提取每个CDS的ID和序列长度。
2. 生成简化的GFF3文件，每一行对应一个CDS，格式示例：

   ##gff-version 3
   fig|203124.6.peg.1  PATRIC  CDS     1       867     .       +       0       ID=fig|203124.6.peg.1;product=Chromosomal replication initiator protein DnaA
   fig|203124.6.peg.2  PATRIC  CDS     1       201     .       -       0       ID=fig|203124.6.peg.2;product=Mobile element protein
   

3. 用修改后的GFF和正确参数跑HTSeq：

   htseq-count -a 0 -t CDS -i ID -m intersection-strict -s no Test.sam modified_cds.gff
   

   这里-i ID指定用CDS的唯一ID作为计数标识，确保和SAM里的参考ID完全对应。

方案3：重新比对到全基因组（适合后续基因组分析） #

如果之后需要做更多基因组层面的分析，建议重新比对到全基因组：

1. 用全基因组fasta建BWA索引：

   bwa index genome.fasta
   

2. 重新比对reads：

   bwa mem genome.fasta your_reads.fastq > genome_aligned.sam
   

3. 用原GFF和正确参数跑HTSeq：

   htseq-count -a 0 -t CDS -i ID -m intersection-strict -s no genome_aligned.sam Tricho.gff
   

   此时SAM里的参考序列ID（accn|NC_008312）和GFF里的序列ID完全匹配，-i ID用CDS的唯一ID计数，就能得到正确结果。

内容的提问来源于stack exchange，提问作者Clayton Tracey