做RNA-seq分析做到第十步终于到了要给每个基因算总账的时候。前面几篇我们聊了质控、比对、转录本组装每一步都在为最后这张表达量清单做准备。这篇笔记要说的就是定量环节里最经典也最容易踩坑的一步用HTSeq-count把比对好的reads数到每个基因头上。很多初学者会以为计数就是拿wc -l数一下比对文件里有多少行对齐记录再把基因ID对上就行。实际完全不是这么回事。一条read落在基因区域内的不同位置和多个注释feature重叠双端读段怎么算一个fragment链特异性文库要不要反转方向这些细节如果处理不对你拿到的count矩阵根本不能用来做差异表达分析。HTSeq-count这个工具虽然命令短、参数少但从输入文件规范到输出结果的解读每一步都藏着决定成败的坑。这篇文章基于我自己跑流程的完整记录从HTSeq-count在RNA-seq流程里的定位讲起把输入文件的准备工作、核心参数的语义、实际运行过程中遇到的报错和排查思路以及从单样本计数到批量生成表达矩阵的衔接方法一次性说清楚。适合正在自己搭建RNA-seq分析流程、或者第一次用HTSeq-count还不太确定参数该怎么配的朋友。1. 计数这一步为什么不能用数read数糊弄过去1.1 HTSeq-count在整个定量链条里的位置RNA-seq的标准分析流水线从原始数据到最终的基因表达量本质上是在做三件事把RNA片段的序列读出来把每个片段放回基因组上的位置然后统计每个基因被多少个片段覆盖。比对这一步解决的是放到哪里而HTSeq-count解决的是数多少个。这里的数并没有看起来那么简单。一个基因不是基因组上一条连续的线而是由多个外显子拼接而成转录组注释文件里记录的是一个个外显子区间。一条read可能落在外显子中间可能跨越外显子和内含子的边界还可能同时和两个基因的区域都有重叠——比如一对双端read的其中一段落在基因A末端另一段落在基因B开头。这些重叠情况如果处理不好计数结果就会失真。HTSeq-count之所以在众多定量工具里仍然被当作标准流程的首选是因为它采用了一种确定性很高、可完全解释的计数规则把每条read或read pair与注释基因的外显子区间做交集判断最常见的union模式下只要read和该基因任一外显子有重叠就计数进去但如果一条read同时和多个基因的重叠就把它归为ambiguous不做硬性分配。这和其他工具的思路有明显区别。比如featureCounts也做类似的外显子交集统计但它在处理多重叠时会用更复杂的优先级规则去强行分配而HTSeq-count的哲学是不确定的就不硬算。再比如Salmon、kallisto这类基于伪比对的定量工具根本不先生成BAM而是在转录组层面直接做丰度估计。HTSeq-count的优势在于它直接建立在比对结果之上出了任何问题都可以回溯检查每一条read的比对状态这在搭建标准化流程时特别重要——可审计这三个字在论文复现和数据质控上的价值远高于一点计算速度的节省。1.2 alignment不是fragment双端数据的计数单位容易被搞错这里要先纠正一个常见的理解偏差。很多人第一次看HTSeq-count的说明文档会默认一条比对记录就是一个计数单位然后就能理解为什么有时候输出文件的数字比samtools flagstat里比对上的read数少了将近一半。实际上对于双端测序数据HTSeq-count不会把两个mate分别计数而是把能够正确配对到同一个片段fragment的两个read合并成一条计数单元。这是从转录组定量的生物学意义出发的一个RNA片段经过PCR扩增和双端测序产生了两个read但它们本质上代表同一个原始分子。如果你把两个mate各算一次等于把每个片段的信号翻倍了不同基因的片段长度差异还会放大这种偏差。HTSeq-count内部做的事情是把pair的两条read取并集区域然后拿这个合并后的区间去和注释基因做overlap判断。这就带来一个直接影响输入HTSeq-count的BAM文件必须按坐标排序而不是按read名排序。因为要合并一个pair的两条read工具需要能在文件里快速找到位于同一坐标区域的mate。如果你的上游比对工具输出的是name-sorted BAM直接用HTSeq-count跑大概率会得到错误的结果甚至直接报排序错误。我实际跑流程时上游用到了某款比对工具时它的默认输出顺序经常是按处理顺序来的并不是坐标排序。所以我每次比对完成后都会固定执行一步samtools sort - 8 -m 4G -o sample.sorted.bam sample.bam samtools index sample.sorted.bam这一小步看起来不起眼但它是HTSeq-count能正确解析双端数据的前提。1.3 数清楚需要的是基因模型不是单纯的重叠计数HTSeq-count真正的工作机制是把比对结果和一套完整的基因模型做碰撞检测。基因模型来自GTF注释文件里面每一行记录一个feature可能是exon、CDS、transcript、gene等类型。默认情况下HTSeq-count只看exon类型的feature因为表达量的判定区域就是外显子。它会把这些exon按第9列里的gene_id属性分组凡是被同一gene_id标记的所有exon区间在计数逻辑上合并为这个基因的可计数区域。这就是为什么如果你拿到的GTF文件里gene_id这一属性写得不规范或者某些行的第9列格式有问题工具会在解析阶段就报错——因为它根本无法可靠地把外显子序列归堆到基因上。在后面的章节里我会把GTF文件第9列的格式坑专门展开。这里只需要记住一个核心结论HTSeq-count的数数规则是read片段与基因外显子区间的overlap判定而不是read与基因组区间位置的覆盖统计。理解了这一点后面看输出文件里那一堆__ambiguous、__no_feature的行就不会一头雾水了。2. 开跑之前先检查这两样东西能避开一半的坑2.1 比对文件的三个硬性要求HTSeq-count对输入的比对文件有三条硬性要求格式必须是BAM或SAM、双端数据必须按坐标排序、每条记录的质量字段必须合法且非空。第三个要求是很多人完全没意识到的也是最容易在深夜跑脚本时被恶心到的点。有些上游工具产出的BAM文件里部分read的比对质量字段第5列的MAPQ可能是255也可能是0。HTSeq-count默认会丢弃MAPQ值低于10的read这个阈值可以通过-a参数调整。如果你发现输出文件里__too_low_aQual这个特殊行的数值异常大先别急着怀疑工具用下面的命令看看你的BAM里整体质量分布samtools view sample.sorted.bam | awk {print $5} | sort -n | uniq -c | head -30如果低质量比对MAPQ10占比很高而你觉得这在生物学上是合理的比方说多比对reads处理策略那么-a 0可以临时把阈值关掉。但从流程规范化的角度我更建议保留默认阈值然后在结论里如实报告这个过滤比例而不是把它悄悄关掉否则下游差异分析的假阳性会明显上升。关于排序这里还有一个细节值得说。HTSeq-count的-r参数接受两个值pos和name。如果你的BAM是name-sorted的理论上可以用-r name让它先按名字把pair归到一起再处理坐标但实际体验下来绝大多数规范的流程都是先用samtools sort转成坐标排序然后固定用-r pos去跑。原因有两个一是坐标排序的BAM在后续IGV可视化、samtools depth等环节都是通用标准二是HTSeq-count处理name-sorted文件时内存占用明显偏高大样本很容易把服务器跑挂。2.2 GTF注释的版本一致性和格式细节GTF文件是计数环节的另一半。这里最重要的一个原则是GTF的基因组版本必须和比对时用的参考基因组版本一致。比如前面比对用的参考是某个版本的UCSC序列文件里面染色体命名可能是chr1而注释文件来自一个只写了1的版本那么HTSeq-count会报染色体名对不上的错误因为所有read的比对位置都在chr1上而GTF里根本没有这个名字的染色体。版本一致性问题很难靠工具自动纠正。我之前踩过一次上游比对用的参考序列是从某镜像站下载的染色体名是chr1这种格式而注释文件是从另一处下载的GTF里面用的是1。跑HTSeq-count时一秒钟就报错了报错信息明确说某个染色体不在注释文件里。解决方法是统一令注解文件中的染色体命名风格。这里建议在比对阶段就和注释文件保持同源最稳妥的做法是直接用注释文件所属项目提供的参考基因组序列或者用sed批量替换染色体名但如果决定替换一定要检查基因注释里偶尔出现的chrM、chrMT等线粒体染色体写法别全省了。GTF文件本身的结构也有强制性要求。第9列是属性列里面的键值对必须是key value;这样的格式每个键值对以分号结束值要加引号。由于HTSeq-count的-i参数默认取gene_id作为分组键所以你的GTF文件里每一行记录都必须包含gene_id这个键否则工具会直接提示属性不存在。我在实际项目中遇到过某个来源的GTF文件部分转录本的属性键用的是geneID而不是gene_id差一个大小写字母整个文件就无法解析。验证GTF文件是否干净可以先试试这个快速检查awk -F \t $3exon annotation.gtf | head -5 awk -F \t {print $9} annotation.gtf | grep -o gene_id | head -3如果第二个命令能正常输出gene_id基本说明格式是标准的。2.3 一个小习惯先跑单个样本的小规模验证在正式对所有样本跑HTSeq-count之前我强烈建议先挑一个样本只取前100万行比对记录跑一次验证命令能完整执行且输出合理。这个习惯帮我省下了无数次批量跑完才发现所有结果都是废的的惨痛教训。具体做法是先把BAM文件的头部保留然后截取一部分比对行samtools view -H sample.sorted.bam header.sam samtools view sample.sorted.bam | head -1000000 temp.sam cat header.sam subset.bam.tmp cat temp.sam subset.bam.tmp samtools sort -o subset.sorted.bam subset.bam.tmp不过这个做法比较绕。更简洁的办法是直接用HTSeq-count处理完整文件但只取前几行或者干脆先跑一个10万行的随机抽取。核心目的不是要一个严谨的小样本而是提前触发那些和格式、排序、注释解析有关的问题。用一个小样本验证完了再放心跑完整批次花不了几分钟。3. 六个核心参数逐个拆解每个都是如果不理解就会翻车的类型3.1-f、-r、-a三个决定文件怎么被读的参数HTSeq-count的完整用法是htseq-count [options] alignment_file annotation_file逐个拆解最关键的参数先把只影响输入读取的三个说清楚。-f或--format指定输入比对文件格式后面跟sam或bam默认是sam。如果你手里是BAM但忘了加-f bam工具会直接报格式解析错误。这个参数本身没有难度但有一个使用习惯问题很多同学写命令时图省事不加-f bam误以为工具会根据文件后缀自动识别——实际上HTSeq-count不会管后缀名它只认参数。我甚至在某个版本上试过后缀是.bam但没加-f bam工具直接按sam格式读产生的报错信息让人一头雾水。-r或--order决定文件内的read排序方式pos或name默认是name。前面已经说了坐标排序的BAM配合-r pos是标准选择。这里再补充一个容易混淆的点如果你的测序数据是单端排序方式并不影响计数逻辑但HTSeq-count仍然会读取这个参数走不同的内部排序归一化逻辑所以惯性用-r pos就好。-a或--minaqual设定最低MAPQ阈值默认10。这个参数的处理很简单低于阈值的比对记录会进入输出文件里的__too_low_aQual行。这个特殊行是评估比对质量的好指标后面讲输出文件时我会专门说怎么看它的占比来判断batch间的差异。3.2-s链特异性参数最容易一错到底、还不报错的坑-s或--stranded是HTSeq-count里最值得单独拉出来讲半天的参数因为它的默认值是yes而很多人并不知道这一点也不知道自己手上的建库试剂盒到底属于哪种链特异性类型。RNA-seq文库在构建过程中可以选择保留转录本的方向信息。经典的链特异性建库方案里第二链合成被专门去掉或做了特殊修饰于是测序得到的read方向能反映它原本是来源于哪条链。对这类文库计数时必须区分正负链对不保留方向信息的非链特异性文库直接按区域个数算。HTSeq-count里有三个值yes表示read方向与参考链相同才计数no表示不区分链方向reverse表示read方向与参考链相反时计数。这个相反很关键——大部分商品化的链特异性文库比如dUTP法测出来的read1实际上对应的是转录本的负链方向所以必须用reverse而不是yes。为什么说这是最坑的参数因为它错了完全不报错而且大多数情况下跑完所有样本、得到全部结果、甚至做完差异分析你都不会发现方向搞反了。直到某天你查看一个已知sense方向的基因发现它的表达量在所有样本里都接近零才意识到整体信号被算到了反义链上。实际运行过程中-s设错导致的结果是有义链基因的reads被打进__no_feature而反义链的reads被打进该基因。如果产出的counts配上有义链的标准基因列表差异倍数可能完全颠倒。怀疑链方向设错时除了翻看建库试剂盒说明书还有一个快速判断方法找同一样本的两个已知强表达基因一个sense、一个antisense方向看HTSeq-count在-s yes和-s reverse下的输出哪个更符合你的预期。或者用专门的工具推断链特异性但通常检查一两个基因就能得出结论。3.3-t、-i、-m决定把read数到谁头上的参数-t或--type指定要对哪种feature类型进行计数默认是exon。这意味着HTSeq-count只会提取GTF里第三列等于exon的行来构建基因模型。如果你的GTF文件是基于转录本水平定义的可能exon行本身没问题但如果你的注释文件用的feature类型叫transcript这个其实挺常见那计数时就要改成-t transcript。RNA-seq定量通常关注的是基因层面的表达量所以绝大多数场景保持exon是合理的但要知道这个参数的存在否则面对某些只定制了非标准GTF文件时无从下手。-i或--idattr指定用GTF第9列的哪个属性作为feature分组的ID默认是gene_id。它和-t配合使用-t exon把外显子行筛出来然后把gene_id相同的所有外显子合并成一个可计数区域。如果你希望统计的是转录本亚型的表达量比如做可变剪接分析时可以把-i改成transcript_id这样每个转录本会被当成独立的计数单元。但从RNA-seq常规差异表达分析看基因水平统计是主流因为亚型定量受read长度和比对准确性影响更大单独用HTSeq-count去算转录本水平并不是太合适后面我会提一句更合适的替代工具。-m或--mode指定overlap判断的模式默认union另外还有intersection-strict和intersection-nonempty。这三种模式的差异在于当一条read和某个基因的一个或多个外显子有部分重叠时如何判定是否计数。union的逻辑最宽容只要read和外显子区间有任何交集就算到该基因头上如果一条read同时重叠到多个基因则归为ambiguous。这个模式在绝大多数流程里都适用它被选为默认是有原因的——宽容的判定减少了因注释边界问题造成的读数流失代价是把部分位于基因间区的小片段也可能被统计进去。intersection-strict则要求read的完整区间必须包含在某一个单个exon内才计数更严格但会漏掉跨外显子连接的read所以不常用。intersection-nonempty介于两者之间。实际使用时除非你有强烈的理由要调整否则用union就好。曾经我在某个低表达基因的计数上对比过三种模式的差异union能多出约10%到15%的reads对检测灵敏度有明显影响。把核心参数整理成一张表方便对照设置参数含义常用值注意事项-f输入文件格式bam不根据后缀自动识别必须显式指定-r内部排序规则pos配-r pos之前必须先samtools sort-s链特异性模式no/reverse默认yes非链特异性文库要改no-a最低MAPQ10低于阈值的read进__too_low_aQual-t计数feature类型exon要符合GTF第三列定义-ifeature ID属性gene_id要存在且规范统一-moverlap判定模式union默认即可一般不用改4. 一次完整运行记录从安装到输出解读4.1 安装HTSeq的两种方式HTSeq是Python包安装最省事的方式是用conda直接建一个独立环境conda create -n htseq_env -c bioconda htseq python3.10 conda activate htseq_env如果你已经有了自己的Python环境也可以用pip安装pip install htseq严格讲HTSeq这个软件包内部包含多个工具命令行里的htseq-count是核心脚本。我习惯用conda管理是因为它会把依赖的numpy等一并装好避免版本冲突。曾经某次我直接pip安装了一个最新版结果和系统原有的numpy版本产生冲突导致所有用到numpy的脚本全部无法运行。conda环境隔离实际上能帮助你规避这类依赖混乱问题。安装完成后验证一下htseq-count --help能正常打印出参数说明就算装好了。4.2 一条实际跑通的完整命令下面是我在实际流程里用过的完整命令处理的是一个双端、链特异性建库的样本htseq-count \ -f bam \ -r pos \ -s reverse \ -t exon \ -i gene_id \ -m union \ --additional-attr gene_name \ sample.sorted.bam \ annotation.gtf sample.count这里面的--additional-attr gene_name是我额外加的它会在输出文件里增加一列gene_name。这个参数在HTSeq 0.12.4之后的版本可用作用是除了gene_id之外把GTF里的gene_name属性也输出到结果里省的后续做注释比对时再来一次映射。不过要注意如果GTF里有些基因没有gene_name属性输出会生成空值不碍事。整个运行过程在服务器上没有任何额外输出工具默认静默运行最后直接生成count文件。想实时看进度的话可以加-p或--progress参数默认是每10万条read输出一个进度点。4.3 输出文件长什么样五类特殊行怎么读HTSeq-count输出的是一个两列的文件第一列是基因ID第二列是原始read计数ENSG00000000003 1520 ENSG00000000005 0 ENSG00000000419 542 ... __no_feature 278431 __ambiguous 22108 __too_low_aQual 1392 __not_aligned 5367 __alignment_not_unique 88421前面是每个基因的计数最后五行是特殊统计行。千万别把这五行当成基因直接交给下游工具——有些流程脚本合并样本count时报错就是因为把__ambiguous这类行也当成了基因名。这些特殊行的意义如下__no_feature整条readfragment和任何基因的外显子都没有重叠。比例过高说明GTF注释不完整或链特异性方向设置错误。__ambiguousread同时和多个独立基因有重叠无法判定归属。这个比例通常不会太高但要注意在基因密度高的区域真实占比可能到5%到8%。__too_low_aQualMAPQ低于阈值被过滤。__not_aligned完全没有比对上的reads。__alignment_not_unique比对到基因组多个位置无法唯一确定的reads。多拷贝基因或重复区域的特异性问题会体现在这里。我判断数据质量的一个重要习惯是计算__ambiguous加上__no_feature在总reads中的占比。正常情况下两者合计应该在15%到25%左右取决于基因注释密度和read长度。如果超过35%优先怀疑GTF版本和比对参考不一致其次怀疑链特异性参数设反了。有一次我跑一个特殊样本时发现__no_feature高达60%排查后发现是GTF文件里只有主要染色体没有线粒体染色体和随机染色体而Ribo-seq数据里有大量线粒体reads把线粒体注释补进去才恢复正常。输出文件的基因间顺序是随GTF文件的读取顺序决定的不是排序后的。所以多样本比对时不要直接用paste拼接两个count文件作为表达矩阵因为不同样本的输出行顺序可能完全相同注释一样但你还是建议统一做一次order校验用sort -k1,1对文件内部按基因名排序后再合并。5. 实战中遇到的报错我把完整排查过程记录下来5.1 GTF解析失败第9列的格式问题第一个典型的报错场景是HTSeq-count在读取GTF文件时直接崩溃报错大致是说attribute ( is not valid here或者类似的解析错误。当时我的GTF文件是从某网站直接下载的用less打开看起来很正常。但真正的问题出在第9列的结尾——GTF标准要求每个键值对以分号结束最后一个属性后面也必须有分号。那个文件里部分行最后一个属性后面没有分号许多宽松的解析器能容忍但HTSeq-count对格式的要求比较严格直接判定解析失败。排查过程很简单我用grep -v尝试找到问题行awk -F \t {if ($3exon) print $9} annotation.gtf | grep -v ;$ | head正常GTF文件的每一行属性列都应该以分号结尾会查出所有不合规的行。我当时发现是某条染色体上所有外显子行的第9列都缺少结尾分号应该是文件在切割或转换过程中被截断了一列。修复方法是补全分号或者直接换成更规范的GTF来源。我的建议是不要自己写脚本去修补格式因为GTF文件的坑往往不止一个。重新下载一份来自官方注释体系、且和参考基因组版本严格配套的GTF远比修补一份来路不明的文件更安全。5.2 染色体命名不一致导致的全军覆没第二个常见的报错场景是HTSeq-count运行了一两秒就结束输出文件里所有基因计数都是0全部reads都被归到__no_feature里。这个场景的排查链路值得写下来因为一开始非常迷惑。我当时先怀疑链特异性设错了但改了-s no和-s reverse之后结果毫无变化。接着怀疑GTF有问题但单独看GTF文件里面明明有几十万个外显子记录。最后我用samtools view抓了几条比对记录再看GTF的染色体名才意识到问题所在samtools view sample.sorted.bam | head -1 | cut -f3输出的是Chr1而GTF文件里写的是1。比对时用的参考序列带上了Chr前缀而GTF里用的是纯数字染色体名。HTSeq-count做overlap判断时染色体名不一致就等于两个完全不同的世界所有read都找不到对应染色体上的基因。处理方式是把GTF里的染色体名统一成和BAM一致。当时用的是基于BAM文件染色体列表做个批量替换的思路samtools view -H sample.sorted.bam | grep ^SQ | cut -f2 | sed s/SN:// awk -F \t {$1Chr$1; print} OFS\t annotation.gtf annotation.fixed.gtf这样把所有染色体名都加上前缀匹配BAM里的命名。在这个问题上花了半小时其实如果最开始就建立了比对参考和GTF必须同源的习惯完全能避免。这个教训在我后面的每次流程搭建中都是必修课。5.3 排序方式引发的报错和误区第三个典型报错是当我直接拿某工具输出的原始比对结果喂给HTSeq-count时报错信息会提到read name order和BAM的sort order不一致。这里的本质是某些比对工具默认输出的BAM就是name-sorted或者没有特定的排序保证。HTSeq-count在-r pos下按坐标读取和合并pair遇到一个pair的mate不在预期位置时它无法把两条read合并成fragment于是处理逻辑可能崩掉。这个场景的排查相对直白看到报错先检查排序状态samtools stats sample.bam | grep bam_sort更简单的判断方法在IGV里加载BAM文件看reads是否按坐标整齐排列。不看IGV的话直接一条命令检查samtools view sample.bam | awk {print $3} | uniq | head如果染色体名交替出现说明文件不是坐标排序的。解决办法就是跑一遍samtools sort再继续。这也是为什么我在前面的章节里一直强调任何工具输入BAM前都先按坐标排序这条铁律。5.4 批量运行时内存被打爆的经验最后一个经验算不上报错但值得记录批量处理几十个样本时如果每个样本的BAM文件特别大并且你用了-r name模式NGS数据是name排序时HTSeq-count内部要把所有read的坐标信息缓存起来等待mate配对内存占用会随文件大小快速膨胀把服务器搞到OOM。我遇到过的情况是一个约80GB的BAM文件用-r pos处理时内存占用稳定在几个GB以内改用-r name后内存直接冲上60GB。原因是pos模式下同一条染色体上的归档数据可以逐段释放而name模式下所有未配对的read必须一直驻留内存直到找到mate。所以不管是单样本还是批量跑老老实实统一samtools sort之后用-r pos既是流程规范的保证也是服务器内存的救星。6. 从单样本计数到批量表达矩阵下游衔接的几个关键细节6.1 批量处理多个样本的bash循环RNA-seq分析很少只跑一个样本。拿到一批样本之后写一个循环批量生成count文件是常规操作cat samples.txt | while read sample; do echo Processing $sample ... htseq-count \ -f bam \ -r pos \ -s reverse \ -t exon \ -i gene_id \ -m union \ ${sample}.sorted.bam \ annotation.gtf ${sample}.count done wait这个循环里需要注意两点。第一是并行任务数量要和服务器核心数匹配直接无脑把所有样本同时丢出去极端情况下会把服务器IO打满反而不如串行快。我给自己的流程加了个简单控制用xargs -P 8限制并发数cat samples.txt | xargs -P 8 -I {} bash -c htseq-count -f bam -r pos -s reverse -t exon -i gene_id \ {}.sorted.bam annotation.gtf {}.count 第二是运行之前检查annotation.gtf的路径是否统一写对了。因为批量脚本里的一个路径错误会导致后续所有样本重新跑一遍——这个坑我踩过两次之后现在都会在批量循环前先跑一个样本确认命令本身没问题。6.2 合并表达矩阵的正确方式每个样本的count文件合到一起才能交给DESeq2或edgeR做差异表达分析。合并矩阵有几种方式我先说我踩过的坑。不要直接用paste -d \t sample1.count sample2.count然后不停地cut掉ID列拼接到一起因为HTSeq-count输出文件的基因顺序和GTF文件读取顺序一致理论上多个样本顺序相同但一旦某个样本的GTF解析出现了微小差异比如来源略有不同版本的注释文件行顺序就可能错位paste直接拼接出的矩阵就是错的而且很难发现。我自己的做法是每个样本先按基因ID排序再用paste合并。更稳妥的做法是用脚本按基因ID做外层joinimport pandas as pd base None for sample in [A, B, C]: df pd.read_csv(f{sample}.count, sep\t, headerNone, names[gene_id, sample]) df df[~df[gene_id].str.startswith(__)] if base is None: base df else: base base.merge(df, ongene_id, howouter) base.to_csv(count_matrix.csv, indexFalse)这段代码用startswith(__)排除了五个特殊统计行用merge按gene_id做完整连接避免出现因行顺序不同导致的错位。实际写脚本时如果样本数不多直接按这个思路去扩展就行。6.3 交给DESeq2之前一定要做的完整性检查合并后的count矩阵在进入DESeq2之前有几步非常快的完整性检查做了能省下后面很多和样本聚类相关的debug时间。首先检查矩阵的行数和GTF里的基因数是否大致一致。一个标准的人类GTF大约有两万个蛋白编码基因加上各种非编码RNA总计数在6万个左右如果你的矩阵里只有几千行大概率是某个样本在计数时把-i参数设成了transcript_id输出的单位是转录本而不是基因行数会多出不少。其次看每个样本的__no_feature占比是否稳定。RNA-seq样本之间的这个比例只要在合理范围内就不会引起注意但如果某个样本明显偏高大量reads散落在外显子区域之外可能是该样本的比对质量差或者文库被污染跑DESeq2之前最好先把这个样本标记出来观察。最后一步把count矩阵的列名改干净再进RsuppressPackageStartupMessages(library(DESeq2)) countData - read.csv(count_matrix.csv, row.names 1) countData - round(countData) condition - factor(c(ctrl, ctrl, treat, treat)) colData - data.frame(row.names colnames(countData), condition) dds - DESeqDataSetFromMatrix(countData, colData, design ~ condition)这里注意count值必须是整数HTSeq-count的输出本来就是整数但如果你做了归一化或者any后续处理一定要还原成整数再交给DESeq2否则DESeqDataSetFromMatrix会报错。6.4 HTSeq-count不是万能的AREA扩展了解一下写到这里特别想说一下HTSeq-count的定位。它适合做基因水平的计数稳健、可解释、结果可复现这是它至今被众多标准化流程选用的根本原因。但如果你要做转录本亚型水平的定量HTSeq-count直接按transcript_id去计数并不理想因为一个转录本的外显子组合包含了多个剪接异构体单纯把reads分到某个转录本的区间里无法处理reads同时来自多个亚型的情况。这种情况可以考虑其他方案比如直接从比对结果做转录本定量的工具不再赘述具体名字了或者如果坚持用HTSeq体系一些扩展版本在输入BAM和GTF后可以产出亚型层面的估计值。我的习惯是基因水平分析固定用HTSeq-count亚型水平分析用专门的转录本定量工具不硬拿一个工具去解决所有问题。7. 系列第11篇之前几个值得记住的经验前面把HTSeq-count从原理到实操完整过了一遍最后整理几条我在实际批次处理中反复验证过的心得希望对正在跑流程的朋友有帮助。第一参数一旦调通不要轻易改动并把它固化到流程管理文件里。HTSeq-count本身是个非常稳定的工具真正的变化往往在环境、GTF版本、上游比对工具参数这些环节。跑一批新数据时优先沿用上一批完全相同的HTSeq-count参数不要因为看到某个教程推荐了不同参数就顺手改掉。RNA-seq分析里流程的可重复性比某个参数的一点优化重要得多。第二在输出count文件的同时把每个样本的比对统计指标也一并留档。包括samtools flagstat的比对率、__no_feature比例、__ambiguous比例。这几个数字连同最后的表达矩阵一起保存遇到差异分析结果不符合预期时顺着这些指标能很快定位到是建库问题、比对问题还是计数问题。第三给小样本做的参数验证在批量跑之前一定再做一次完整文件的验证。小样本能提前暴露格式和兼容性问题但有些和内存、排序相关的问题只在完整文件上才会出现。批量跑之前先完整跑一个样本确保时间和服务器资源不被浪费。这篇文章的核心价值就是让HTSeq-count的使用者少走我走过的弯路。从输入文件的规范到参数的语义从报错的排查思路到下游矩阵的衔接每一步都有明确的为什么。读完你不需要再反复查手册去理解一个参数的含义——真正上手跑一个样本把这些细节过一遍后面的整个流程就会顺畅很多。