RNA-seq数据上游处理完整流程

原创

生信医道

发布于 2025-07-01 18:59:31

3710

文章被收录于专栏：RNA-seqRNA-seq Linux上游分析

生信医道

前言

今天卡卡在处理RNA-seq的数据，整理出来一份上游的处理流程，分享给大家。从环境搭建到最终的表达量定量，一步步带你掌握RNA-seq数据分析的核心技能。

环境准备

创建conda环境

# 创建名为rnaseq的conda环境（基于Python 3.8）
conda create -n rnaseq python=3.8 -y 

# 激活环境 
conda activate rnaseq

# 安装STAR比对工具
conda install -c bioconda star -y

# 安装featureCounts（通过subread包）
conda install -c bioconda subread -y

# 安装MultiQC用于质量控制报告
conda install -c bioconda multiqc -y

参考基因组准备

下载参考文件

下载网址： https://asia.ensembl.org/Mus_musculus/Info/Index

需要下载两个文件：

基因组序列文件：选择 dna.primary_assembly.fa.gz
基因注释文件：选择 GRCm39.114.gtf.gz

解压文件

# 解压基因组fasta文件 
gunzip Mus_musculus.GRCm39.dna.primary_assembly.fa.gz 

# 解压GTF注释文件 
gunzip Mus_musculus.GRCm39.114.gtf.gz

构建STAR索引

# 创建索引目录
mkdir -p star_index

# 构建STAR基因组索引
STAR --runThreadN 6 \
     --runMode genomeGenerate \
     --genomeDir star_index \
     --genomeFastaFiles Mus_musculus.GRCm39.dna.primary_assembly.fa \
     --sjdbGTFfile Mus_musculus.GRCm39.114.gtf \
     --sjdbOverhang 149

参数说明：

--runThreadN 6：使用6个线程
--runMode genomeGenerate：构建基因组索引模式
--genomeDir：索引输出目录
--genomeFastaFiles：基因组序列文件
--sjdbGTFfile：基因注释文件
--sjdbOverhang 149：测序读长减1（通常为150bp读长）

STAR比对详解

STAR（Spliced Transcripts Alignment to a Reference）是目前最流行的RNA-seq比对工具之一，特别擅长处理跨越剪接位点的读段。

核心比对命令

STAR --runThreadN ${THREADS} \
     --genomeDir ${INDEX_DIR} \
     --readFilesIn ${R1} ${R2} \
     --readFilesCommand zcat \
     --outFileNamePrefix ${OUT_DIR}/${sample}/ \
     --outSAMtype BAM SortedByCoordinate \
     --quantMode GeneCounts \
     --outFilterMultimapNmax 20 \
     --alignSJoverhangMin 8 \
     --alignSJDBoverhangMin 1 \
     --outFilterMismatchNmax 999 \
     --outFilterMismatchNoverReadLmax 0.04 \
     --alignIntronMin 20 \
     --alignIntronMax 1000000 \
     --alignMatesGapMax 1000000

关键参数解释

基础设置参数：

--runThreadN：并行线程数，根据服务器配置调整
--genomeDir：预构建的STAR索引目录
--readFilesIn：输入的FASTQ文件（R1和R2）
--readFilesCommand zcat：处理压缩文件的命令

输出格式参数：

--outSAMtype BAM SortedByCoordinate：直接输出按坐标排序的BAM文件
--quantMode GeneCounts：同时生成基因计数表
--outFileNamePrefix：输出文件名前缀

比对质量控制参数：

--outFilterMultimapNmax 20：每条读段最多允许20个比对位置
--outFilterMismatchNmax 999：最大错配数（999表示不限制）
--outFilterMismatchNoverReadLmax 0.04：错配率不超过4%

剪接位点参数：

--alignSJoverhangMin 8：跨越剪接位点的最小悬垂长度
--alignSJDBoverhangMin 1：注释剪接位点的最小悬垂长度

内含子和配对参数：

--alignIntronMin 20：最小内含子长度
--alignIntronMax 1000000：最大内含子长度
--alignMatesGapMax 1000000：配对读段间最大距离

featureCounts定量分析

featureCounts是Subread包中的核心工具，用于对比对结果进行基因水平的读段计数。

基本原理

featureCounts根据GTF注释文件，统计落在每个基因区域内的读段数量，是RNA-seq表达量分析的关键步骤。

核心命令

