8.5KB
转录组数据分析报告
分析人：
复旦大学生命科学学院PGx group石乐明老师课题组
1. 简介
项目目的
实验设计
分析内容
采用FastQC、FastQ Screen、QualiMap对数据进行质控分析，并使用MultiQC进行可视化；
使用 Trimmomatic + HISAT2+ StringTie + ballgown 构建的分析流程对转录组数据进行分析获得表达谱；
根据获得的表达谱进行表达分析、差异基因分析，并进行功能分析（GO/KEGG）;
主要结论
- 基于Ensembl注释，我们获得了个样本个基因的表达图谱；
- 主成分分析PCA和分层聚类HCA分析发现??聚类情况？
- 差异基因分析，分析方法为 P 值小于 0.05， 倍数变化大于 2 或小于 0.5，共获得个差异表达基因，进行功能富集后发现一系列与表达相关的富集通路
- 整体数据质量情况：。原始数据中存在少量adapter残留，通过序列修剪后可以进行后续分析；
- 样品特征回溯是数据质量控制的方法之一。
2. 数据质量控制
2.1 原始数据质量
通过FastQC软件获得的结果可知，本次转录组测序的原始数据量范围在（22~32）M reads之间，平均值为27 M reads；其碱基质量值在35之间，表明测序质量较好；GC含量、碱基含量、数据重复率、重复序列出现比例等指标均未见异常；通过FastQC检测，发现所得的原始数据内包含有较多的接头未去除，因此使用了Trimmomatic软件先对原始文件进行去接头操作，之后数据分析的结果都使用的为去接头之后的数据（clean data）。
（注：FastQC软件结果说明参照 fatqc introduction）
@data-table-js(dataUrl=‘./data/multiqc_general_stats.csv’)
FastQC软件所获得的结果为每个样品所得到的结果，我们使用MultiQC将本次测序所获得的所有样本数据质控情况集合在一个报告文件中进行可视化展示。
（注：MultiQC软件结果说明参照：https://multiqc.info/docs/#fastqc）
@multiqc(analysisDir=‘./data/fastqc’)
FastQ Screen允许针对多个物种进行比对，以确定测序文库是否与你预期的一样，可用于检查是否存在序列污染。同样我们使用multiqc进行可视化汇总。
通过所获得的结果可知，序列比对到的基因组为人类的基因组序列。
（注：FastQ Screen软件说明FastQ Screen Documentation）
@multiqc(analysisDir=‘./data/fastqscreen’)
2.2 比对结果质量
在对原始数据进行比对之后，我们获得了比对结果（即bam文件），可以通过使用QualiMap软件对比对后的结果进行质控，主要关注数据的 mapping 上参考基因组的比例。
通过比对质控的结果可知， 比对到人基因组读段比例在。
（注：QualiMap软件结果说明：Qualimap）
@multiqc(analysisDir=‘./data/qualimap’)
2.3 样品特征回溯
根据基因表达谱特征，我们对样品的物种、性别和器官来源进行了回溯，结果如下：
@data-table-js(dataUrl=‘./data/sample_predict.csv’)
3. 表达谱结果
3.1 表达谱结果描述
基因表达量的计算使用FPKM (Fragments Per Kilobase of transcript per Million mapped reads)，FPKM能消除基因长度和测序量差异对计算基因表达的影响。我们对各个样本的基因表达量以 2 为底求对数值。
3.2 整体聚类情况
在RNA组学研究中利用主成分分析（PCA），将样本所包含的上万个维度的信息（上万个基因的表达量），降维为数个维度的综合指标（主成分），以便于进行样本间的比较，同时保证原始数据中包含信息尽可能多地被保留 。
@scatter-plot(dataFile=‘./data/rnaseq_pca.rds’, dataType=‘rds’, xAxis=‘PC1’, xTitle=“PC1”,yAxis=‘PC2’,yTitle=“PC2”,colorAttr=“group1”)
我们基于基因表达量，对样本和基因间的关系进行层级聚类，并使用热图来呈现聚类结果。对不同样品和基因进行层级聚类分析，图中每列代表一个样品，每行代表一个基因，基因在不同样品中的表达量用不同颜色表示。
@heatmap-d3(dataFile=‘./data/rnaseq_cor.rds’, dataType=‘rds’, labCol=‘TRUE’)
3.3 差异基因
差异表达基因，我们采用统计方法 t 检验对本次分析样本的表达谱进行了比较。差异基因选定的阈值是 P 值小于 0.05， 倍数变化大于 2 或小于 0.5，过滤低表达基因（所有基因表达量均低于FPKM 0.1）后，发现差异基因1064个，其中组P1-6（A）较组P7-13（B）高的基因504个，A较B显著低的基因560个，数量如下图所示：
@stack-barplot-r(dataFile=‘./data/rnaseq_degs_stats.rds’, dataType=‘rds’, xAxis=‘versus’, yAxis=‘number’,labelAttr=‘type’,barPos=‘stack’)
两组间的差异表达基因如下所示：
@scatter-plot(dataFile=‘./data/rnaseq_AvsB_degs.rds’, dataType=‘rds’, xAxis=‘logfc’, xTitle=“log2FC”,yAxis=‘log10p’,yTitle=“-log10 (p)”, colorAttr=“sigene”)
差异基因清单如下：
@data-table-js(dataUrl=‘./data/rnaseq_degs_acrossgroups.csv’)
3.4 基因功能富集分析
得到差异表达基因之后，我们对差异表达基因做GO功能分析和KEGG Pathway分析。
GO功能分析
Gene Ontology（简称GO）是一个国际标准化的基因功能分类体系，提供了一套动态更新的标准词汇表（controlled vocabulary）来全面描述生物体中基因和基因产物的属性。GO总共有三个ontology（本体），分别描述基因的分子功能（molecular function）、细胞组分（cellular component）、参与的生物过程（biological process）。GO的基本单位是term（词条、节点），每个term都对应一个属性。 GO功能分析一方面给出差异表达基因的GO功能分类注释；另一方面给出差异表达基因的GO功能显著性富集分析。
KEGG Pathway分析
在生物体内，不同基因相互协调行使其生物学，基于Pathway的分析有助于更进一步了解基因的生物学功能。KEGG是有关Pathway的主要公共数据库。Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著性富集的Pathway。
(1) 采用GSEA分析基因通路分析
利用fgsea包根据基因的表达水平进行基因功能分析，通路如下所示：
富集的通路结果：
@data-table-js(dataUrl=‘./data/rnaseq_gsea_curatedgenesets.csv’)
富集的GO功能：
@data-table-js(dataUrl=‘./data/rnaseq_gsea_go.csv’)
(2) 采用clusterprofiler分析差异基因富集的KEGG和GO通路分析
采用clusterprofiler分析差异基因富集的KEGG和GO通路分析，结果如下所示：
KEGG通路：
@data-table-js(dataUrl=‘./data/rnaseq_KEGGenrich.csv’)
GO功能：
@data-table-js(dataUrl=‘./data/rnaseq_GOenrich.csv’)
4. 附录：数据分析流程介绍
利用 Trimmomatic软件对测序数据进行去接头引物和修剪低质量序列，利用 HISAT2将高质量序列比对到人的参考基因组hg38上，基于Ensembl基因模型，利用StringTie进行转录本重构和定量，利用Ballgown对结果进行基因水平转化。此外，我们采用 FastQC、FastQ_Screen和MultiQC等软件对测序数据进行质量评估，利用 QualiMap对比对数据进行质量评估。主成分分析、相关性分析、聚类分析、差异基因和功能富集分析利用R/Bioconducter进行分析。