企业官网建设流程全解析

简介：这是一份面向生物信息学初学者的RNA-seq转录组数据分析完整流程项目代码，覆盖从原始测序数据到功能富集的核心环节。内容以两个不同组织的六个样本为例，串联FastQC质控、Trimmomatic剪接、HISAT2比对、featureCounts定量、DESeq2差异表达分析，并包含GO与KEGG富集分析及差异基因火山图绘制方法。压缩包共3个文件，以可运行的inscode项目文件为主，附带了项目说明html及gitignore配置文件，整体约8KB，便于直接查看代码结构和启动演示。已有518人学习下载，适合希望系统掌握转录组上游分析流程并动手实操的入门学习者。资源中的数据经过特殊修改，可在学习过程中放心使用，帮助理解每一步工具参数设置、结果解读与常见分析思路，为后续独立开展类似项目打下基础。

1. RNA-seq转录组分析：为什么完整流程值得跑通

做转录组的人应该都遇到过这种场景：拿着测序公司返回的FASTQ文件，脑子里很清楚要的是差异表达基因列表，但真到动手时才发现，从原始数据到火山图之间隔着质控、过滤、比对、定量、差异分析、富集整整六道工序。我去年帮某课题组复现一套RNA-seq分析时，发现他们三个人分别跑了三套流程，比对软件一个用STAR一个用Hisat2，定量一个用featureCounts一个用salmon，最后差异基因列表对上的不到一半。问题不是工具不行，而是流程里每一步的选型和参数会层层传导，任何一处不严谨都会让结果变成一纸空谈。这份项目代码的价值就在于把一整套可靠的RNA-seq转录组分析流程固化成了可复现的脚本，新手拿过去能按顺序跑通，熟手也能在关键步骤上看到有据可依的参数配置。

2. 全流程选型与路径：从FASTQ到差异基因列表的六个环节

2.1 流程的六个环节与每步的产出物

RNA-seq分析无论数据来自人、小鼠还是其他物种，主干路径大致是固定的。搞清楚每个环节输入什么、输出什么、验证什么，比急着敲命令更重要。我按实际执行顺序把这条链路拆成六个环节：

我在实际复现这份项目代码时，它把前四个环节固化成了Shell脚本，第五和第六环节用R脚本承接。这个拆分思路我认为是对的：上游是重计算任务，适合写在Shell里丢给服务器跑；下游是交互式探索，放在R里调整阈值和绘图更方便。值得注意的一个细节是，每个环节之间通过中间文件衔接——过滤后的FASTQ、比对完的BAM、定量的count矩阵，这样任何一个环节出问题都可以只重跑该环节，不需要从头再来。

2.2 软件选型：为什么是STAR加featureCounts加DESeq2

这套流程选了STAR做比对、featureCounts做定量、DESeq2做差异分析，这个组合反映了目前最主流且相对稳定的技术路线。STAR的核心优势是剪接比对能力强，它把基因组切分成suffix array索引，搜索reads时能同时处理跨外显子的junction read，对真核生物转录组来说这是刚需。

featureCounts的优势体现在两点：一是速度极快，二是它对多映射read的处理策略是默认丢弃，这符合大多数差异分析对保守性的要求。DESeq2选它不是因为功能最花哨，而是它基于负二项分布建模，对生物学重复数较少的实验设计（比如三组重复）有比较稳健的表现，而且它对低表达基因的收缩估计能有效减少假阳性。

这套组合里需要预先确认的是STAR的索引构建需要参考基因组FASTA和GTF注释文件，featureCounts只接受BAM格式输入，DESeq2的输入需要是原始整数计数而非归一化后的值。如果手头数据来自链特异性建库，三条工具都要注意链方向参数的一致性，STAR里是--outSAMstrandField，featureCounts里是-s，DESeq2不直接感知链信息，但定量结果已经受前两者影响。这套流程把这几个参数的默认选择写在配置里，跑之前花十分钟确认一遍，比跑完再返工划算得多。

3. 上游实操：质量控制、STAR比对与featureCounts定量的参数细节

3.1 环境准备与数据质量控制：先用conda把环境固定住

做生信分析的人最怕的就是软件版本漂移。两个月前能跑的脚本，两个月后因为某个依赖库升级就报错，这是最常见的翻车方式。这套流程代码里用conda环境锁定了所有依赖的版本，我复现时按它提供的environment.yml直接建环境：

conda env create -f environment.yml conda activate rnaseq

