微科盟DIA无参蛋白质组结题报告

一、概述

  Astral DIA是基于新一代Orbitrap Astral高分辨质谱仪开发的蛋白质组学。2023年6月4日，赛默飞在美国质谱年会（ASMS）上全球首发新一代Orbitrap Astral高分辨质谱仪。Orbitrap Astral高分辨质谱实现了DIA采集窗口2Th，与常规DDA蛋白质组采集窗口相同，但不同的是2Th窗口采集的DIA无论是鉴定数目还是定量CV<20%的占比都优于DDA。它结合了四极杆质量分析器（Quadrupole）、轨道阱质量分析器（Orbitrap）以及全新非对称轨道无损质量分析器（Astral），具有更高通量、更深覆盖、更高灵敏度的特点，可实现大规模蛋白质组学研究。下图为Astral DIA定量技术原理图：

图 1.1.1 Astral DIA定量技术原理图
图 1.1.1 Astral DIA定量技术原理图

  Astral DIA蛋白质组学是一种LFQ数据非依赖分析（DIA）的多肽定量方法，它基于质谱分析技术来实现定量测定。该技术通过直接比较样品中多肽的质谱特征，来推断其在不同样品中的相对含量。

二、实验流程

  从组织样品到最终数据获得的过程中，蛋白的提取、定量、检测、酶切与除盐1、馏分分离和质谱检测2每一个环节都会对数据质量和数量产生影响，而数据质量又会直接影响后续信息分析的结果。为了从源头上保证测序数据的准确性、可靠性，微科盟对每一个实验步骤都严格把控，从根本上确保了高质量数据的产出，流程图如下：

图 2.1.1 Astral DIA实验流程图
图 2.1.1 Astral DIA实验流程图

  在Astral DIA蛋白质定量技术中，样品中的多肽首先经过消化，产生肽段。然后，这些肽段通过液相色谱（LC）进行分离，使其逐一进入质谱仪。在质谱仪中，肽段被离子化形成离子，并通过质谱分析生成质谱峰。 质谱峰是质谱仪输出的一系列信号，用于表示肽段的特征。通过分析质谱峰的位置、强度和面积等信息，可以推断出样品中多肽的相对含量。这种定量方法无需引入标记剂，避免了标记剂对样品的干扰和影响，同时具备较高的准确性和灵敏度。

三、分析流程

  Astral DIA无参蛋白质组分析云流程是微科盟针对缺少公共数据库注释的物种，或自建库推出的基于蛋白质序列功能注释的蛋白质组学分析云流程。该流程可以使用自建库，或其他自定义蛋白质数据（如多物种合并蛋白质库）作为数据源。基于质谱检测得到的Raw文件，进行对应数据库的搜索，然后基于数据库搜索的结果进行蛋白质鉴定，同时进行肽段、蛋白和母离子质量容差分布分析来评定质谱检测数据的质量；使用eggNOG-mapper对鉴定到的蛋白进行常见功能数据库注释，包括GO、KEGG、COG等注释信息；接下来进行蛋白质的定量分析，包括鉴定到的蛋白质表达量基础统计、差异蛋白的筛选、表达模式聚类分析；最后针对筛选出来的差异蛋白进行GO、KEGG功能富集分析和互作网络分析等一系列的差异蛋白功能分析。通过WGCNA分析描述不同样品之间蛋白质关联模式，分析高度协同变化基因集和表型之间的关系。对于多样本数据（单个分组30个以上样本），补充了ANOVA单因素方差分析用于寻找组间特异性表达的蛋白。Astral DIA无参蛋白质组云流程解决了蛋白质组学分析中：自建库，物种合并库等分析困难的问题，以及KEGG，Stringdb等数据库部分物种未被收录导致下游分析无法进行的问题。

图 3.1.1 鉴定蛋白分析流程图
图 3.1.1 鉴定蛋白分析流程图

四、鉴定蛋白基本分析

点击打开文件夹

点击查看文件夹结构 
## ./0.Sample_description
## ├── ecdf-intensities.png
## ├── log-intensities.png
## ├── pearson.png
## ├── sum-intensites.png
## └── 蛋白基础分析readme.png

