基于体型RNA测序数据的转录组分析

Summary

本方案建立了一套完整的流程，用于分析从原始数据到功能富集分析的体质RNA测序过程。

Abstract

非酒精性脂肪肝（NAFL）通常被认为是良性疾病;然而，一旦发展为非酒精性脂肪性肝炎（NASH），患者罹患末期肝病的风险显著增加。许多研究正试图阐明从NAFL向NASH转变的分子机制。高通量测序技术（如体型RNA-seq）通过检查转录组，揭示了分子表达、信号通路激活及其他与疾病进展相关的因素，为研究人员提供了更深入的理解。研究人员有大量开源数据可供分析，以识别潜在的疾病治疗靶点。然而，相关研究受限于缺乏高效可靠的转录组上游分析流程。这里提供了高度可重复且用户友好的上游分析及后续的差异基因分析流程，以实现对私密或公共数据的标准化处理和深度解析。该流程分为四个步骤：（1）数据质量控制;（2）基因定位;（3）差异基因分析;以及（4）泛函分析。该过程旨在揭示疾病转化的分子机制，并通过分析Bulk RNA-seq数据，协助研究人员筛选潜在药物靶点和治疗方法。

Introduction

非酒精性脂肪肝（NAFLD）是全球最常见的慢性肝病，影响超过四分之一的人口。近几十年来，其发病率急剧上升，^1,2,3。日益增长的疾病负担，尤其是其更为严重的形式——非酒精性脂肪肝炎（NASH），构成了全球重大健康挑战和沉重的经济负担^4.NAFLD的第一阶段是非酒精性脂肪肝（NAFL），伴有炎症和纤维化，可能进展为NASH。后者显著增加肝病末期发展的风险，包括肝硬化和肝细胞癌（HCC）^5,6,7。HCC的发生率和死亡率与^NASH的增加相关^，预计到2030年，NAFLD/NASH将成为肝移植的主要适应^症。然而，NAFLD的临床进展^极为异质，严重阻碍了相关药物的开发，因此精确探究相关分子机制尤为重要。

基于RNA测序的体型细胞组成信息获取可以显著阐明多种疾病的发病机制。近几十年来，在模型生物和人类中进行了大量体质RNA-测序研究，旨在阐明NASH13^、14、15进展中的基因表达差异，以确定新的干预治疗靶点。基于整体RNA-seq分析，熊等发现肝脏中的非实质细胞（NPCs）参与细胞外基质形成和细胞粘附等过程，这些过程对NASH¹⁶的进展有贡献。Li等人证明肝细胞中的Wilms肿瘤1结合蛋白（WTAP）调控异位脂质积累和炎症，从而促进NASH形成¹⁷。尽管体质RNA测序分析是阐明NASH机制的有力工具，但其结果对上游数据质量极为敏感。上游实验作和分析过程的异质性会严重降低数据的可靠性，从而掩盖真实的生物学信息，干扰后续分析的准确性。因此，建立一套标准化的上游分析程序非常重要。

与单细胞RNA测序（scRNA-seq）相比，体型RNA-seq在实验设计和实际应用中具有多项显著优势。虽然scRNA-seq能够在单细胞层面识别细胞异质性，并精确分析细胞类型特异性转录特征，但其成本高昂、数据处理复杂，且检测低丰度转录本的灵敏度^{有限，因此}存在相关性。相比之下，大宗RNA测序提供了更高的测序深度、更低的成本和更高的样本通量，使其特别适合群体层面的差异基因表达分析和分子机制的^探索。因此，在标准化分析流程指导下，体型RNA测序依然是一种高效、经济且稳健的方法，用于研究复杂疾病的分子基础。

该方案专为来自高RNA完整性（RIN ≥ 7.0）且输入RNA充足（每样本≥500 ng）的人类组织的散体RNA测序数据集设计。为确保对齐和量化步骤的可靠执行，建议配备至少10核CPU、32GB内存和至少200GB可用磁盘空间的本地工作站。基于这些需求，该协议提供了高效且用户友好的分析流程，包括详细的作说明和标准化参数配置，以满足分析大规模转录组数据的研究人员需求。