environment.yml里把fastp、STAR、featureCounts、samtools、multiqc这些工具的精确版本号固定住了，这样做的好处是整个环境可复现，换一台机器也能得到一致结果。我一般还会在激活环境后运行一次conda list | grep -E "star|featurecounts|fastp"确认关键工具确实加载的是指定版本。

质控环节我重点看fastp参数。下面这个命令是这套流程里的核心过滤逻辑：

fastp \ -i sample_1.fastq.gz \ -I sample_2.fastq.gz \ -o clean_1.fastq.gz \ -O clean_2.fastq.gz \ --detect_adapter_for_pe \ --trim_poly_g \ --cut_front \ --cut_front_window_size 1 \ --cut_front_mean_quality 3 \ --cut_tail \ --cut_tail_window_size 1 \ --cut_tail_mean_quality 3 \ --cut_right \ --cut_right_window_size 4 \ --cut_right_mean_quality 15 \ --qualified_quality_phred 15 \ --length_required 36 \ --compression 6 \ --thread 8 \ --json sample.json \ --html sample.html

逻辑说明：--detect_adapter_for_pe让fastp自动识别并切除双端测序的接头序列，这比手动指定接头序列更省心，因为不同建库试剂盒的接头可能存在差异。--trim_poly_g针对NextSeq平台常见的polyG尾进行切除，如果测序平台是HiSeq或Novaseq，环境不同这条建议保留。三段式质量裁剪分别处理reads的头部、尾部和滑动窗口，之所以分成三部分是因为reads不同位置的碱基质量分布差异很大。

参数说明：qualified_quality_phred 15表示碱基质量低于Q15会被当作低质量处理，这是比较常见的阈值选择，如果数据整体质量很好可以放宽到Q20，但你会在reads保留率上付出代价。length_required 36过滤掉裁剪后长度不足36bp的reads，这个阈值跟后续STAR比对的最小匹配长度有关，设得太高会损失真实现有转录本的信息。跑完后我通常会看一眼HTML报告里的reads通过率和平均碱基质量，通过率低于70%的话建议回到原始数据去确认建库质量是否有问题。

3.2 STAR比对：索引构建与两遍比对的核心参数

STAR是这套流程里最吃内存也最需要细调参数的一步。第一步是构建参考基因组索引，代码如下：

STAR \ --runMode genomeGenerate \ --genomeDir ./star_index \ --genomeFastaFiles genome.fa \ --sjdbGTFfile genes.gtf \ --sjdbOverhang 149 \ --runThreadN 16 \ --genomeSAindexNbases 14 \ --limitGenomeGenerateRAM 31000000000

sjdbOverhang的值应该设置为读长减1，这套流程的测试数据是150bp双端测序，所以配置里写149。这个值影响剪接位点数据库的构建，如果读长变了必须跟着改。genomeSAindexNbases 14是基因组suffix array索引的规模参数，基因组越大这个值越大，人基因组默认14够用，基因组特别大时可以适当调大到15，但我不建议随便改，因为改得不好会显著增大内存占用。limitGenomeGenerateRAM设到31GB是为了确保索引构建过程不被内存限制中断。

比对阶段用的是两遍比对策略，第一遍生成剪接位点信息，第二遍利用这些信息做更精准的比对：

STAR \ --runMode alignReads \ --genomeDir ./star_index \ --readFilesIn clean_1.fastq.gz clean_2.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM SortedByCoordinate \ --outSAMunmapped Within \ --outSAMattributes Standard \ --outFilterMismatchNmax 2 \ --outFilterMultimapNmax 10 \ --outFilterScoreMinOverLread 0.33 \ --outFilterMatchNminOverLread 0.33 \ --alignSJoverhangMin 8 \ --alignSJDBoverhangMin 1 \ --runThreadN 16 \ --outFileNamePrefix sample_

逻辑说明：outSAMtype BAM SortedByCoordinate直接输出按坐标排序的BAM文件，省掉后续再单独跑一次samtools sort。readFilesCommand zcat告诉STAR输入文件是gzip压缩格式，需要先解压读取，如果不加这个参数STAR会把压缩文件当普通文本处理，导致比对率惨不忍睹。outFilterMismatchNmax 2限制每个read最多只能有2个错配，这个阈值对一致序列数据来说是合理的；如果样本是肿瘤或含有较高变异率的细胞系，可以放宽到3到4，但需要警惕错配容忍度过高带来的多定位问题。