4.1 鉴定总蛋白数柱状图

  采用柱状图来展示各样本中鉴定到的蛋白总定量值。

图 4.1.1 鉴定总蛋白柱状图
图 4.1.1 鉴定总蛋白柱状图

注：X轴表示蛋白定量值的和，Y轴表示不同样本。

4.2 鉴定蛋白数定量值累积分布图

  对蛋白定量值从高到低进行排序后做出累积分布图（图 4.2.1），可以用来展示每种蛋白质定量值对总蛋白质的累积贡献，一般按照样本组进行分组画图分析。

图 4.2.1 蛋白定量值累积分布图
图 4.2.1 蛋白定量值累积分布图

注：X轴表示蛋白定量值从高到低排序后蛋白对应坐标位置，Y轴表示 蛋白定量值从高到低排序后累积占比，可观察出定量值高的蛋白对整体蛋白量的贡献度，图中曲线上升得越快表示蛋白质定量值对总蛋白累计贡献度越高；颜色表示不同样本。

4.3 鉴定蛋白数 Rank 图

  对蛋白定量值从高到低进行排序后画图。可以用来展示鉴定总蛋白数、定量值数量级、以及不同样本之间定量值密度分布程度。

图 4.3.1 Protein Rank图
图 4.3.1 Protein Rank图

注：X轴表示蛋白定量值从高到低排序后蛋白对应坐标位置（即按定量值排序后第n个蛋白），Y轴表示蛋白定量值进行log10计算，本项目实现了近5个数量级定量值跨度的设别；颜色表示不同样本。

五、蛋白定量分析

点击打开文件夹

点击查看文件夹结构 
## ./1.Basic_statistics
## ├── boxplot
## │   ├── Category1_boxplot.svg
## │   ├── Category1_data_matrix.xls
## │   └── acc_symbol.txt
## ├── heatmap
## │   ├── Category1_correlation_heatmap.svg
## │   ├── Category1_correlation_matrix.xls
## │   └── Category1_p_value_matrix.xls
## ├── pca
## │   ├── Category1_PCA.svg
## │   ├── Category1_PCA_3D.svg
## │   └── Category1_pca_points_ordinates.xls
## ├── violin
## │   ├── Category1_data_matrix.xls
## │   └── Category1_violin.svg
## └── 蛋白定量分析readme.png

5.1 数据标准化

  为展示各个样本的蛋白表达量的整体的分布情况，我们对所有样本的蛋白丰度值，用R4.3语言程序作“箱线图”和“小提琴图”。为避免图中数据点的跨度太大，我们对蛋白质定量值，做了log2(定量值+1)的缩小变换处理。

图 5.1.1 蛋白原始定量值分布
图 5.1.1 蛋白原始定量值分布

注：X轴表示样本名，Y轴表示蛋白定量值进行log2计算，本项目实现了近5个数量级定量值跨度的设别；颜色表示不同样本；距离箱图较远的点表示离群值。

图 5.1.2 鉴定总蛋白小提琴图
图 5.1.2 鉴定总蛋白小提琴图

注：小提图琴是个密度图，越宽的地方说明对应于纵轴表达量的基因越多。

5.2 PCA 分析

  主成分分析(PCA)也常用来评估组间差异及组内样本重复情况，PCA采用线性代数的计算方法，对数以万计的基因变量进行降维及主成分提取。我们对所有样本的蛋白质定量值进行PCA分析，如下图所示。理想条件下，PCA图中，组间样本应该分散，组内样本应该聚在一起。

图 5.2.1 样品分组PCA分析二维平面图
图 5.2.1 样品分组PCA分析二维平面图

注：图中横坐标为第一主成分，纵坐标为第二主成分，不同的颜色表示不同的分组。样本点之间的距离近似于样本之间各蛋白丰度差异的总和。百分数：表示成分的贡献率。

图 5.2.2 样品分组PCA分析三维立体图
图 5.2.2 样品分组PCA分析三维立体图

注：图中X坐标为第一主成分，Y坐标为第二主成分，Z坐标为第三主成分，不同的颜色表示不同的分组。样本点之间的距离近似于样本之间各蛋白丰度差异的总和。百分数：表示成分的贡献率

