从零开始构建转录组分析流程:Hisat2-Stringtie实战指南

发布时间:2026/5/19 22:29:23

从零开始构建转录组分析流程:Hisat2-Stringtie实战指南 1. 转录组分析流程概述转录组分析是研究基因表达的重要技术手段它能帮助我们了解细胞在特定条件下的基因活动情况。对于刚接触生物信息学的同学来说构建一个完整的分析流程可能会觉得无从下手。今天我就来分享一个经过实战检验的流程Hisat2-Stringtie组合方案。这个流程特别适合处理真核生物的RNA-seq数据从原始测序数据开始经过质量控制、序列比对、转录本组装最终获得基因表达量FPKM值。我去年在实验室处理一批番茄胁迫实验数据时就是用的这套方法效果相当稳定。整个流程可以分为几个关键阶段首先是数据获取和质控然后是基因组比对接着是转录本重构最后是表达量计算。每个阶段都有对应的工具来完成特定任务比如Hisat2负责高效比对Stringtie擅长转录本组装。下面我会详细介绍每个环节的具体操作。2. 环境准备与软件安装2.1 基础环境配置我强烈建议使用conda来管理分析环境这能避免软件依赖冲突的问题。首先安装miniconda然后配置国内镜像源加速下载wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh创建专属的转录组分析环境conda create -n rna-seq python3.8 conda activate rna-seq2.2 核心软件安装在配置好环境后用mamba安装所需软件mamba是conda的加速版conda install -y mamba mamba install -y hisat2 stringtie samtools subread fastp fastqc这里解释下各个工具的作用Hisat2高效的序列比对工具Stringtie转录本组装和定量工具Samtools处理SAM/BAM格式文件Subread包含featureCounts用于基因计数Fastp/FastQC数据质控双雄我曾经遇到过软件版本不兼容的问题所以特别建议固定版本号。比如mamba install hisat22.2.1 stringtie2.2.13. 数据获取与预处理3.1 原始数据下载假设我们要分析公共数据库中的RNA-seq数据以NCBI的SRA数据为例。首先安装sra-toolsmamba install -y sra-tools下载数据时推荐使用aspera加速ascp -QT -l 300m -i ~/asperaweb_id_dsa.openssh era-faspfasp.sra.ebi.ac.uk:/vol1/srr/SRR692/009/SRR6929579 .3.2 数据质控与过滤拿到原始数据后先用fastp进行质控过滤fastp -i SRR6929579_1.fastq.gz -I SRR6929579_2.fastq.gz \ -o clean_1.fq.gz -O clean_2.fq.gz \ -q 20 -u 30 -n 5 -l 50 \ -w 16 -h report.html -j report.json关键参数说明-q 20质量阈值20-u 30最多允许30%的低质量碱基-n 5最多5个N碱基-l 50保留长度≥50的reads质控报告要用FastQC检查fastqc -t 8 clean_1.fq.gz clean_2.fq.gz -o qc_report/我处理过一批数据原始质量Q20只有85%经过过滤后提升到98%有效去除了接头污染和低质量reads。4. 序列比对与转录本组装4.1 基因组索引构建使用Hisat2前需要先构建索引。以拟南芥基因组为例wget ftp://ftp.ensemblgenomes.org/pub/plants/release-48/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz hisat2-build -p 16 Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz athal_index索引构建很耗时建议在服务器上运行。我曾经构建人类基因组索引16线程跑了近2小时。4.2 序列比对使用构建好的索引进行比对hisat2 -x athal_index -1 clean_1.fq.gz -2 clean_2.fq.gz \ --dta --rna-strandness RF \ -p 16 -S aligned.sam 2 align_stats.txt重要参数解析--dta为Stringtie优化输出--rna-strandness RF链特异性设置-p 16使用16线程比对完成后转换为BAM格式samtools sort - 16 -o aligned.sorted.bam aligned.sam samtools index aligned.sorted.bam4.3 转录本组装现在轮到Stringtie出场了stringtie aligned.sorted.bam -G Arabidopsis_thaliana.TAIR10.48.gtf \ -o transcripts.gtf -p 16 -l SRR6929579这里需要注意参考注释文件GTF需要与基因组版本匹配-l参数设置样本前缀建议保留默认的FPKM计算方式我曾经对比过不同组装参数的效果发现过度调整-c(min_read_coverage)参数反而会丢失低表达转录本。5. 表达定量与分析5.1 合并多个样本结果当有多个样本时需要先合并GTF文件ls *.gtf gtf_list.txt stringtie --merge -G reference.gtf -o merged.gtf gtf_list.txt5.2 基因计数使用featureCounts进行基因水平计数featureCounts -T 16 -p -t exon -g gene_id \ -a merged.gtf \ -o gene_counts.txt *.bam5.3 FPKM计算Stringtie已经计算了FPKM如果需要重新计算# CountToFPKM.pl脚本示例 while(IN){ chomp; Fsplit; $total$F[1]; } seek(IN,0,0); while(IN){ next if $.1; Fsplit; $fpkm$F[1]*1e9/($total*$F[5]); print $F[0]\t$fpkm\n; }实际分析时我通常会同时保留原始counts和FPKM值因为不同下游分析需求不同。DESeq2等差异表达分析更适合用counts数据而样本间比较则可以用FPKM。6. 常见问题排查6.1 比对率低怎么办如果Hisat2比对率低于70%可能需要检查基因组版本是否匹配测序数据是否有污染是否忘记设置链特异性参数我遇到过因为链特异性参数设置错误导致比对率只有40%的情况调整后提升到85%。6.2 Stringtie组装结果异常当发现组装出的转录本数量异常多时检查是否启用了--conservative模式适当提高-c参数值确认不是DNA污染导致有一次我将-c从默认的2.5提高到5转录本数量从6万降到了3万更接近预期。6.3 计算资源优化对于大型数据集Hisat2比对可以使用--parallel参数Stringtie的-p参数不要超过实际CPU核心数排序BAM文件时控制内存使用-m 4G在处理100个样本的大项目时我编写了批量提交脚本配合任务队列系统将总运行时间从一周缩短到两天。

相关新闻