参数说明：outFilterMultimapNmax 10允许一个read最多比对到基因组的10个位置，超过10个位置的read会被标记为多映射而丢弃，这个阈值实际上决定了最终比对率里包含还是排除重复序列区域的信号。alignSJoverhangMin 8要求剪接位点两侧至少各有8bp的比对锚定长度，这个值直接影响junction read的识别质量。两遍比对的好处是，第一遍把新发现的junction位点合并进索引，第二遍时这些位点已经可以作为已知剪接位点用于比对，对非模式物种或注释不完整的基因组来说差距尤其明显。

3.3 featureCounts定量：从BAM到计数矩阵

featureCounts这一步相对简单，但它决定了最终输入DESeq2的计数矩阵长什么样。命令里有个值得注意的细节是-s参数，它表示链特异性模式：

featureCounts \ -a genes.gtf \ -o counts.txt \ -t exon \ -g gene_id \ -p \ -s 2 \ -T 8 \ -B \ -C \ sample_1.bam sample_2.bam sample_3.bam

逻辑说明：-t exon -g gene_id指定了计数单元——exon是GTF注释文件第三列的特征类型，gene_id是第九列的属性名，意思是把所有比对到同一个基因所有外显子上的read累加起来作为该基因的表达量。-p表示输入的是paired-end数据，featureCounts会把成对的read当作一个片段来计数，避免同一片段的正反两条read被重复计算。-B要求片段的两端都能成功比对到参考基因组上，-C要求一个片段的两端比对位置一致，不允许一端比对到一条染色体、另一端比对到另一条染色体，这两个选项都用于提高计数质量。

参数说明：-s 2适用于stranded reverse文库（比如dUTP法建库），这是目前主流链特异性建库方式；如果数据是非链特异性文库，这里应该改成s 0，设错会导致一半的read计数丢失，表现是最终计数矩阵里大量基因表达量接近0。-T 8是线程数，这个参数和STAR的线程数保持一致即可。跑完后我建议看一眼counts.txt.summary文件，它会列出每个样本成功分配的read数量，如果所有样本的分配率都非常低，先回头检查BAM的染色体命名是否与GTF一致，这是最常见的静默错误来源。

4. 差异表达分析实操：DESeq2的输入、设计公式与结果解读

4.1 DESeq2的输入格式与样本设计

拿到featureCounts输出的计数矩阵后，接下来桌面切换到R环境。DESeq2对输入数据有一个硬性要求：必须是整数计数，不能是FPKM、TPM这些已经归一化过的值。这是因为DESeq2内部会对原始计数做自己的归一化处理，输入归一化值等于进行二次归一化，结果会失真。

读取数据并构建DESeq2对象的代码如下：

library(DESeq2) count_data <- read.table("counts.txt", header=TRUE, sep="\t", row.names=1, skip=1) count_matrix <- count_data[, 6:ncol(count_data)] colnames(count_matrix) <- c("ctrl_1", "ctrl_2", "ctrl_3", "treated_1", "treated_2", "treated_3") sample_info <- data.frame( row.names = colnames(count_matrix), condition = factor(c("ctrl", "ctrl", "ctrl", "treated", "treated", "treated"), levels = c("ctrl", "treated")) ) dds <- DESeqDataSetFromMatrix( countData = count_matrix, colData = sample_info, design = ~ condition )

逻辑说明：skip=1是因为featureCounts输出文件的第一行是命令参数信息，真正的内容从第二行开始。count_data[, 6:ncol(count_data)]的作用是把前五列注释信息（Geneid、Chr、Start、End、Strand）去掉，只留下各样本计数列。设计公式~ condition告诉DESeq2我们只关心处理组与对照组之间的差异，如果实验设计包含批次信息，可以在design里加上批次变量，比如写~ batch + condition。

参数说明：factor(c(...), levels = c("ctrl", "treated"))这一行的顺序决定了后续差异分析中的比较方向。DESeq2默认把factor的第一个水平作为参考组，这里ctrl排在第一个，所以最终结果里log2FoldChange的正负号表示处理组相对于对照组的变化方向。这个顺序非常关键，写反了会导致所有基因的上调和下调方向互换，而且不会产生任何报错，属于最危险的错误之一。

4.2 差异基因筛选标准与结果表的列解读

构建好DESeq2对象后，一行命令完成差异分析：

dds <- DESeq(dds) res <- results(dds, alpha = 0.05, lfcThreshold = 0) res_df <- as.data.frame(res) res_df$gene <- rownames(res_df) sig_genes <- subset(res_df, padj < 0.05 & abs(log2FoldChange) > 1)

