1. 项目概述从“数数”到基因表达定量如果你刚接触高通量测序数据分析尤其是RNA-seq可能会被一堆“定量”、“计数”、“差异表达”这些词搞得晕头转向。别慌其实核心的第一步就是搞清楚我们测出来的海量序列片段到底有多少落在了我们感兴趣的基因上。这个过程专业上叫“基因表达定量”而今天要聊的featureCounts就是干这个活儿的“金牌计数员”。简单来说featureCounts是一个专门用来统计比对到基因组上各个特征比如基因、外显子的序列片段数量的工具。它隶属于Subread软件包由澳大利亚Walter and Eliza Hall Institute的杨云浩Yunshun Chen等人开发以其速度快、内存占用低、结果准确可靠而闻名已经成为许多生物信息学分析流程中的标准配置。无论是处理来自Illumina、BGI还是其他平台的测序数据只要你有序列比对后的SAM/BAM文件和一个基因注释文件GTF/GFFfeatureCounts就能帮你高效地完成计数矩阵的构建。对于新手而言直接看官方文档里一长串的参数可能会感到畏惧。这篇内容的目的就是带你从零开始手把手完成featureCounts的安装、配置并运行一个完整的计数流程。更重要的是我会结合自己多年处理真实数据的经验把那些官方文档里不会写、但实际操作中大概率会踩到的“坑”和报错信息掰开揉碎了讲清楚让你不仅能跑通代码更能理解每一步背后的逻辑遇到问题知道去哪儿找答案。2. 环境准备与工具安装在开始“数数”之前我们得先把“计数员”请到我们的工作环境里。featureCounts的安装非常灵活你可以通过多种方式获取它。2.1 安装方式选择与实操最推荐的方式是通过Conda进行安装。Conda是一个强大的包管理和环境管理工具它能帮你自动解决软件依赖问题避免“装了这个缺那个”的窘境。如果你主要做生物信息分析还没用上Conda我强烈建议你从现在开始使用它。首先你需要安装Miniconda或Anaconda。之后创建一个专门用于RNA-seq分析的环境是个好习惯这能保证环境的纯净避免不同软件包之间的版本冲突。# 创建一个名为rna-seq的新环境并指定python版本非必须但推荐 conda create -n rna-seq python3.9 # 激活这个环境 conda activate rna-seq # 在激活的环境中通过bioconda频道安装subreadfeatureCounts包含在其中 conda install -c bioconda subread安装完成后在终端输入featureCounts -v如果能看到版本号例如featureCounts v2.0.3恭喜你安装成功。注意Bioconda频道需要提前配置。如果你首次使用可能需要运行conda config --add channels bioconda和conda config --add channels conda-forge。另外网络状况可能导致下载缓慢或失败可以考虑配置国内镜像源。如果你无法使用Conda也可以选择编译安装。这种方式适合对系统控制要求高或者需要特定版本的情况。下载源码前往Subread项目的SourceForge页面或GitHub仓库下载最新的源码包通常是.tar.gz格式。解压与编译tar -zxvf subread-2.0.3-source.tar.gz cd subread-2.0.3-source/src make -f Makefile.Linux # 对于Linux系统 # 对于Mac可能是 make -f Makefile.MacOS编译成功后在bin目录下会生成可执行文件你可以将其路径加入系统的PATH环境变量或者直接使用绝对路径调用。第三种方式是通过包管理器例如在Ubuntu上你可以尝试sudo apt-get install subread但这种方式提供的版本可能不是最新的。2.2 核心输入文件准备安装好工具只是第一步接下来要准备featureCounts工作的“原材料”。主要有两类比对文件SAM/BAM这是你的测序数据经过比对软件如HISAT2, STAR, BWA与参考基因组比对后产生的文件。它记录了每一条测序读段read在基因组上的位置信息。BAM是SAM的二进制压缩格式更节省空间featureCounts两者都支持但处理BAM通常更快。注释文件GTF/GFF这是基因组的“说明书”以文本格式详细描述了基因组上各个特征基因、转录本、外显子、启动子等的坐标和属性。最常用的是GTF格式。你需要从Ensembl、NCBI或UCSC等数据库下载与你所用参考基因组版本完全一致的GTF文件。版本一致性至关重要用GRCh38的基因组比对结果却搭配GRCh37的注释文件结果肯定是混乱的。假设你的工作目录结构如下这是一个清晰的好习惯/project/ ├── data/ │ ├── sample1.sorted.bam │ ├── sample2.sorted.bam │ └── ... (其他样本的BAM文件) ├── ref/ │ ├── genome.fa │ └── annotation.gtf └── scripts/ └── run_featureCounts.sh3. 核心参数解析与基础命令运行现在工具和材料齐备我们可以开始运行featureCounts了。一个最基础、最常用的命令长这样featureCounts \ -T 4 \ # 使用4个CPU线程 -t exon \ # 对GTF中的“exon”特征进行计数 -g gene_id \ # 使用“gene_id”属性将计数汇总到基因水平 -a /path/to/annotation.gtf \ # 注释文件路径 -o counts.txt \ # 输出计数矩阵文件 /path/to/sample1.bam /path/to/sample2.bam # 输入BAM文件让我们拆解这个命令里的关键参数理解它们为什么这么设置-T指定线程数。多线程能极大提升处理速度通常设置为可用CPU核心数。但要注意如果同时运行多个任务不要占满所有核心。-t指定要对注释文件中的哪类特征feature进行计数。对于标准的基因表达定量我们关心的是外显子exon因为成熟的mRNA是由外显子拼接而成的测序读段只有比对到外显子区域才可能来自该基因的转录本。-g指定用于将计数汇总aggregate的属性。一个外显子可能属于多个转录本但通常我们最终需要的是基因水平的表达量。gene_id属性就是用来告诉featureCounts“把属于同一个gene_id的所有外显子的计数加起来作为这个基因的计数”。如果你需要转录本水平的计数可以设置-g transcript_id。-a注释文件路径。-o输出文件路径。输出的counts.txt是一个制表符分隔的文本文件可以直接用Excel或R打开查看。运行成功后你会得到两个主要文件counts.txt计数矩阵和counts.txt.summary计数统计摘要。摘要文件非常有用它快速告诉你每个样本总共有多少比对上的读段其中有多少被成功分配给了基因。4. 高级功能与参数深度优化掌握了基础命令你已经能完成90%的工作。但要想让计数结果更精准适应更复杂的分析场景就需要了解下面这些高级参数和功能。4.1 处理双端测序与链特异性数据现在的RNA-seq实验绝大多数都是双端测序Paired-end并且很多是链特异性建库Strand-specific。featureCounts对此有专门参数进行优化。-p这个参数告诉featureCounts输入的是双端测序数据。开启后它会将一对读段mate pairs作为一个整体来考虑而不是两个独立的片段。这能提高计数的准确性特别是对于跨越外显子边界的读段。只要数据是双端的就一定要加-p。-s指定链特异性。这是新手最容易出错的地方之一。它的取值有0(默认)非链特异性。读段可以来自正义链或反义链。1链特异性且读段来自转录本链即与mRNA同向。例如常见的Illumina TruSeq Stranded Total RNA文库。2链特异性且读段来自转录本链的互补链。如何确定-s的值这完全取决于你的建库试剂盒和方案。你必须查阅实验记录或文库制备手册。设置错误会导致计数大幅减少甚至为零因为软件会在错误的链上寻找读段。一个处理链特异性双端数据的完整命令示例featureCounts -T 8 -p -s 1 -t exon -g gene_id \ -a annotation.gtf \ -o stranded_pe_counts.txt \ *.bam4.2 精细化计数策略允许多定位与分数计数高通量测序读段有时会比对到基因组的多个位置这被称为“多定位读段”Multi-mapping reads。此外一个读段可能同时覆盖多个基因例如重叠基因或反义转录本。featureCounts提供了精细的策略来处理这些“模糊”的读段。-M允许对多定位读段进行计数。默认情况下featureCounts会忽略不计数那些比对到多个位置的读段。加上-M后这些读段会被保留并计数。-O允许一个读段分配给多个特征如多个重叠的基因。默认情况下一个读段只会分配给一个特征-O关闭。--fraction当使用-M或-O时此参数启用分数计数。例如一个读段比对到3个位置那么每个位置分配到的计数是1/3而不是1。这比简单计数更符合生物学事实是推荐的做法。何时使用这些参数这取决于你的分析目的。如果你追求高精度的基因表达差异分析通常建议使用默认设置即忽略多定位和重叠读段因为它们的归属不确定引入噪声。但如果你在研究重复序列、假基因或需要最全面的计数时可以开启-M和--fraction。一个常见的折中方案是-M --fraction -O对多定位和重叠读段进行分数计数。4.3 定制化输出与元信息提取除了基础的计数矩阵featureCounts还能提供更多信息。-f在计数时统计特征如外显子水平的计数而不是直接汇总到基因水平。这通常与-O一起使用用于更精细的分析。--extraAttributes除了默认的gene_id你还可以在输出文件中包含注释文件里的其他属性列例如gene_name,gene_biotype等。这能让你的结果表格更易读。featureCounts -t exon -g gene_id -a annotation.gtf \ --extraAttributes gene_name,gene_biotype \ -o counts_with_attr.txt sample.bam-R输出详细的读段分配报告为每个输入BAM文件生成一个.featureCounts文件里面列出了每条读段被分配给了哪个基因。这个文件非常大主要用于深度调试和验证。5. 实战流程从原始BAM到计数矩阵让我们通过一个完整的实战案例将上述所有知识点串联起来。假设我们有6个样本的双端、链特异性RNA-seq数据已经用STAR比对并排序好了BAM文件。5.1 编写可复现的批处理脚本在scripts/目录下创建一个Shell脚本run_featureCounts.sh#!/bin/bash # 脚本run_featureCounts.sh # 描述使用featureCounts对多个BAM文件进行基因计数 # 设置工作目录和文件路径 PROJECT_DIR/path/to/your/project BAM_DIR${PROJECT_DIR}/data/bam_sorted GTF_FILE${PROJECT_DIR}/ref/GRCh38.gtf OUTPUT_DIR${PROJECT_DIR}/results/counts LOG_DIR${PROJECT_DIR}/logs # 创建输出目录和日志目录 mkdir -p ${OUTPUT_DIR} ${LOG_DIR} # 获取所有BAM文件列表按样本名排序 BAM_FILES($(ls ${BAM_DIR}/*.sorted.bam | sort)) echo 找到 ${#BAM_FILES[]} 个BAM文件待处理。 # 设置featureCounts参数 THREADS8 STRANDEDNESS1 # 根据你的建库方案修改1或2 FEATUREexon ATTRIBUTEgene_id EXTRA_ATTRgene_name,gene_biotype # 构建输出文件名 OUTPUT_MATRIX${OUTPUT_DIR}/gene_counts_matrix.txt SUMMARY_FILE${OUTPUT_DIR}/gene_counts_matrix.txt.summary LOG_FILE${LOG_DIR}/featureCounts_$(date %Y%m%d_%H%M%S).log # 执行featureCounts命令 echo 开始运行featureCounts... | tee -a ${LOG_FILE} echo 输入文件: ${BAM_FILES[]} | tee -a ${LOG_FILE} echo 输出文件: ${OUTPUT_MATRIX} | tee -a ${LOG_FILE} featureCounts \ -T ${THREADS} \ -p \ # 双端数据 -s ${STRANDEDNESS} \ # 链特异性 -t ${FEATURE} \ -g ${ATTRIBUTE} \ --extraAttributes ${EXTRA_ATTR} \ -a ${GTF_FILE} \ -o ${OUTPUT_MATRIX} \ ${BAM_FILES[]} \ 21 | tee -a ${LOG_FILE} # 将标准输出和错误都重定向到日志文件 # 检查命令是否成功执行 if [ $? -eq 0 ]; then echo featureCounts运行成功 | tee -a ${LOG_FILE} echo 计数矩阵: ${OUTPUT_MATRIX} | tee -a ${LOG_FILE} echo 统计摘要: ${SUMMARY_FILE} | tee -a ${LOG_FILE} # 可选查看一下摘要文件的前几行 echo 计数统计摘要前5行 | tee -a ${LOG_FILE} head -n 5 ${SUMMARY_FILE} | tee -a ${LOG_FILE} else echo featureCounts运行失败请检查日志: ${LOG_FILE} | tee -a ${LOG_FILE} exit 1 fi给脚本添加执行权限并运行chmod x scripts/run_featureCounts.sh ./scripts/run_featureCounts.sh。5.2 输出结果解读与后处理运行成功后我们重点查看两个文件计数矩阵 (gene_counts_matrix.txt)用head -n 10 gene_counts_matrix.txt查看前10行。文件开头会有几行以#开头的注释信息接着是表头和数据。第一列是基因ID随后各列是每个样本的原始计数。这个文件可以直接导入R中用于后续的差异表达分析例如使用DESeq2、edgeR等包。统计摘要 (gene_counts_matrix.txt.summary)这个文件没有表头第一列是状态后面各列是每个样本处于该状态的读段数。关键行包括Assigned: 成功分配给某个基因的读段数。这是我们最关心的比例越高说明数据质量越好通常期望在70%-90%以上。Unassigned_NoFeatures: 比对上了基因组但不在任何注释特征外显子区域的读段。可能是内含子、基因间区或注释不全。Unassigned_Ambiguity: 无法明确分配给单一基因的读段如重叠基因区域当不使用-O参数时这些读段会被归为此类。一个健康的样本Assigned的比例应该占主导。如果Unassigned_NoFeatures比例异常高可能需要检查注释文件版本是否与基因组匹配或者考虑是否包含了链特异性参数-s。后处理小技巧原始的计数矩阵文件包含注释行和长长的基因ID列有时我们需要一个“干净”的纯数字矩阵。可以用一些简单的命令处理# 提取纯计数矩阵跳过注释行取基因ID之后的所有列 tail -n 2 gene_counts_matrix.txt | cut -f 7- gene_counts_clean.txt # 将第一列基因ID作为行名保存为CSV格式方便R读取 tail -n 2 gene_counts_matrix.txt | awk BEGIN{OFS,} {print $1, $7, $8, $9, $10, $11, $12} gene_counts_for_r.csv6. 报错大赏与疑难排解实录即使命令看起来正确在实际运行中你仍可能遇到各种报错。下面是我总结的几个最常见的问题及其解决方法。6.1 常见报错代码与原因分析报错Error: cannot find the input file \xxx.bam\原因文件路径错误或文件不存在。这是最典型的“粗心错误”。解决使用绝对路径。在脚本中用$(realpath file.bam)获取绝对路径。用ls -lh确认文件是否存在且有读权限。报错Error: failed to find the gene identifier attribute in the provided GTF file原因GTF文件中找不到-g参数指定的属性默认是gene_id。不同来源的GTF文件其属性字段的命名可能略有差异。解决查看你的GTF文件前几行找到基因水平的标识符。使用head -n 5 your.gtf和grep -v \^#\ your.gtf | head -n 1 | tr \t \n | nl来查看第九列属性的格式。你可能需要将-g改为gene_name或其他标识符。报错Error: read is paired but no correct paired-end alignment is found for it(当使用-p参数时)原因BAM文件中的读段标记为双端但featureCounts无法找到其配对信息或者配对的读段没有正确比对。解决检查BAM文件是否经过正确排序和索引。使用samtools view -H your.bam查看文件头确认PG行中的比对软件是否正确。确保使用samtools sort按坐标排序并使用samtools index建立索引生成.bai文件。有时比对软件产生的BAM文件本身配对信息就有问题需要回顾比对步骤。报错Segmentation fault (core dumped)原因内存访问越界。可能由于软件bug、版本不兼容或者输入文件尤其是GTF格式异常、损坏。解决首先检查输入文件samtools quickcheck your.bam检查BAM文件grep -v \^#\ your.gtf | head检查GTF格式。尝试更新featureCounts到最新版本。如果GTF文件过大尝试使用-Q参数设置一个最小的比对质量阈值过滤掉低质量读段减少内存压力。警告Warning: 10% of reads are multi-mapping. Use -M to count them.原因这不是错误而是一个重要提示。有相当比例的读段是多定位的。解决根据你的分析目标决定。如果进行差异表达分析通常忽略此警告坚持默认设置不计数多定位读段。如果你需要统计所有可能来源可以加入-M --fraction参数但务必在报告中注明。6.2 结果异常排查清单有时候程序能跑完但结果看起来不对劲问题Assigned比例极低 50%检查1链特异性参数-s。这是头号嫌疑犯用错链方向会导致绝大部分读段被归为Unassigned_NoFeatures。务必核实建库方法。检查2GTF文件版本。确保GTF与比对所用的参考基因组版本完全一致例如都是GRCh38.p13。检查3特征类型-t。确保你计数的是exon而不是gene或CDS。检查4BAM文件内容。用samtools view your.bam | head -20看看读段是否真的比对到了染色体上参考序列名不是*以及比对质量如何。问题所有样本的计数几乎一样检查很可能你在命令中错误地多次指定了同一个BAM文件或者文件路径列表生成有误。检查你的输入文件列表。问题输出文件为空或只有表头检查输入BAM文件路径可能全部错误或者GTF文件路径错误导致程序无法开始计数。查看日志文件.log的前几行错误信息。6.3 性能优化与内存管理处理大型项目如上百个样本、全基因组测序时性能和内存是关键。速度优化使用-T设置多线程这是最有效的加速手段。输入BAM文件时尽量提供排序好的BAM文件列表而不是使用*.bam通配符这有助于内部优化。使用BAM格式而非SAM格式。内存管理featureCounts内存占用主要与GTF文件的复杂度和特征数量有关。如果内存不足可以尝试使用-Q过滤低质量比对。分批次运行每次处理部分样本最后合并结果featureCounts支持多次运行同一注释然后合并计数。确保服务器有足够的物理内存避免使用Swap分区否则会急剧变慢。7. 上下游流程衔接与最佳实践featureCounts不是分析的终点而是连接比对和差异表达的桥梁。这里分享一些与上下游工具衔接的经验。7.1 与比对工具的配合不同的比对工具STAR, HISAT2, TopHat2产生的BAM文件在格式上可能有些细微差别但featureCounts兼容性很好。关键在于排序务必使用samtools sort或比对工具自带的排序功能按基因组坐标排序coordinatesort。这是featureCounts高效工作的前提。索引排序后立即建立索引samtools index生成.bai文件。虽然featureCounts不强制要求索引但某些情况下如并行处理有索引会更稳健。标记重复在ChIP-seq分析中标记PCR重复使用samtools markdup或Picard是重要步骤。但在RNA-seq基因表达定量中通常不建议在featureCounts之前标记或去除重复因为高表达基因区域本身就会产生大量“重复”读段这是真实的生物学信号而非技术重复。7.2 为下游差异表达分析准备输入DESeq2和edgeR是两大主流的差异表达分析R包。它们需要的输入就是featureCounts产生的原始计数矩阵。对于DESeq2你需要准备一个数据框行是基因列是样本值是整数计数。就是我们在5.2节后处理中得到的“干净”矩阵。在R中使用read.table(“gene_counts_clean.txt”, headerTRUE, row.names1)即可读入。对于edgeR读入方式类似。它们都要求是原始计数raw counts不要使用FPKM、TPM等标准化后的值。标准化是这些包内部自己完成的。一个重要的准备工作是检查计数摘要。在将矩阵导入R之前养成习惯先看.summary文件。如果某个样本的Assigned比例远低于其他样本或者总比对读段数Total异常低这个样本可能在下游分析中被视为质量不佳的离群值需要考虑剔除。7.3 版本控制与可复现性生物信息学分析中可复现性至关重要。记录下你使用的每一个软件版本和关键参数。# 记录软件版本 featureCounts -v 21 | tee software_versions.log samtools --version software_versions.log # 将最终运行的命令完整保存 echo featureCounts -T 8 -p -s 1 -t exon -g gene_id -a ref.gtf -o counts.txt *.bam run_command.log更好的做法是使用Conda环境导出功能conda env export -n rna-seq environment.yaml。这样别人就能用这个文件完全重建你的分析环境。最后关于“零基础看懂运行代码”我想说生物信息学的代码并不可怕它们大多是由一系列简单的命令和参数组成的。遇到报错首先不要慌仔细阅读错误信息它通常已经指明了方向。从最简单的命令开始逐步添加参数每步都验证输出这是最稳妥的学习方式。featureCounts作为一个高效稳健的工具掌握它你的转录组数据分析之路就打下了一块坚实的基石。
原理、配置与工业级应用全解析)