Protocol

为示范目的，Lan Bai等人生成的公开数据集被用来说明上下游分析的每一步PRJNA1023502²⁰。由于该数据集来源于开放访问的NCBI SRA数据库，无需额外的权限或伦理审批。请参阅 材料表 以核实所有必需的软件和R包版本。公开数据集PRJNA1023502包括6个非NASH、6个NAFL和6个NASH肝RNA测序样本。该方案中，该数据集展示了体型RNA测序工作流程的所有步骤，包括从SRA数据库检索数据、质量控制（fastp）、比对（HISAT2）、定量（featureCounts）以及下游差异表达和功能富集分析。

1. SRA工具包安装

1. 请访问SRA工具包官方网站并下载3.2.1版本。

2. 公共数据下载

1. 获取SRA号码。
   1. 可通过文章补充²⁰、数据可用性部分访问，或在NCBI SRA数据库中搜索关键词访问。
   2. 在终端输入预取即可下载。

3. 生成基因计数矩阵

1. 请定义如下所述的路径和设置：
   REFERENCE=~/reference/human/GRCh38/GRCh38.primary_assembly.genome.fa
GTF=~/reference/human/GRCh38/gencode.v44.annotation.gtf
INDEX=~/reference/human/GRCh38/GRCh38_index
FASTQ_DIR=~/SRA_tutorial/fastq
OUT_FASTP=~/RNAseq/fastp
OUT_HISAT2=~/RNAseq/hisat2
OUT_COUNTS=~/RNAseq/counts
​mkdir -p $FASTQ_DIR $OUT_FASTP $OUT_HISAT2 $OUT_COUNTS
   1. 从官方GENCODE数据库（https://www.gencodegenes.org/human/）下载人类参考基因组（GRCh38初级组装）及对应的基因注释文件（gencode.v44，综合基因注释集）。
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
​for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
2. 重命名SRA文件。
   1. 确保所有SRA文件都使用.sra扩展名，以便下游工具准确识别和处理。
   2. 执行以下命令：
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
3. 执行以下命令将SRA文件转换为FASTQ格式：
   for f in *.sra; do fasterq-dump "$f" --split-files -O $FASTQ_DIR - e 20; done
4. 从参考基因组构建HISAT2指数。
   hisat2-build $REFERENCE $INDEX
