之不可乎骤得 托遗响于悲风❝RNA-seq转录组分析是生信分析的高频需求但手动敲命令不仅繁琐易出错还难复现、难扩展。今天分享一套基于Snakemake的自动化RNA-seq流程涵盖「质控→比对→定量」全核心步骤代码可直接复用新手也能快速上手conda环境搭建# 创建一个conda环境 conda create -n rnaconda createconda的核心命令用于创建新的虚拟环境-n rna-n是--name的简写指定环境名称为rna该行执行后会提示确认安装基础依赖后续加-y可自动确认。# 激活新创建的环境 conda activate rnaconda activate激活指定的conda环境rna要激活的环境名激活后终端前缀会变成(rna)表示后续命令都在该环境中执行。#安装所需的软件 conda install fastp hisat2 samtools subread r-argparse r-jsonlite r-dplyr r-tidyr r-tibble bioconductor-rsubread -yconda install在当前激活的环境中安装软件后续是要安装的软件列表fastp测序数据质控工具hisat2RNA-seq序列比对工具samtools处理SAM/BAM格式文件的核心工具subread包含featureCounts基因表达定量工具r-argparse/r-jsonlite/r-dplyr/r-tidyr/r-tibble分析所需的R包参数解析、数据处理bioconductor-rsubreadR版本的subread用于调用featureCounts-y自动确认所有安装提示无需手动敲yes。Snakemake脚本解析Snakemake 是流程管理工具核心逻辑是「定义目标文件→自动推导执行步骤」分模块实现功能1. 全局配置与样本管理import pandas as pd导入Python的pandas库命名为pd作用用于读取样本列表文件处理表格数据。# 指定配置文件参考基因组路径等全局参数 configfile: config.yamlconfigfile:Snakemake的关键字用于指定外部配置文件路径config.yaml配置文件名称YAML格式里面存放不随样本变化的全局参数比如参考基因组路径、链特异性类型、样本列表路径后续所有规则都能通过config[参数名]调用避免硬编码。# 从配置文件指定的样本列表文件中读取样本名生成样本列表SAMPLES SAMPLES pd.read_csv(config[samples_file], sep\t, headerNone, usecols[1]).squeeze(columns).tolist()这行是读取样本列表并提取样本名sep\t文件分隔符为制表符样本列表通常是TSV格式headerNone文件无表头如果样本列表有表头需改为header0usecols[1]只读取第2列索引从0开始假设样本名存在第2列config[samples_file]调用配置文件中samples_file参数样本列表文件的路径比如sample_list.txtpd.read_csv(...)pandas读取文本文件参数说明.squeeze(columns)将读取的单列DataFrame转换为Series简化数据结构.tolist()将Series转换为Python列表最终SAMPLES是所有样本名的列表比如[sample1, sample2, sample3]。print(SAMPLES)打印SAMPLES列表作用调试用运行流程时可确认样本名是否正确读取新手建议保留正式运行可注释。rule all:ruleSnakemake的核心关键字用于定义规则一个规则对应一个分析步骤all规则名自定义通常用all表示“流程最终目标”作用rule all是流程的“总目标”Snakemake会反向推导要生成rule all的输入文件需要执行哪些规则。input:input:规则的关键字定义该规则需要的输入文件对于rule allinput就是流程最终要生成的所有文件只要这些文件不存在/有更新Snakemake就会执行对应规则。expand(03.aligned/{sample}.bam, sampleSAMPLES), # 所有样本比对BAM文件expand()Snakemake的核心函数用于批量生成文件路径03.aligned/{sample}.bam文件路径模板{sample}是“通配符”sampleSAMPLES用SAMPLES列表中的每个样本名替换通配符举例如果SAMPLES [s1, s2]会生成[03.aligned/s1.bam, 03.aligned/s2.bam]04.quanti/genes.counts.matrix, # 基因表达矩阵指定流程最终要生成的“基因count表达矩阵”文件路径该文件是所有样本的基因计数合并后的表格用于后续差异分析。04.quanti/genes.TMM.TPM.matrix, # TPM标准化矩阵指定“TPM标准化表达矩阵”文件路径TPM是RNA-seq中常用的表达量标准化方法消除测序深度和基因长度的影响。01.clean_data/fastp_summary.tsv, # fastp质控汇总报告指定“所有样本的fastp质控结果汇总表”路径包含每个样本的原始读长、过滤后读长、Q30比例等关键质控指标。03.aligned/hisat2_summary.tsv # hisat2比对汇总报告指定“所有样本的hisat2比对结果汇总表”路径包含每个样本的总读长、比对率、唯一比对率等关键比对指标。2. 模块1测序数据质控trim_fastq规则rule trim_fastq:定义名为trim_fastq的规则作用是“对单个样本的原始fastq文件做质控”。input: # 输入原始双端fastq文件00.rawdata/input:该规则的输入文件即“需要质控的原始数据”注释说明输入文件的位置00.rawdata/目录。r100.rawdata/{sample}_R1.fq.gz,r1给输入文件起“别名”方便后续引用00.rawdata/{sample}_R1.fq.gz原始数据的R1端文件路径双端测序的正向读长{sample}通配符Snakemake会自动匹配SAMPLES列表中的每个样本名。r200.rawdata/{sample}_R2.fq.gzr2R2端文件别名反向读长路径对应原始数据的R2端文件双端测序必须同时输入R1和R2。output: # 输出质控后fastq fastp报告output:该规则的输出文件即“质控后生成的文件”。r101.clean_data/{sample}_R1.fastp.fq.gz,质控后的R1端文件路径存放在01.clean_data/目录后缀.fastp.fq.gz标识是fastp处理后的文件。r201.clean_data/{sample}_R2.fastp.fq.gz,质控后的R2端文件路径。report01.clean_data/{sample}.fastp.html,fastp生成的可视化质控报告HTML格式可直接用浏览器打开包含质控结果图表。json01.clean_data/{sample}.fastp.jsonfastp生成的结构化质控报告JSON格式方便后续用脚本汇总所有样本的结果。threads: 4 # 4线程threads:指定该规则运行时使用的CPU线程数4线程表示并行处理可根据服务器性能调整比如8、16线程越多速度越快但占用资源越多。log: 02.logs/{sample}.fastp.log # 日志文件log:指定该规则的日志文件路径作用记录fastp运行的所有输出信息包括警告、错误如果流程失败可查看日志排查问题1 {log} 21后续shell中会把所有输出重定向到日志文件。shell: # fastp命令质控过滤低质量序列shell:Snakemake的关键字用于执行shell命令。 fastp -i {input.r1} -I {input.r2} \ -o {output.r1} -O {output.r2} \ -h {output.report} -j {output.json} \ -w {threads} 1 {log} 21 fastp的核心命令fastp调用fastp工具-i {input.r1}指定输入的R1文件{input.r1}引用该规则input中r1的路径-I {input.r2}指定输入的R2文件注意是大写I区分R1/R2-o {output.r1}指定输出的R1文件-O {output.r2}指定输出的R2文件大写O-h {output.report}指定HTML报告路径-j {output.json}指定JSON报告路径-w {threads}指定使用的线程数引用该规则的threads参数1 {log}将命令的“标准输出stdout”重定向到日志文件21将命令的“标准错误stderr”重定向到标准输出最终也写入日志反斜杠\用于换行避免命令行过长不影响命令执行。3. 模块2汇总fastp质控结果combine_fastp_report规则rule combine_fastp_report:定义名为combine_fastp_report的规则作用是“汇总所有样本的fastp JSON报告为TSV表格”。input: # 输入所有样本的fastp json报告该规则的输入是“所有样本的fastp JSON报告”。expand(01.clean_data/{sample}.fastp.json, sampleSAMPLES)用expand()生成所有样本的JSON报告路径列表比如[01.clean_data/s1.fastp.json, 01.clean_data/s2.fastp.json]。output: # 输出所有样本的质控汇总表 01.clean_data/fastp_summary.tsv输出汇总后的TSV表格制表符分隔可直接用Excel/Excel打开。shell: unique_files$(echo {input} | tr \n | sort -u | tr \n ) Rscript scripts/combine_fastp_report.R -i 01.clean_data/ -o 01.clean_data/ unique_files$(echo {input} | tr \n | sort -u | tr \n )echo {input}输出所有输入的JSON文件路径空格分隔tr \n将空格替换为换行每行一个文件路径sort -u去重避免重复文件导致汇总错误tr \n 将换行换回空格unique_files$(...)将处理后的文件路径列表赋值给变量unique_files后续脚本可调用Rscript scripts/combine_fastp_report.R -i 01.clean_data/ -o 01.clean_data/Rscript调用R脚本scripts/combine_fastp_report.R自定义的R脚本路径需提前编写作用是解析JSON报告并汇总-i 01.clean_data/指定输入目录存放所有JSON报告-o 01.clean_data/指定输出目录生成汇总表。4. 模块3构建HISAT2索引build_hisat2_index规则rule build_hisat2_index:定义名为build_hisat2_index的规则作用是“为参考基因组构建HISAT2索引”比对前必须做只需运行一次。input: # 输入配置文件中指定的基因组fasta文件 genomeconfig[genome]genomeconfig[genome]调用配置文件中genome参数参考基因组FASTA文件的路径比如/ref/hg38/hg38.fa。output: # 输出HISAT2索引文件 idxexpand(03.aligned/genome.hisat2.{idx}.ht2, idxrange(1, 9))expand(...)HISAT2索引会生成8个文件后缀.ht2编号1-8用range(1,9)生成1-8的数字输出路径为03.aligned/genome.hisat2.1.ht2到03.aligned/genome.hisat2.8.ht2这是HISAT2索引的标准命名格式。log: 02.logs/build_hisat2_index.log记录构建索引的日志索引构建耗时较长日志可查看进度/错误。shell: # hisat2-build构建索引 hisat2-build {input.genome} 03.aligned/genome.hisat2 1{log} 21 hisat2-buildHISAT2的索引构建工具{input.genome}输入的参考基因组FASTA文件03.aligned/genome.hisat2索引文件的前缀最终生成的索引文件会在该前缀后加.1.ht2等1{log} 21将所有输出重定向到日志文件。5. 模块4HISAT2序列比对align_hisat2规则rule align_hisat2:定义名为align_hisat2的规则作用是“将质控后的reads比对到参考基因组”。input: # 输入质控后fastp HISAT2索引 r101.clean_data/{sample}_R1.fastp.fq.gz, r201.clean_data/{sample}_R2.fastp.fq.gz, idxexpand(03.aligned/genome.hisat2.{idx}.ht2, idxrange(1, 9))输入包含两部分质控后的R1/R2文件来自trim_fastq规则的输出HISAT2索引文件来自build_hisat2_index规则的输出。output: # 输出SAM文件 比对日志 samtemp(03.aligned/{sample}.sam), log03.aligned/{sample}.hisat2.logsamtemp(...)temp()是Snakemake的特殊函数标记该文件为“临时文件”作用SAM文件体积大比对后通常几十G后续会转换为BAM压缩格式Snakemake会在后续规则执行完成后自动删除该临时文件节省磁盘空间log记录单个样本的比对日志包含比对率、读长分布等信息。params: libtypeconfig[libtype]params:规则的“参数”非文件输入而是自定义参数libtypeconfig[libtype]调用配置文件中libtype参数链特异性类型比如RF/FRRNA-seq文库的关键属性。threads: 4比对使用4线程HISAT2比对是计算密集型线程越多越快。shell: hisat2 -x 03.aligned/genome.hisat2 -1 {input.r1} -2 {input.r2} \ --new-summary --rna-strandness RF \ -S {output.sam} \ 1 {output.log} 21 HISAT2比对命令hisat2调用HISAT2比对工具-x 03.aligned/genome.hisat2指定索引文件前缀无需写.1.ht2HISAT2会自动识别-1 {input.r1}指定R1端输入文件-2 {input.r2}指定R2端输入文件--new-summary生成新格式的比对汇总信息方便后续脚本解析--rna-strandness RF指定RNA-seq的链特异性RF是最常见的文库类型可根据实际情况改为FR-S {output.sam}指定输出的SAM文件路径1 {output.log} 21将比对的汇总信息/错误信息写入日志文件。6. 模块5汇总HISAT2比对结果combine_hisat2_report规则rule combine_hisat2_report:定义名为combine_hisat2_report的规则作用是“汇总所有样本的比对日志为TSV表格”。input: # 输入所有样本的hisat2比对日志 expand(03.aligned/{sample}.hisat2.log, sampleSAMPLES)输入是所有样本的比对日志文件列表。output: # 输出比对结果汇总表比对率等 03.aligned/hisat2_summary.tsv输出汇总后的比对结果表包含每个样本的总读长、比对数、比对率、唯一比对率等。shell: # 调用自定义Python脚本汇总报告 python scripts/combine_hisat2_report.py -i 03.aligned/ -o 03.aligned/ python scripts/combine_hisat2_report.py调用自定义Python脚本需提前编写-i 03.aligned/指定输入目录存放所有比对日志-o 03.aligned/指定输出目录生成汇总表。7. 模块6SAM转BAMsamtools_sort_idx规则rule samtools_sort_idx:定义名为samtools_sort_idx的规则作用是“将SAM文件排序并转换为BAM格式同时构建BAM索引”。input: # 输入SAM文件 sam03.aligned/{sample}.sam输入是align_hisat2规则生成的SAM文件临时文件。output: # 输出排序后的BAM BAM索引 bamprotected(03.aligned/{sample}.bam), idx03.aligned/{sample}.bam.baibamprotected(...)protected()是Snakemake的特殊函数标记该文件为“受保护文件”作用防止Snakemake误删除该文件BAM是后续定量的核心文件必须保留idx03.aligned/{sample}.bam.baiBAM文件的索引后缀.bai用于快速查询BAM文件中的特定区域。threads: 4samtools排序是IO密集型使用4个线程。log: 02.logs/{sample}.sort.log记录排序过程的日志排序耗时较长可查看进度。shell: # samtools排序建索引 samtools sort - {threads} -o {output.bam} {input.sam} 1{log} 21 samtools index {output.bam} 第一行samtools sortSAM排序并转BAM- {threads}指定线程数-o {output.bam}指定输出的BAM文件路径{input.sam}输入的SAM文件第二行samtools index {output.bam}为BAM文件构建索引生成.bai文件。8. 模块7链特异性参数转换自定义函数hisat2_strandness config[libtype]从配置文件中读取链特异性参数比如RF/FR赋值给变量hisat2_strandness。def convert_strandness(hisat2_strandness):定义名为convert_strandness的自定义函数作用HISAT2的链特异性参数RF/FR和featureCounts的参数1/2格式不同需要转换。if hisat2_strandness FR: return 1如果HISAT2的参数是FR返回1对应featureCounts的-s 1。elif hisat2_strandness RF: return 2如果HISAT2的参数是RF返回2对应featureCounts的-s 2。else: raise ValueError(Invalid HISAT2 strandness value)如果参数不是FR/RF抛出错误避免后续定量出错。9. 模块8featureCounts定量quanti_featureCounts规则rule quanti_featureCounts:定义名为quanti_featureCounts的规则作用是“基于BAM文件计算基因的read count”。input: # 输入排序后的BAM 基因组gtf注释文件 bam03.aligned/{sample}.bam, gtfconfig[gtf]输入包含排序后的BAM文件来自samtools_sort_idx规则配置文件中的GTF注释文件基因结构注释比如/ref/hg38/hg38.gtf。output: # 输出每个样本的基因表达量文件 expr04.quanti/{sample}.count, log04.quanti/{sample}.logexpr单个样本的基因count文件每行一个基因包含count数log定量过程的日志记录定量参数、基因数等。params: strandSpecificconvert_strandness(config[libtype]), prefix04.quanti/{sample}strandSpecificconvert_strandness(...)调用自定义函数转换链特异性参数prefix输出文件的前缀featureCounts会生成多个文件前缀统一。shell: Rscript software/RunFeatureCounts/featurecounts.R \ -b {input.bam} -g {input.gtf} -s {params.strandSpecific} \ -o {params.prefix} Rscript software/RunFeatureCounts/featurecounts.R调用R版本的featureCounts脚本rsubread包-b {input.bam}指定输入的BAM文件-g {input.gtf}指定GTF注释文件-s {params.strandSpecific}指定链特异性参数1/2-o {params.prefix}指定输出前缀。10. 模块9合并表达矩阵abundance_estimates_to_matrix规则rule abundance_estimates_to_matrix:定义名为abundance_estimates_to_matrix的规则作用是“将所有样本的count文件合并为表达矩阵”。input: # 输入所有样本的count文件 exprsexpand(04.quanti/{sample}.count, sampleSAMPLES)输入是所有样本的count文件列表。output: # 输出合并后的计数矩阵 TPM矩阵 counts_mat04.quanti/genes.counts.matrix, tpm_mat04.quanti/genes.TMM.TPM.matrix输出两个矩阵genes.counts.matrix原始count矩阵行基因列样本genes.TMM.TPM.matrixTPM标准化矩阵TMM是边缘修剪均值标准化方法。log: 02.logs/abundance_estimates_to_matrix.log记录矩阵合并过程的日志。run:run:Snakemake的关键字用于执行Python代码区别于shell:执行shell命令。unique_exprs list(set(input.exprs))input.exprs是所有样本的count文件列表set()去重避免重复文件再转回列表。with open(quant_files.txt, w) as f: for expr in unique_exprs: f.write(f{expr}\n)打开quant_files.txt文件写入模式将去重后的count文件路径逐行写入作用后续Perl脚本需要读取该文件获取所有count文件的路径。shell( perl software/RunFeatureCounts/abundance_estimates_to_matrix.pl \ --est_method featureCounts \ --out_prefix 04.quanti/genes \ --quant_files quant_files.txt 1{log} 21 sed -i 1s/^/gene_id/ {output.counts_mat} sed -i 1s/^/gene_id/ {output.tpm_mat} )perl software/RunFeatureCounts/abundance_estimates_to_matrix.pl调用Perl脚本用于合并count文件--est_method featureCounts指定定量方法为featureCounts--out_prefix 04.quanti/genes输出矩阵的前缀最终生成genes.counts.matrix和genes.TMM.TPM.matrix--quant_files quant_files.txt指定包含count文件路径的文本文件sed -i 1s/^/gene_id/ {output.counts_mat}sed文本处理工具-i直接修改文件1s/^/gene_id/在第一行开头插入gene_id为矩阵添加“基因ID”表头方便后续分析sed -i 1s/^/gene_id/ {output.tpm_mat}给TPM矩阵也添加gene_id表头。关键点configfile:调用外部配置文件expand()批量生成文件路径temp()/protected()管理文件生命周期核心分析步骤fastp质控→HISAT2建索引→HISAT2比对→SAM转BAM→featureCounts定量→合并表达矩阵链特异性参数转换RF/FR→1/2是RNA-seq定量的关键需与文库类型匹配。