逻辑说明：DESeq(dds)内部会依次完成大小因子估计、离散度估计和负二项分布模型拟合，这一步不需要人工干预，但如果样本数量太少会在此处报错，至少需要两个生物学重复。results()里设置alpha = 0.05表示显著性阈值，lfcThreshold = 0表示不设置表达变化倍数的最小门槛。筛选差异基因时我习惯在padj < 0.05之外加一个abs(log2FoldChange) > 1，即表达量变化超过两倍的标准，这个组合是当前文献中最常见的默认标准。

结果表中每一列的含义需要解释清楚：baseMean是该基因在所有样本中的平均归一化计数，log2FoldChange是变化倍数以2为底的对数，正值表示处理组上调，负值表示处理组下调。lfcSE是log2FoldChange的标准误，stat是Wald检验的统计量，pvalue是未校正的P值，padj是经过Benjamini-Hochberg方法校正后的P值。实际筛选时永远以padj为准而不是pvalue，因为当基因数量在两万左右时，单纯按pvalue < 0.05筛选会出来上千个假阳性基因。

排序和导出差异基因：

res_ordered <- sig_genes[order(sig_genes$padj), ] write.csv(res_ordered, file = "DE_genes_padj0.05_log2FC1.csv", quote = FALSE)

我通常会导出一份中英文注释信息附在后面，方便直接用Excel查看。还有个习惯是顺手跑一下plotMA(res, main="MA plot")看整体分布，如果图形呈对称喇叭状说明数据质量可靠，如果出现大面积偏移或形状异常，说明定量环节可能有问题。

5. 避坑排查：五个让分析重跑的常见问题与处理办法

5.1 参考基因组版本与GTF不匹配的沉默错误

现象：STAR比对率正常，featureCounts也正常跑完，但最终计数矩阵里大量基因的count值为0或异常低。 原因：这是我见到的频率最高的翻车现场。参考基因组FASTA和GTF注释文件来自不同版本或不同来源，比如基因组是ENSEMBLrelease 110的，GTF却是GENCODE的，两者对同一基因的ID命名方式不同，featureCounts按gene_id聚合时自然匹配不上。 解决：下载参照文件时把基因组和GTF的版本号记录下来，尽量从同一版本同一来源下载。用下面的命令快速核对染色体命名是否一致：

grep "^>" genome.fa | head -5 cut -f1 genes.gtf | sort -u | head -5

如果一边是chr1另一边是1，说明命名不兼容。需要一个脚本把GTF里的染色体名换成基因组里的写法，或者反过来。从那以后我拿到新参考文件第一件事就是跑这两行命令，花不了一分钟，但能省掉一整天返工。

5.2 双端测序比对率低的真实原因

现象：STAR比对率只有60%出头，看着明显偏低但又不完全算失败。 原因：这种情况最常见的原因是两条reads的序列方向信息在比对时没有被正确处理，或者reads里残留了未切除干净的接头序列。另一个容易忽视的原因是--outFilterMismatchNmax设置得太严格，对变异率较高的样本来说错配容忍度不够，导致大量reads因为两三个错配就被判定为未比对。 解决：先检查clean数据里接头残留是否已经清除，fastp的HTML报告里有明确指标。再把outFilterMismatchNmax从2放宽到4，如果比对率从60%提到85%以上，说明原数据变异率比预期高，放宽阈值是合理的。最后用samtools stats检查BAM文件里paired properly的比例，这个数字应该在90%以上才算正常。

5.3 STAR内存不足导致的任务卡死

现象：--runMode genomeGenerate或比对阶段运行到一半进程被OOM kill直接杀掉。 原因：STAR对内存的需求超出预期，尤其是构建索引时--genomeSAindexNbases设置过大，或者--limitGenomeGenerateRAM低于实际需要。人基因组索引构建通常需要30GB左右内存，而比对阶段每个线程也要消耗几GB。 解决：在提交任务时直接申请足够的内存，不要用服务器默认配置。如果基因组很大，把genomeSAindexNbases降到12到13也可以显著降低内存占用，代价是比对速度稍慢。还可以用samtools view -H检查一下BAM文件头部里的STAR版本信息，确认比对时用的索引和当前基因组目录一致。

5.4 featureCounts计数丢失：多映射与配对选项检查