5.3 Pearson 相关性分析

  为了反映样本间蛋白表达的相关性，本研究计算了每两个样品之间所有定量蛋白表达量的Pearson相关系数，并将这些系数以热图的形式反映出来，如下图 5.3.1 所示。右上三角为样本之间Pearson相关性值颜色表示，值参考右边条带图。左下三角为相关性系数值，值越高说明相似程度越高（最大值为1）。

图 5.3.1 Pearson 相关性
图 5.3.1 Pearson 相关性

注：X、Y轴均代表每个样品。颜色代表相关性系数（Pearson’s correlation coefficients），颜色越蓝代表正相关性越高，颜色越红代表负相关性越高；圆圈的大小表示变量重要性（variable importance）

5.4 层次聚类

  利用Pearson相关系数反映样本间蛋白表达量的相关性，并采用欧式距离进行层次聚类分析（图 5.4.1），该图能够直接反映每个样本之间的关系。

图 5.4.1 样本层次聚类
图 5.4.1 样本层次聚类

注：热图中的单元格颜色越深，表示对应的两个样本之间的相关性系数越高。同一个小组内的样本应该具有较高的相关性系数。

六、差异表达蛋白分析

点击打开文件夹

点击查看文件夹结构 
## ./2.Differential_expression
## ├── Category1_Model.vs.Contrl
## │   ├── Category1_Model.vs.Contrl_All.xls
## │   ├── Category1_Model.vs.Contrl_All_symbol_id.xls
## │   ├── Category1_Model.vs.Contrl_DEG.xls
## │   ├── Category1_Model.vs.Contrl_DEG_symbol_id.xls
## │   ├── Category1_Model.vs.Contrl_Volcano.svg
## │   ├── Category1_Model.vs.Contrl_deeseq-density.png
## │   ├── heatmap.svg
## │   ├── input.xls
## │   └── table.xls
## ├── Vennupset
## │   └── venn_only1_group
## └── 差异蛋白分析readme.png

6.1 差异蛋白筛选（group1）

  蛋白质差异性分析能够极大地推动发现新的生物标志物，提升生物标志物鉴定的精准度，对分子作用机理、生物标志物、疾病早诊、分子分型、预后以及临床诊断等均具有重要价值。蛋白差异表达的输入数据为蛋白表达丰度表数据。分析时我们采用基于负二项分布的DESeq2进行分析3，具体使用基于BiocManager，getopt，ggplot2，DESeq2等R包的R语言脚本；对各个样本分组进行两两之间的比较，找到在不同分组中表达差异的蛋白。 差异蛋白的筛选标准是非常重要的，我们给出的标准|log 2 (FoldChange)| > 1 & p值 < 0.05是常用的经验值，也是我们分析流程的默认值。在实际项目中可以根据情况灵活选择，例如，差异倍数可以选择1.5倍，也可以选择3倍，p值常用的阈值包括0.01、0.05、0.1等。若按照以上标准筛选得到的差异蛋白过少，很有可能导致后面的功能富集分析没有显著性结果。若项目实验只关注某几个基因的表达情况(如基因敲除)， 不在意富集结果，从下面的差异分析表格中筛选关注的那几个基因即可。反之， 如果得到的差异蛋白数目过多，不利于后续目标基因的筛选，这个时候可使用更严格的阈值标准进行筛选。生科云的云流程提供了调整筛选阈值的途径。蛋白表达差异分析结果见下表 6.1。差异蛋白火山图见下图6.1.1，差异倍数密度图见下图6.1.2。

注：
    1. Features: 蛋白的 Accession ID
    2. baseMean: 样本标准化counts的均值
    3. log2FoldChange: 处理组与对照组蛋白质表达水平的比值，再经过差异分析软件收缩模型处理，最后以2为底取对数
    4. IfcSE: 标准误，标准误越小，样本均值和总体均值差距越小；标准误越大，样本均值和总体均值差距越大
    5. stat: wald检验的统计量，它实际上是logFC除以lfcse
    6. pvalue: p值
    7. padj: BH方法校正后的p值