featureCounts -T ${THREADS} \
              -a ${GTF_FILE} \
              -o ${COUNT_DIR}/counts.txt \
              -g gene_id \
              -p \
              -B -C \
              --primary \
              $(find ${OUT_DIR} -name "Aligned.sortedByCoord.out.bam")

参数详解

基础参数：

-T：线程数
-a：GTF注释文件路径
-o：输出文件名
-g gene_id：使用基因ID作为特征标识符

配对末端参数：

-p：配对末端模式，将一对读段作为一个片段计数
-B：只计数两端都成功比对的读段对
-C：不计数嵌合比对的读段

比对质量参数：

--primary：只计数主要比对结果，忽略多重比对

输出结果处理

# 创建简化的基因计数矩阵（去除注释行，只保留基因ID和计数列）
grep -v '^#' ${COUNT_DIR}/counts.txt | cut -f1,7- > ${COUNT_DIR}/gene_counts_matrix.txt

完整分析流程脚本

将所有步骤整合成一个完整的自动化脚本：

#!/bin/bash
# RNA-Seq数据处理完整流程

# ========== 参数配置 ==========
INDEX_DIR="/home/kaka/newpan/kaka/project/base/gtf/mouse/star_index"
GTF_FILE="/home/kaka/newpan/kaka/project/base/gtf/mouse/Mus_musculus.GRCm39.114.gtf"
WORK_DIR="/home/kaka/newpan/kaka/project/25-6-30-bulk/work"
OUT_DIR="${WORK_DIR}/alignment"
COUNT_DIR="${WORK_DIR}/counts"
THREADS=6

# ========== 创建目录结构 ==========
mkdir -p ${OUT_DIR} ${COUNT_DIR} ${QC_DIR}

# ========== STAR比对分析 ==========
echo "开始STAR比对分析..."
for sample_dir in ${WORK_DIR}/fq/*; do
    sample=$(basename ${sample_dir})
    echo "正在处理样本: ${sample}"
    
    R1=${sample_dir}/${sample}.R1.clean.fastq.gz
    R2=${sample_dir}/${sample}.R2.clean.fastq.gz
    
    mkdir -p ${OUT_DIR}/${sample}
    
    STAR --runThreadN ${THREADS} \
         --genomeDir ${INDEX_DIR} \
         --readFilesIn ${R1} ${R2} \
         --readFilesCommand zcat \
         --outFileNamePrefix ${OUT_DIR}/${sample}/ \
         --outSAMtype BAM SortedByCoordinate \
         --quantMode GeneCounts \
         --outFilterMultimapNmax 20 \
         --alignSJoverhangMin 8 \
         --alignSJDBoverhangMin 1 \
         --outFilterMismatchNmax 999 \
         --outFilterMismatchNoverReadLmax 0.04 \
         --alignIntronMin 20 \
         --alignIntronMax 1000000 \
         --alignMatesGapMax 1000000
done

# ========== featureCounts基因定量 ==========
echo "开始featureCounts基因定量..."

# 重新创建计数目录
rm -rf ${COUNT_DIR}
mkdir -p ${COUNT_DIR}

# 执行基因计数
featureCounts -T ${THREADS} \
              -a ${GTF_FILE} \
              -o ${COUNT_DIR}/counts.txt \
              -g gene_id \
              -p \
              -B -C \
              --primary \
              $(find ${OUT_DIR} -name "Aligned.sortedByCoord.out.bam")

# 生成简化的基因计数矩阵
grep -v '^#' ${COUNT_DIR}/counts.txt | cut -f1,7- > ${COUNT_DIR}/gene_counts_matrix.txt

# ========== 分析完成 ==========
echo "=========================================="
echo "RNA-Seq数据分析流程完成！"
echo "=========================================="
echo "输出结果目录："
echo "比对结果: ${OUT_DIR}"
echo "计数结果: ${COUNT_DIR}"
echo "=========================================="

输出文件说明

比对结果文件

Aligned.sortedByCoord.out.bam：排序后的比对文件
Log.final.out：比对统计信息
ReadsPerGene.out.tab：STAR生成的基因计数表

计数结果文件

counts.txt：详细的featureCounts输出
gene_counts_matrix.txt：简化的基因表达矩阵
counts.txt.summary：计数统计摘要

结束语

通过这个完整的流程，你就可以从原始的FASTQ文件得到用于下游分析的基因表达矩阵了！

参考资料：

https://zhuanlan.zhihu.com/p/362727395

https://mp.weixin.qq.com/s/RblPPDhH2BCwux6sfF5jsQ

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

linux

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

linux

登录后参与评论

0 条评论

热度