5. 通过可视化进展监控处理FASTQ文件。
   for fq in $FASTQ_DIR/*.fastq; do
​sample=$(basename "$fq" .fastq)
   1. 对于配对端序列数据，执行以下命令
      for fq1 in $FASTQ_DIR/*_1.fastq; do
sample=$(basename "$fq1" _1.fastq)
​fq2=$FASTQ_DIR/${sample}_2.fastq
   2. 使用fastp进行质量控制和过滤，并将结果输出到$OUT_FASTP目录。对于单端测序数据，执行以下命令：
      fastp \
-i "${fq}" \
-o $OUT_FASTP/${sample}.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
-w 20
      对于配对端序列数据，执行以下命令：
      fastp \
-i "${fq}" \ -I "$fq2" \
-o $OUT_FASTP/${sample}_1.clean.fastq \
-O $OUT_FASTP/${sample}_2.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
​-w 20
   3. 对于每个样本，检查HTML报告以进行视觉质量控制，并参考JSON报告中的结构化指标，包括读质量、GC分布、读冗余和适配器污染。确保所有样品达到以下质量门槛：Q30≥85%，适配器污染率<5%，且无异常GC模式。将所有输出文件存储在 $OUT_FASTP 目录中。
   4. 使用HISAT2进行序列比对，生成SAM文件，并将其输出到$OUT_HISAT2目录。对于单端测序数据，执行以下命令：
      hisat2 -p 20 \ -x $INDEX \-U $OUT_FASTP/${sample}.clean.fastq \
-S $OUT_HISAT2/${sample}.sam
      对于配对端序列数据，执行以下命令：
      hisat2 -p 20 \-x $INDEX \-1 $OUT_FASTP/${sample}_1.clean.fastq \
-2 $OUT_FASTP/${sample}_2.clean.fastq \
​-S $OUT_HISAT2/${sample}.sam
   5. 使用由GRCh38主要组装FASTA文件构建的索引，对经过质量过滤的读段进行HISAT2比对至人类参考基因组。提供单端和配对端测序数据的示例命令。
   6. 将SAM转换为BAM，排序并索引。
      samtools view -@ 20 -bS $OUT_HISAT2/${sample}.sam \
| samtools sort -@ 20 -o $OUT_HISAT2/${sample}.sorted.bam
samtools index $OUT_HISAT2/${sample}.sorted.bam
​done
   7. 为每个采样生成一个SAM文件，并用SAMtools将其转换为排序和索引的BAM文件。对于高质量的人类转录组样本（RNA完整性数，RIN ≥ 8.0），确保整体比对率超过85%;对于典型的散体RNA测序样本（RIN ≥ 7.0），比对率≥可接受70%。
   8. 使用featureCount进行基因定量。
      featureCounts -T 20 -p -s 0 \
-a $GTF \
-o $OUT_COUNTS /${sample}.counts.txt \
$OUT_HISAT2/${sample}.sorted.bam
Done
   9. 检查由 featureCounts 生成的制表符分隔输出文件（*.counts.txt）及其摘要报告（*.counts.txt.summary）。确保读段分配率达到人类整体RNA测序的典型阈值≥70%;明显较低的比率可能表示链方向不匹配、注释问题或比对质量差。对于这个非链特异性RNA测序数据集，使用参数-s 0。对于特定链库，命令中用 -s 1 或 -s 2 替代 -s 0。
6. 生成如下所述的基因计数矩阵。
   1. 用基因ID和第一个样本的计数初始化计数矩阵。
      cut -f1 $(ls $OUT_COUNTS/*.counts.txt | head -1) > all_counts.txt
   2. 将每个样本的计数反复加到单一矩阵中。
      for f in $OUT_COUNTS/*.counts.txt; do
cut -f7 "$f" | paste all_counts.txt - > tmp && mv tmp
all_counts.txt
​done
   3. 在计数矩阵中添加带有样本ID的头部行。
      samples=$(ls *.counts.txt | sed 's/.counts.txt//' | paste -sd "\t")
echo -e "Geneid\t$samples" | cat - all_counts.txt > counts_matrix.txt
   4. 从GTF文件中提取基因长度（每个基因外显子长度的总和）。
      awk '$3=="exon"{match($0,/gene_id "([^"]+)"/,a); if(a[1]!=""){len=$5-$4+1; gene_len[a[1]]+=len}} END{print "GENE_ID\tLENGTH"; for(g in gene_len) print g"\t"gene_len[g]}' \$GTF > gene_length.txt

4. 原始计数矩阵处理与基因注释

1. R环境初始化和加载相关R包。
   1. 使用 install.packages（） 函数安装 tidyverse、ggplot2、ggrepel、RColorBrewer、ggridges、FactoMineR。tidyverse 包用于数据作和绘图;GGPLOT2 软件包用于可视化;Ggrepel包用于不重叠的文本标签;色彩调色板采用了RColorBrewer包;Ggridges套餐用于山脊地块;FactoMineR 软件包用于主成分分析（PCA）和多变量分析。
   2. 使用 BiocManager：：install（） 函数来安装 biomaRt、DESeq2、clusterProfiler。biomaRt 软件包用于 Ensembl 的基因注释;DESeq2 软件包用于差分表达式分析;clusterProfiler 包用于功能丰富分析。
   3. 读取步骤2生成的原始计数矩阵，使用biomaRt将Ensembl ID映射到HGNC基因名进行后续分析，去除重复的Geneid列，生成清理后的计数矩阵。输入原始计数矩阵（counts_matrix.csv），行用Ensembl基因ID，列用样本。将清理后的计数矩阵保存为输出（clean_counts_SRA.csv），HGNC基因符号为行，样本为列。
      mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
id_map <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = exprSet$GeneID,
mart = mart)
exprSet <- exprSet %>%
left_join(id_map, by = c("GeneID" = "ensembl_gene_id")) %>%
filter(!is.na(hgnc_symbol), hgnc_symbol != "") %>%
distinct(hgnc_symbol, .keep_all = TRUE) %>%
column_to_rownames("hgnc_symbol")

5. 基因表达定量

注意：详细脚本请参阅 补充文件1 。

1. 执行以下命令，计算每千碱基转录本每百万映射读取（FPKM）和每百万转录本（TPM）片段数。
   counts <- read.csv("output/clean_counts_SRA.csv", header=TRUE, row.names=1)
gene_len <- read.delim("data/gene_length.txt", header=FALSE, col.names=c("gene_symbol","length"))
gene_len <- gene_len %>% distinct(gene_symbol, .keep_all=TRUE)
rownames(gene_len) <- gene_len$gene_symbol
gene_len <- gene_len[match(rownames(counts), gene_len$gene_symbol),]
length_bp <- gene_len$length
fpkm <- (counts / length_bp) * 1e9 / colSums(counts)
write.csv(fpkm, "output/clean_fpkm_SRA.csv")
tpm <- (counts / length_bp) / colSums(counts / length_bp) * 1e6
write.csv(tpm, "output/clean_tpm_SRA.csv")

6. 样本聚类与差分可视化

1. 通过PCA检查样本聚类来评估数据质量。详细脚本请参见 补充文件1 。
2. 对于PCA可视化，请使用下面的代码。
   gene.pca <- PCA(exprSet, ncp = 2, scale.unit = TRUE, graph = FALSE)
ggplot(pca_sample, aes(x = Dim.1, y = Dim.2)) +
geom_point(aes(color = group)) +
labs(x = paste('PC1:', pca_eig1, '%'),
y = paste('PC2:', pca_eig2, '%'))

7. 差异表达分析与结果可视化

注意：详细脚本请参阅 补充文件1 。

1. 执行以下命令构建DESeq2数据集并设定阈值。建模前仅保留了所有样本总计数>1的基因。利用调整后的p值0.05和|log2FC| 的显著阈值<识别差异表达基因（DEGs）。> 0.5。
   dds <- DESeq(DESeqDataSetFromMatrix(countData = exprSet, colData = colData, design = ~group)); sizeFactors(dds); res <- results(dds); dds <- dds[rowSums(counts(dds)) > 1,]
dd1 <- results(dds, contrast = contrast, alpha = 0.05)
dd2 <- lfcShrink(dds, contrast = contrast, res = dd1, type = "ashr")
2. 用火山图可视化差异表达的基因。
   ggplot(data = data, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(aes(color = group), alpha = 1, size = 1.2) +
geom_hline(yintercept = -log10(0.05), lty = 4) +
geom_vline(xintercept = c(-0.5, 0.5), lty = 4) +
geom_text_repel(data = subset(data, abs(log2FoldChange) >= 1.5 & padj < 0.05),
aes(label = gene_id))

8. 进行功能丰富分析与可视化

注意：详细脚本请参阅 补充文件1 。

1. 可视化KEGG浓缩结果。
   EGG <- enrichKEGG(gene = gene$ENTREZID, organism = 'hsa',
pvalueCutoff = 0.05, qvalueCutoff = 0.05)
ggplot(symboldata, aes(richFactor, Description)) +
geom_point(aes(color = p.adjust, size = Count))
2. 可视化 GO 的丰富效果。通过根据实验规格修改弧分辨率（n=500）、色彩映射和面排列，定制可视化。
   ego <- enrichGO(gene = gene$ENTREZID, OrgDb = "org.Hs.eg.db", ont = "ALL",
pvalueCutoff = 0.05, qvalueCutoff = 0.05, pAdjustMethod = "BH")
ggplot(df) +
ggforce::geom_link(aes(x = 0, y = Description, xend = -log10(p.adjust),
yend = Description, color = ONTOLOGY), n = 500, show.legend = FALSE) +
facet_wrap(~ONTOLOGY, scales = "free", ncol = 1)
3. 可视化基因集富集分析（GSEA）结果。调整通路数量、密度脊尺度和层透明度，以适应不同的效应大小分布。
   genelist <- sort(res$log2FoldChange, decreasing = TRUE)
names(genelist) <- rownames(res)
hallmarks <- read.gmt('resource/h.all.v2023.2.Hs.symbols.gmt')
y <- GSEA(genelist, TERM2GENE = hallmarks, pvalueCutoff = 0.05)
gsearesult <- yd %>% arrange(desc(NES)) %>% slice_head(n = 10)
ggplot(gsearesult, aes(x = logFC, y = Description, fill = -log10(pvalue))) +
geom_density_ridges(alpha = 0.8, scale = 0.8) +
geom_point(aes(size = abs(NES), x = -0.4, color = NES)) +
scale_fill_distiller(palette = 'Spectral') +
scale_color_distiller(palette = 'Reds') +
scale_size_continuous(range = c(2, 6))

Results

体型RNA测序的上游分析流程如 图1A所示。该工作流程在 Linux 平台上依次执行以下关键步骤：首先，使用 fastp 对原始测序数据进行严格质量控制，以去除低质量读段和适配器序列;随后，HISAT2 将高质量读段比对到参考基因组，Samtools 对比对文件进行转换和排序;最后，FeatureCounts执行基因级定量，生成基因表达矩阵，为后续分析提供高质量输入。后续对生成表达式矩阵的处理和统计分析在R环境中进行，相关工作流程和所需软件包如 图1B所示。分析中选取了一项已发表研究的数据，包括6个非NASH、6个NAFL和6个NASH样本²⁰ （见图1C）。非NASH对照样本来自不符合肝移植标准且无NAFLD或NASH的个体。

需要注意的是，由于Bai等人公开的数据未提供明确的批次信息（如批次测序或文库制备日期），也未提供显著表达基因的可下载列表，本研究无...

Discussion

体质RNA测序数据分析被定义为一项跨学科任务，融合了基因组学、生物信息学、统计学和计算机科学。完整的分析流程涵盖多个上下游步骤，包括原始数据预处理、质量控制、序列比对、基因级定量、数据归一化、差异表达分析和生物解读。在这些步骤中，准确将原始测序读段转换为高质量的基因表达矩阵尤为关键，因为上游处理过程中引入的错误可能传递到所有下游生物学结论中。因此，建立透明且标准化的上游分析工作流程对于提高转录组研究的可重复性至关重要。

该协议提供了一个简化的、完全基于脚本的工作流程，集成了广泛使用的工具，如fastp（用于读取裁剪和质量控制）、HISAT2（用于剪接感知比对）和featureCounts（用于基因级定量）。这些工具已被广泛应用于成熟的RNA-seq分析框架中——包括基于Tuxedo的流程和常用的协议工作流程——并经过了准确性和效率的严格验证（^21,22）。在这些基础方法的基础上，工作流程通过在 Li...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
生物体	生物导体	2.64.0	来自Ensembl的基因注释
clusterProfiler	生物导体	4.16.0	功能富集分析
DESeq2	生物导体	1.48.1	差分表达分析
事实矿山R	农业巴黎科技	2.11.0	主主成分分析与多变量分析
FASTP	开放基因	1.0.1	FASTQ数据的质量控制与过滤
功能计数	沃尔特与伊丽莎白·霍尔医学研究所生物信息学部	2.0.0	统计每个基因映射的读段数量以定量基因表达
GGPLOT2	假设	3.5.2	数据可视化
格雷佩尔	卡米尔·斯沃维科夫斯基	0.9.6	不重叠的文本标签
格里奇斯	克劳斯·O·威尔克	0.5.6	创建山脊线地块
HISAT2	约翰斯·霍普金斯大学	2.2.1	将过滤后的高质量读段与参考基因组比对
R	R核心团队	4.5.0	一个用于数据计算、分析和可视化的环境
RColor酿酒师	埃里希·诺伊维尔特	1.1.3	用于绘图的色彩调色板
samtools	大规模基因组学工作流	1.22.0	转换和处理SAM文件，以实现高效的检索和访问
SRA工具包	国家生物技术信息中心	3.2.1	从NCBI SRA数据库获取并预处理原始测序数据