图 6.1.1 差异蛋白火山图
图 6.1.1 差异蛋白火山图
注：
    1. 横坐标表示两个分组之间的差异倍数，其绝对值越大说明某基因在两组之间的表达差异也越大。该值为正时，表示差异上调；该值为负时，表示差异下调
    2. 纵坐标代表蛋白表达量变化的统计学显著程度
    3. 图中的散点代表各个蛋白，灰色圆点表示无显著性差异的蛋白，红色圆点表示显著上调的差异蛋白，绿色圆点表示显著下调的差异蛋白

图 6.1.2 差异倍数密度图
图 6.1.2 差异倍数密度图

注： X轴为log2(FC)，Y轴表示蛋白在组间的倍数密度，即该差异倍数下的蛋白数与总数的比例。理论上绝大部分蛋白是不显著差异，所以FC峰值位置应位于0附近，并呈现正态分布。

图 6.1.3 差异蛋白聚类热图
图 6.1.3 差异蛋白聚类热图

注： 图中用颜色代表表达量，越红表示蛋白质表达量越高，越蓝表示表达量越低，数据跨度太大会导致差异无法区分，我们对表达量进行标准化处理（减去平均值，除以标准差）；聚类到一起的样本表达模式比较接近。

6.2 差异蛋白韦恩图(venn-group1)

  差异蛋白韦恩图是对差异分析结果进行进一步分析展示，summary_up_down.svg(图6.2.1)汇总了各差异蛋白集上调(处理组表达量高于对照组)，下调数目(处理组表达量低于对照组)。如果有多个比较方案，多个差异蛋白集，我们将用Venn图和UpSet图（DEG_venn.svg|图6.2.2，DEG_upset.svg|图6.2.3）展示各个差异蛋白集之间的交集大小(共有的差异蛋白数目)，以及各个差异蛋白集特有的差异蛋白。Venn图和UpSet图表达的信息是一样的，只是用了不同的数据可视化形式。Venn图最多能展示五个差异蛋白集的关系，而UpSet图能展示的差异蛋白集个数能够更多。 该部分计算并绘制了所选分组所有差异分析结果两两组合的结果，如需3个及以上差异结果的venn图，可在分步骤中操作获取。

图 6.2.1 各差异蛋白集上调，下调数目汇总柱形图
图 6.2.1 各差异蛋白集上调，下调数目汇总柱形图

注： 图中红色表示表达量上调的蛋白(处理组表达量高于对照组)，绿色表示表达量下调的蛋白。

点击查看差异蛋白韦恩图和Upset图（当输入包含多个Category时，才输出该部分的结果）
图 6.2.2 差异蛋白韦恩图
图 6.2.2 差异蛋白韦恩图

注： 图中每个圈代表一个组间的差异蛋白，重叠部分代表多个组间共同的差异蛋白，非重叠部分代表相应组间特异性的差异蛋白。

图 6.2.3 各差异蛋白集之间的UpSet图
图 6.2.3 各差异蛋白集之间的UpSet图
注： 图中柱形图表示差异蛋白的数量，下面的黑点，一个表示该差异蛋白集特有的蛋白（补集中该组特有的蛋白），两个点表示两个差异蛋白集共有的蛋白数（交集）。

七、差异蛋白功能富集分析

  根据对鉴定到的所有差异蛋白，我们进行了GO、KEGG、结构域等富集分析，目的是检测差异表达蛋白是否在某些功能类型上有显著性的富集趋势。

7.1 差异表达蛋白 GO 富集分析

7.1.1 GO 富集分析结果

点击打开文件夹