现象：featureCounts跑完，日志显示某个样本成功分配的片段数只有二三十万，明显少于BAM文件里的总片段数。 原因：这种情况通常是一个或多个选项没设对。-s链方向设反会让一半reads因链信息不符被丢弃，-C选项会把那些一端比对到线粒体另一端比对到核基因组的嵌合片段全部排除，这些片段在很多样本里占比不小。 解决：先把-s 2改成-s 0跑一次，如果计数大幅回升说明是链方向问题。如果改了-s还是没改善，再检查BAM里比对到线粒体上的比例，用samtools idxstats sample.bam | grep MT看数值，线粒体reads比例超过5%时需要单独处理，常见做法是在比对前过滤掉线粒体序列或者专门加一条线粒体染色体参数。

5.5 DESeq2设计公式与因子水平的方向陷阱

现象：差异分析结果里log2FoldChange方向整体和实验预期相反。 原因：factor的levels顺序决定参考组，而参考组决定了变化方向。比如实验是处理组vs对照组，对照组必须是参考组，但很多人写代码时习惯把处理组写在前面，导致每个基因的上下调方向全部反转。这个错误不报错、不警告，如果不拿已知marker基因验证方向，很容易一路错到论文里。 解决：运行差异分析之前用下面两行代码主动确认参考组：

colData(dds)$condition <- relevel(colData(dds)$condition, ref = "ctrl") results(dds)$log2FoldChange

在写文章之前，选几个你明确知道该上调的基因（比如经典marker基因），看一眼log2FoldChange是不是正的。从那次之后我形成了个条件反射——每一组新的分析跑完DESeq2，先看参考组定义对不对再往下走，绝不动摇。

6. 进阶用法：把流程封装成脚本并验证批次效应

6.1 流程封装与参数配置分离

跑通一遍只是开始，真正能提升效率的是把整套流程封装成一份可用配置文件驱动的脚本。我复现这套流程代码后做的最重要一件事，是把所有变量抽出来放进config.sh，脚本主体只负责按顺序执行。配置文件大致长这样：

GENOME_DIR="/path/to/star_index" GTF_FILE="/path/to/genes.gtf" FASTQ_DIR="/path/to/fastq_files" OUTPUT_DIR="/path/to/output" THREADS=16 READ_LENGTH=150 STRANDNESS=2 REFERENCE_GROUP="ctrl" TREATMENT_GROUP="treated"

然后在主流程脚本里通过source config.sh方式导入这些变量，后续所有命令都用$THREADS、$STRANDNESS这样的变量代替硬编码。这样做的价值在于，换一批数据时只需要改一个配置文件，而不是钻进脚本里逐行替换路径。对批次效应明显的项目，我还会在config.sh里额外加两个字段：BATCH变量和CONDITION变量，DESeq2的design就可以写~ batch + condition，把批次信息纳入统计模型。

6.2 用PCA验证样本整体关系

拿到计数矩阵后不要急着跑差异分析，我现在的习惯是先做一次PCA快速看看全局。这段R代码是这套流程里验证批次效应最直接的方法：

library(DESeq2) vsd <- vst(dds, blind = TRUE) pca_data <- plotPCA(vsd, intgroup = c("condition"), returnData = TRUE) percentVar <- round(100 * attr(pca_data, "percentVar")) ggplot(pca_data, aes(PC1, PC2, color = condition)) + geom_point(size = 3) + xlab(paste0("PC1: ", percentVar[1], "% variance")) + ylab(paste0("PC2: ", percentVar[2], "% variance")) + theme_bw()

逻辑说明：vst()对计数矩阵做方差稳定变换，目的是让高表达基因和低表达基因的方差变得可比，这一步对PCA和样本间相关性计算很重要。blind = TRUE表示变换时不考虑分组信息，避免在探索性分析时引入实验设计的先验影响。plotPCA返回PC1和PC2坐标，把它们画出来就能直观看出各组的聚类情况。

判断维度有三条：第一，同一组的样本应该聚在一起，如果对照组里有一个样本明显游离在外，说明它可能是批次异常样本，需要回到原始数据检查；第二，处理组和对照组之间沿着PC1方向分开，说明处理效应非常强，差异分析结果大概率干净可靠；第三，如果PCA图上样本完全按批次聚类而不是按处理组聚类，说明批次效应占据了主导地位，这种情况下DESeq2的design里必须加入batch项，否则结果会被批次效应污染。这是我从一次失败项目中吸取的教训——当时有一组RNA-seq数据由于建库时间不同，PCA图里两个批次分别聚成两团，起初没在意，后来差异基因列表里全是批次差异而非处理差异。从那以后我每次拿到数据都强制把PCA走一遍，确定样本关系正常才继续往下做，希望帮到你。

本文还有配套的精品资源，点击获取