点击查看文件夹结构 
## ./3.go_enrichment
## ├── Category1_Model.vs.Contrl
## │   ├── Category1_Model.vs.Contrl_ALL_barplot.svg
## │   ├── Category1_Model.vs.Contrl_ALL_cnetplot.svg
## │   ├── Category1_Model.vs.Contrl_ALL_dotplot.svg
## │   ├── Category1_Model.vs.Contrl_BP_dagplot.png
## │   ├── Category1_Model.vs.Contrl_BP_idagplot.svg
## │   ├── Category1_Model.vs.Contrl_CC_dagplot.png
## │   ├── Category1_Model.vs.Contrl_CC_idagplot.svg
## │   ├── Category1_Model.vs.Contrl_GO_Bubble.svg
## │   ├── Category1_Model.vs.Contrl_MF_dagplot.png
## │   ├── Category1_Model.vs.Contrl_MF_idagplot.svg
## │   ├── Category1_Model.vs.Contrl_circ.csv
## │   ├── Category1_Model.vs.Contrl_ego.csv
## │   ├── Category1_Model.vs.Contrl_enrichlist.csv
## │   ├── Category1_Model.vs.Contrl_go_result.txt
## │   ├── Category1_Model.vs.Contrl_go_symbolid.csv
## │   └── Category1_Model.vs.Contrl_relist.csv
## └── 蛋白GO富集分析readme.png

  GO 功能显著性富集分析给出与所有鉴定到的蛋白质背景相比，差异蛋白质中显著富集的GO功能条目，从而给出差异蛋白质与哪些生物学功能显著相关。GO分为分子功能(Molecular function)、细胞组分(Cellular component)和生物过程(Biological process)三个部分。该分析首先把所有差异蛋白质向Gene Ontology数据库 (https://www.geneontology.org/)的各个term映射，计算每个term的蛋白质数目，然后应用超几何检验，找出与所有蛋白质背景相比，在差异蛋白质中显著富集的GO条目。其计算公式:

GO富集显著性计算公式
GO富集显著性计算公式

其中N为所有蛋白中具有GO注释信息的蛋白数目，n为N中差异蛋白的数目，M为所有蛋白中注释到某个GO条目的蛋白数目，x为注释到某个GO条目的差异蛋白数目。计算得到p值，以p值小于0.05为阈值，满足此条件的GO term定义为在差异蛋白质中显著富集的GO term。通过GO显著性分析能确定差异蛋白行使的主要生物学功能。差异蛋白GO富集结果见下表。

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html
注：
    1. GO ID: Gene Ontology数据库中唯一的标号信息
    2. Description: Gene Ontology功能的描述信息
    3. GeneRatio: 差异基因中与该Term相关的基因数与整个差异基因总数的比值
    4. BgRatio: 背景基因中与该Term相关的基因数与所有基因的比值
    5. pvalue: 富集分析统计学显著水平，一般情况下，p值小于0.05时该功能为显著富集项
    6. padjust: 校正后的p值
    7. qvalue: 对p值进行统计学检验的q值
    8. Count: 差异基因中与该Term相关的基因数

图 7.1.1 候选蛋白GO富集柱状图
图 7.1.1 候选蛋白GO富集柱状图

注： 图中显示的是p值极显著的前10个子功能，纵坐标是GO三个大类的下一层级的GO term，横坐标为注释到该term下（包括该term的子term）的候选蛋白个数。

图 7.1.2 候选蛋白GO富集气泡图
图 7.1.2 候选蛋白GO富集气泡图
注：
    1. geneRatio，是通过计算我们的目标基因富集到目的通路基因个数占目标基因总数（包含基因集总基因）的比值。
    2. 纵坐标轴是细胞生物学过程或者是一些发挥功能的通路。表示基因所对应的生物学过程或者分子功能是什么。
    3. Count是基因的个数，圈越大说明基因数目越多，反之越少。

图 7.1.3 候选蛋白GO cnetplot网络图
图 7.1.3 候选蛋白GO cnetplot网络图

注： 灰色的点代表蛋白，灰点下的字母数字表示该蛋白的Uniprot Accession(Uniprot上的蛋白编号，或者其他蛋白编号)，黄色的点代表富集到的GO terms，默认画top5富集到的GO terms，GO节点的大小对应富集到的蛋白个数

7.1.2 有向无环图

  有向无环图（Directed Acyclic Graph，DAG）为候选蛋白GO富集分析结果的图形化展示方式，分支代表包含关系，从上至下所定义的功能范围越来越小，一般选取GO富集分析的结果前10位作为有向无环图的主节点，并通过包含关系，将相关联的GO Term一起展示，矩形代表富集到的前10个GO terms，颜色从黄色过滤到红色，对应p值从大到小。我们的项目中分别绘制生物过程（biological process）、分子功能（molecular function）和细胞组分（cellular component）的候选蛋白DAG图（Directed Acyclic Graph）。

图 7.1.4 GO富集有向无环图
图 7.1.4 GO富集有向无环图
图 7.1.5 GO富集有向无环图（idagplot风格）
图 7.1.5 GO富集有向无环图（idagplot风格）

注： 每个方框或圆圈代表一个GO term，放大之后其中内容从上到下，代表的含义依次为:GO term的id、GO的描述、GO富集的p值、该GO下候选蛋白的数目/该GO下背景基因的前10的GO，颜色的深浅代表富集程度，颜色越深就表示富集程度越高。（两张图为一组表示两种风格的有向无环图（dagplot和idagplot），上图展示的是与分子功能（molecular function）相关的GO条目），isa表示父类与子类之间的关系。

  Cellular component 解释的是基因存在在哪里，在细胞质还是在细胞核？如果存在细胞质那在哪个细胞器上？如果是在线粒体中那是存在线粒体膜上还是在线粒体的基质当中？这些信息都叫Cellular component。 Biological process是在说明该基因参与了哪些生物学过程，比如，它参与了rRNA的加工或参与了DNA的复制，这些信息都叫Biological process。Molecular function在讲该基因在分子层面的功能是什么？它是催化什么反应的？

7.2 候选蛋白KEGG通路分析

点击打开文件夹

点击查看文件夹结构 
## ./4.kegg_enrichment
## ├── Category1_Model.vs.Contrl
## │   ├── Category1_Model.vs.Contrl_deg_upordown.csv
## │   ├── Category1_Model.vs.Contrl_enrichlist.csv
## │   ├── Category1_Model.vs.Contrl_kegg_barplot.svg
## │   ├── Category1_Model.vs.Contrl_kegg_dotplot.svg
## │   ├── Category1_Model.vs.Contrl_kegg_keggid.csv
## │   ├── Category1_Model.vs.Contrl_kegg_result.csv
## │   ├── Category1_Model.vs.Contrl_kegg_symbolid.csv
## │   ├── Category1_Model.vs.Contrl_keggid_list.csv
## │   ├── Category1_Model.vs.Contrl_symbolid_deg.csv
## │   ├── Chordal_kegg.png
## │   ├── kegg_chord.csv
## │   ├── kegg_circ.csv
## │   ├── kegg_circlize.svg
## │   ├── kegg_selete.csv
## │   ├── kegg_selete_all.csv
## │   └── pathway
## └── 蛋白KEGG富集分析readme.png

  在生物体内，不同基因相互协调行使其生物学功能，通过Pathway显著性富集能确定候选蛋白参与的最主要生化代谢途径和信号转导途径。KEGG（Kyoto Encyclopedia of Genesand Genomes）是有关Pathway的主要公共数据库4。Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出与整个基因组背景相比，在候选蛋白中显著性富集的Pathway。该分析的计算公式5：

KEGG通路富集显著性计算公式
KEGG通路富集显著性计算公式

其中，N为所有基因中具有Pathway注释的基因数目；n为N中候选蛋白的数目；M为所有基因中注释为某特定Pathway的基因数目；m为注释为某特定Pathway的候选蛋白数目。用BH的方法对p值进行校正，得到的校正后的p值越小代表越显著。这里将值小于0.05的Pathway定义为在候选蛋白中显著富集的Pathway。候选蛋白KEGG富集散点图是KEGG富集分析结果的图形化展示方式。

7.2.1 KEGG分析结果

注：
    1. ID: KEGG 数据库中通路唯一的编号信息。
    2. Description: Gene Ontology功能的描述信息
    3. GeneRatio: 差异基因中与该Term相关的基因数与整个差异基因总数的比值
    4. BgRatio: 背景基因中与该ID相关的基因数与所有基因的比值
    5. pvalue: 富集分析统计学显著水平，一般情况下，p值小于0.05时，该功能为富集项
    6. padjust: 校正后的p值
    7. qvalue: 对p值进行统计学检验的q值
    8. Count: 差异基因中与该Term相关的基因数

图 7.2.1 候选蛋白KEGG富集柱状图
图 7.2.1 候选蛋白KEGG富集柱状图

注： 图中显示的是p值极显著的前10个KEGG通路（颜色越红p值越小，越蓝p值越大），纵坐标是通路的描述，横坐标为注释到该通路下的候选蛋白个数。