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

一、概述

  TMT(Tandem Mass Tags，TMT)技术是由Thermo SCIENTIFIC公司研发的一种体外同重同位素标记的相对与绝对定量技术1。该技术利用多种同位素试剂标记蛋白多肽N末端或赖氨酸侧链基团，经高精度质谱仪串联分析，可同时比较多达16种样品之间的蛋白表达量，是近年来定量蛋白质组学常用的高通量筛选技术。下图为TMT定量技术原理图：

图 1.1.1 TMT定量技术原理图
图 1.1.1 TMT定量技术原理图

  标签TMT技术是一种使用标记剂的多肽定量方法，它基于质谱分析技术来实现定量测定。该技术通过直接比较样品中多肽的质谱特征，来推断其在不同样品中的相对含量。

二、实验流程

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

图 2.1.1 TMT实验流程图
图 2.1.1 TMT实验流程图

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

三、分析流程

  TMT无参蛋白质组分析云流程是微科盟针对缺少公共数据库注释的物种，或自建库推出的基于蛋白质序列功能注释的蛋白质组学分析云流程。该流程可以使用自建库，或其他自定义蛋白质数据（如多物种合并蛋白质库）作为数据源。基于质谱检测得到的Raw文件，进行对应数据库的搜索，然后基于数据库搜索的结果进行蛋白质鉴定，同时进行肽段、蛋白和母离子质量容差分布分析来评定质谱检测数据的质量；使用eggNOG-mapper对鉴定到的蛋白进行常见功能数据库注释，包括GO、KEGG、COG等注释信息；接下来进行蛋白质的定量分析，包括鉴定到的蛋白质表达量基础统计、差异蛋白的筛选、表达模式聚类分析；最后针对筛选出来的差异蛋白进行GO、KEGG功能富集分析和互作网络分析等一系列的差异蛋白功能分析。通过WGCNA分析描述不同样品之间蛋白质关联模式，分析高度协同变化基因集和表型之间的关系。对于多样本数据（单个分组30个以上样本），补充了ANOVA单因素方差分析用于寻找组间特异性表达的蛋白。TMT无参蛋白质组云流程解决了蛋白质组学分析中：自建库，物种合并库等分析困难的问题，以及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_PRO.vs.CTR
## │   ├── Category1_PRO.vs.CTR_All.xls
## │   ├── Category1_PRO.vs.CTR_All_symbol_id.xls
## │   ├── Category1_PRO.vs.CTR_DEG.xls
## │   ├── Category1_PRO.vs.CTR_DEG_symbol_id.xls
## │   ├── Category1_PRO.vs.CTR_Volcano.svg
## │   ├── Category1_PRO.vs.CTR_deeseq-density.png
## │   ├── heatmap.svg
## │   ├── input.xls
## │   └── table.xls
## ├── Vennupset
## │   └── venn_only1_group
## └── 差异蛋白分析readme.png

6.1 差异蛋白筛选（group1）

  蛋白质差异性分析能够极大地推动发现新的生物标志物，提升生物标志物鉴定的精准度，对分子作用机理、生物标志物、疾病早诊、分子分型、预后以及临床诊断等均具有重要价值。蛋白差异表达的输入数据为蛋白表达丰度表数据。分析时我们采用基于负二项分布的DESeq2进行分析4，具体使用基于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_PRO.vs.CTR
## │   ├── Category1_PRO.vs.CTR_ALL_barplot.svg
## │   ├── Category1_PRO.vs.CTR_ALL_cnetplot.svg
## │   ├── Category1_PRO.vs.CTR_ALL_dotplot.svg
## │   ├── Category1_PRO.vs.CTR_BP_dagplot.png
## │   ├── Category1_PRO.vs.CTR_BP_idagplot.svg
## │   ├── Category1_PRO.vs.CTR_CC_dagplot.png
## │   ├── Category1_PRO.vs.CTR_CC_idagplot.svg
## │   ├── Category1_PRO.vs.CTR_GO_Bubble.svg
## │   ├── Category1_PRO.vs.CTR_MF_dagplot.png
## │   ├── Category1_PRO.vs.CTR_MF_idagplot.svg
## │   ├── Category1_PRO.vs.CTR_circ.csv
## │   ├── Category1_PRO.vs.CTR_ego.csv
## │   ├── Category1_PRO.vs.CTR_enrichlist.csv
## │   ├── Category1_PRO.vs.CTR_go_result.txt
## │   ├── Category1_PRO.vs.CTR_go_symbolid.csv
## │   ├── Category1_PRO.vs.CTR_relist.csv
## │   ├── Category1_PRO.vs.CTR_symbol2accession.csv
## │   ├── Category1_PRO.vs.CTR_symbolid_deg.csv
## │   ├── GO_selete.csv
## │   ├── GO_selete_all.csv
## │   └── go_circlize.svg
## └── 蛋白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富集结果见下表。

注：
    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网络图

注： 灰色的点代表基因，灰点下的数字表示该基因的Entrez ID(NCBI上的基因编号)，黄色的点代表富集到的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_PRO.vs.CTR
## │   ├── Category1_PRO.vs.CTR_deg_upordown.csv
## │   ├── Category1_PRO.vs.CTR_enrichlist.csv
## │   ├── Category1_PRO.vs.CTR_kegg_barplot.svg
## │   ├── Category1_PRO.vs.CTR_kegg_dotplot.svg
## │   ├── Category1_PRO.vs.CTR_kegg_keggid.csv
## │   ├── Category1_PRO.vs.CTR_kegg_result.csv
## │   ├── Category1_PRO.vs.CTR_kegg_symbolid.csv
## │   ├── Category1_PRO.vs.CTR_keggid_list.csv
## │   ├── Category1_PRO.vs.CTR_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的主要公共数据库5。Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出与整个基因组背景相比，在候选蛋白中显著性富集的Pathway。该分析的计算公式6：

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值越大），纵坐标是通路的描述，横坐标为注释到该通路下的候选蛋白个数。

图 7.2.2 候选蛋白KEGG富集散点图
图 7.2.2 候选蛋白KEGG富集散点图

在此图中，KEGG富集程度通过Gene Ratio、padjust和富集到此通路上的基因个数来衡量。其中Gene Ratio指差异表达的基因中位于该Pathway条目的基因数目与差异基因中具有KEGG注释的基因总数的比值，Count表示富集到该通路的差异基因数量。Gene Ratio越大，表示富集的程度越大。padjust是校正后的p值。我们挑选了富集前20位的Pathway条目在该图中进行展示，若富集的Pathway条目不足20条，则全部展示。

7.2.2 KEGG富集通路代谢图

（结题报告篇幅有限，仅展示一个通路的代谢图，其余代谢图见附件）

注：
    1. K+num（基因ID号，表示在所有同源物种中具有相似结构或功能的一类同源蛋白）。如K01012=>生物素合成酶（备注：K建议大写）
    2. ko+num（代谢通路名称，表示一个特定的生物路径）如:ko0078=>生物素代谢 通路（备注：ko小写）
    3. M + num(模块名称）M00123=>生物素合成模块 C+num（化合物名）
    4. E -.-.-.-(酶名）EC2.1.1.116=>生物素合成酶（其实也就是k01012）
    5. R+num(反应名)
    6. RC+num(反应类型）
    7. RP+num（反应物质对）
    8. 图中红色表示根据log2FC计算得到的上调基因，蓝色表示下调基因
    9. 图中符号的详细解释参照KEGG官网（https://www.genome.jp/kegg/document/help_pathway.html）

7.3 差异蛋白互作分析

点击打开文件夹

点击查看文件夹结构 
## ./5.PPI_network
## ├── Category1_PRO.vs.CTR
## │   ├── cog_annotation.csv
## │   ├── cog_ppi.csv
## │   └── cog_ppi.html
## └── 蛋白PPI互作网络分析readme.png

  利用 StringDB7蛋白质互作数据库 (https://string-db.org/)对筛选的差异蛋白进行互作分析。在数据库中找到相应的物种，默认相互作用得分（combined）设置为400，得出相应的互作信息后，对分析结果进行网络图构建。差异蛋白互作结果见下表及下图。如果物种未被stringdb数据库收录将只返回蛋白对应的COG_PPI网络分析图。 ppi网络分析的动态显示涉及到浏览器渲染，请耐心等待，并尽可能将差异蛋白数量控制在200个以下避免图形过于复杂。

  COG，即Clusters of Orthologous Groups of proteins（同源蛋白簇）。是由NCBI创建并维护的蛋白数据库，根据细菌、藻类和真核生物完整基因组的编码蛋白系统进化关系分类构建而成。通过比对可以将某个蛋白序列注释到某一个COG中，每一簇COG由直系同源序列构成，从而可以推测该序列的功能。我们根据蛋白序列进行COG注释，使用差异分析结果进行筛选，然后映射到Stringdb数据库中COG互作网络中。

注：
    1. group1: 互作蛋白1
    2. group2: 互作蛋白2
    3. combined_score: 数据库中COG互作关系得分

点击查看差异蛋白COG注释信息
图 7.3.1 差异蛋白同源蛋白簇（COG）互作网络
点击查看PPI网络分析结果（Stringdb收录的物种有该部分结果） 注：
    1. from: 互作蛋白1
    2. to: 互作蛋白2
    3. combined_score: 数据库中通过实验方式证明互作关系得分

图 7.3.2 差异蛋白互作网络

注：当需要使用力导向图，或者调整网络图到窗口中部时可勾选右侧编辑栏中的enabled。

  将差异分析结果导入蛋白质互作网络分析中近一步得到蛋白质上下调互作网络：

图 7.3.3 差异蛋白上下调互作网络
注：图中红色圈表示表达量显著上调的蛋白，绿色表示显著下调的蛋白。

八、蛋白质结构域分析

点击打开文件夹

点击查看文件夹结构 
## ./6.Domain_analysis
## ├── Category1_PRO.vs.CTR
## │   ├── A0A0H2WVD4_domain.csv
## │   ├── A0A0H2WVD4结构域示意图.svg
## │   ├── A0A0H2WVE4_domain.csv
## │   ├── A0A0H2WVE4结构域示意图.svg
## │   ├── A0A0H2WWB3_domain.csv
## │   ├── A0A0H2WWB3结构域示意图.svg
## │   ├── A0A0H2WWP5_domain.csv
## │   ├── A0A0H2WWP5结构域示意图.svg
## │   ├── A0A0H2WWU3_domain.csv
## │   ├── A0A0H2WWU3结构域示意图.svg
## │   ├── A0A0H2WWU9_domain.csv
## │   ├── A0A0H2WWU9结构域示意图.svg
## │   ├── A0A0H2WX97_domain.csv
## │   ├── A0A0H2WX97结构域示意图.svg
## │   ├── A0A0H2WYR2_domain.csv
## │   ├── A0A0H2WYR2结构域示意图.svg
## │   ├── A0A0H2WZW1_domain.csv
## │   ├── A0A0H2WZW1结构域示意图.svg
## │   ├── A0A0H2WZW8_domain.csv
## │   ├── A0A0H2WZW8结构域示意图.svg
## │   ├── Q5HDV9_domain.csv
## │   ├── Q5HDV9结构域示意图.svg
## │   ├── Q5HEN0_domain.csv
## │   ├── Q5HEN0结构域示意图.svg
## │   ├── Q5HFA0_domain.csv
## │   ├── Q5HFA0结构域示意图.svg
## │   ├── Q5HFV0_domain.csv
## │   ├── Q5HFV0结构域示意图.svg
## │   ├── Q5HIM4_domain.csv
## │   └── Q5HIM4结构域示意图.svg
## └── 蛋白结构分析readme.png

  根据p值从小到大对蛋白质差异分析结果进行排序，选取前20个具有显著差异的蛋白。利用其Accession ID在uniprot上检索得到其结构域数据（表：8.1），在R4.3下利用drawProteins包对其进行可视化（图8.1.1）。结果展示了蛋白质结构域以及磷酸化位点（如果有磷酸化位点）。分析时需要在线检索，请确保你的蛋白质Accession被Uniprot收录。

注：
    1. type: 蛋白质结构域类型
    2. description: 对蛋白质结构域的描述，如该区域是某种酶的前体
    3. begin: 该结构域的开始位置位于第几个氨基酸
    4. end: 该结构域的结束位置位于第几个氨基酸
    5. length: 该结构域的长度
    6. accession: 蛋白的Accession ID
    7. entryName: 蛋白ID的简要名字
    8. taxid: 物种的编号

图 8.1.1 蛋白质结构示意图
图 8.1.1 蛋白质结构示意图
注：
    1. circles: 图中黄色圆圈表示磷酸化位点（如果有磷酸化位点）
    2. source: 数据来源于uniprot官网 (“https://www.uniprot.org/”)
    3. 横坐标表示: 氨基酸的位置
    4. 灰色矩形表示: 肽链
    5. 彩色方框表示: 该肽链上的结构域

九、差异蛋白理化性质分析

点击打开文件夹

点击查看文件夹结构 
## ./7.physicochemical_analysis
## ├── Category1_PRO.vs.CTR
## │   ├── Physicochemical_properties.csv
## │   ├── cog_bar.svg
## │   ├── cog_target.fasta
## │   ├── java_echart
## │   └── target.fasta
## └── 蛋白理化性质分析readme.png

  差异蛋白理化分析,根据差异蛋白id从uniprot数据库中下载相应的蛋白质序列数据。然后使用R语言包：Peptides基于序列数据进行理化性质预测。理化性质分析结果包括4个方面的内容：序列长度lengthpep，分子量mw，疏水性hydrophobicity，和等电点pI。分析结果以交互式散点图进行展示。

  COG（Clusters of Orthologous Groups ）注释是差异蛋白功能注释的一种方法。COG是由NCBI创建并维护的蛋白数据库，是对基因产物进行同源分类，为较早的识别直系同源基因的数据库.差异蛋白cog注释结果见下图：

图 9.1.2
图 9.1.2
注：
    1. 图例: COG功能注释分类和缩写
    2. 纵坐标: 蛋白质数量
    3. 横坐标: COG注释缩写

十、ANOVA单因素方差分析

点击打开文件夹

点击查看文件夹结构 
## ./8.Anova_analysis
## ├── Category1
## │   ├── ANOVA_bar_of_A0A0H2WW23.svg
## │   ├── ANOVA_bar_of_A0A0H2WWQ0.svg
## │   ├── ANOVA_bar_of_A0A0H2WX45.svg
## │   ├── ANOVA_bar_of_A0A0H2X0F6.svg
## │   ├── ANOVA_bar_of_Q5HEN0.svg
## │   ├── Category1_all_feature_anova_results.xls
## │   ├── Category1_all_significant_feature_barplot_of_duncan.svg
## │   ├── Category1_all_significant_feature_duncan_results.xls
## │   ├── Category1_anova_deg_result.xls
## │   ├── Category1_anova_deg_sig.xls
## │   ├── Domain
## │   │   ├── Q5HEN0_domain.csv
## │   │   └── Q5HEN0结构域示意图.svg
## │   ├── GO_enrich
## │   │   ├── Category1_ALL_barplot.svg
## │   │   ├── Category1_ALL_cnetplot.svg
## │   │   ├── Category1_ALL_dotplot.svg
## │   │   ├── Category1_BP_dagplot.png
## │   │   ├── Category1_BP_idagplot.svg
## │   │   ├── Category1_CC_dagplot.png
## │   │   ├── Category1_CC_idagplot.svg
## │   │   ├── Category1_MF_dagplot.png
## │   │   ├── Category1_MF_idagplot.svg
## │   │   ├── Category1_enrichlist.csv
## │   │   ├── Category1_go_result.txt
## │   │   ├── Category1_go_symbolid.csv
## │   │   ├── Category1_symbol2accession.csv
## │   │   ├── Category1_symbolid_deg.csv
## │   │   ├── GO_selete.csv
## │   │   ├── GO_selete_all.csv
## │   │   └── go_circlize.svg
## │   ├── KEGG_enrich
## │   │   ├── Category1_enrichlist.csv
## │   │   ├── Category1_kegg_symbolid.csv
## │   │   └── pathway
## │   ├── PPI_network
## │   │   ├── cog_annotation.csv
## │   │   ├── cog_ppi.csv
## │   │   └── cog_ppi.html
## │   ├── Physicochemical
## │   │   ├── Physicochemical_properties.csv
## │   │   ├── cog_bar.svg
## │   │   ├── cog_target.fasta
## │   │   ├── java_echart
## │   │   └── target.fasta
## │   ├── final_data_for_anova.xls
## │   └── heatmap
## │       ├── Category1_anova_deg_result.xls
## │       └── input.xls
## └── 蛋白单因素方差分析readme.png

  ANOVA单因素方差分析用于分析单一控制变量影响下的多组样本的均值是否存在显著性差异。通常用于单个分组样本数大于30，具有两组即两组以上分组的比较。ANOVA单因素方差分析要求数据符合正态分布，分析前默认使用scale(center=T,scale=T)（软件版本：R4.2）对蛋白质丰度数据进行正态转换。分析结果使用阈值padj<0.05进行筛选后绘制蛋白差异表格柱状图，差异蛋白GO、KEGG富集分析、PPI网络分析、结构域分析等下游结果详见结果文件夹(如有需要也可在云流程分步骤中挑选任意蛋白，作为差异蛋白)。

图 10.1.1
图 10.1.1
注：
    1. Group: 分组名
    2. Mean 纵坐标: 蛋白平均表达量，为分析过程中的因变量
    3. facter（Group）横坐标: 分组，因子变量，又称自变量

十一、WGCNA分析

点击打开文件夹

点击查看文件夹结构 
## ./9.WGCNA_analysis
## ├── Category1
## │   ├── All_hub_proteins.csv
## │   ├── black_hub_proteins.csv
## │   ├── brown_hub_proteins.csv
## │   ├── coral_hub_proteins.csv
## │   ├── darkgrey_hub_proteins.csv
## │   ├── darkred_hub_proteins.csv
## │   ├── darksalmon_hub_proteins.csv
## │   ├── figures
## │   ├── firebrick_hub_proteins.csv
## │   ├── gainsboro_hub_proteins.csv
## │   ├── lightcoral_hub_proteins.csv
## │   ├── lightgrey_hub_proteins.csv
## │   ├── maroon_hub_proteins.csv
## │   ├── mistyrose_hub_proteins.csv
## │   ├── red_hub_proteins.csv
## │   ├── rosybrown_hub_proteins.csv
## │   ├── salmon_hub_proteins.csv
## │   ├── sample_info.csv
## │   ├── snow_hub_proteins.csv
## │   ├── tomato_hub_proteins.csv
## │   ├── tpm.csv
## │   ├── whitesmoke_hub_proteins.csv
## │   └── 模块表达量矩阵.csv
## └── 蛋白WGCNA分析readme.png

  加权基因共表达网络分析 (WGCNA, Weighted correlation network analysis)是用来描述不同样品之间基因关联模式的系统生物学方法，可以用来鉴定高度协同变化的基因集, 并根据基因集的内连性和基因集与表型之间的关联鉴定候补生物标记基因或治疗靶点。

  相比于只关注差异表达的基因，WGCNA利用全部基因的信息识别感兴趣的基因集，并与表型进行显著性关联分析。一是充分利用了信息，二是把数千个基因与表型的关联转换为数个基因集与表型的关联，免去了多重假设检验校正的问题。

  共表达网络定义为加权基因网络。点代表基因，边代表基因表达相关性。加权是指对相关性值进行幂次运算(幂次的值也就是软阈值 (power, model_fit.svg展示了确定合适的power的过程))。无符号(unsigned)网络的边属性计算方式为abs(cor(genex, geney)) ^ power；有符号(signed)网络的边属性计算方式为(1+cor(genex, geney)/2) ^ power; sign hybrid的边属性计算方式为cor(genex, geney)^power if cor>0 else 0。这种处理方式强化了强相关，弱化了弱相关或负相关，使得相关性数值更符合无标度网络特征，更具有生物意义。

  WGCNA分析将共表达网络转换为TOM拓扑重叠矩阵 (Topological overlap matrix)，以降低噪音和假相关，获得的新距离矩阵。当两个基因在共表达网络中拥有较多相同的”邻居”时，我们说这两个基因拓扑重叠

  WGCNA分析利用TOM拓扑重叠矩阵对基因进行聚类，将拓扑重叠的基因聚为一个模块。Module(模块)指的是高度內连的基因集。在无符号(unsigned)网络中，模块内是高度相关的基因。在有符号网络中，模块内是高度正相关的基因。WGCNA的分析过程可参考图11.1.1。

图 11.1.1 WGCNA分析流程示意图
图 11.1.1 WGCNA分析流程示意图

11.1 鉴定模块

  WGNCA分析使用Python的PyWGCNA模块进行8，默认使用蛋白质表达量丰度数据protein_count.text作为输入。数据通过sklearn包的preprocessing.MinMaxScaler()函数进行标准化。通过WGCNA流程对基因进行聚类：1.选择软阈值2.计算邻接矩阵3.计算TOM相似矩阵4.对dissTOM的基因进行聚类5.动态合并相似聚类。   图11.1.3展示了不同power(幂次)下，WGCNA拓扑模型的拟合度，拟合度绝对值高于0.9时的power作为最终加权幂次，以及不同power(幂次)下，加权共表达网络中基因的平均连接度。连接度太高，不利于模块的区分，可能导致最终结果里只有一个模块，发生这种现象的原因主要有三个:   (1) 有离群样本   (2) 各分组间样本差异太大(批次效应)   (3) 样本数太少。样本少时，基因之间倾向于有较高的相关系数值。举个极端例子，当只有两个样本时，任意两个基因之间的相关系数都是1，无论怎么分，只可能有一个模块。通常样本数多于15个时，才能有较好的模块鉴定结果。

图 11.1.2 不同power下拓扑模型的拟合度和基因平均连接度
图 11.1.2 不同power下拓扑模型的拟合度和基因平均连接度

  图11.1.3展示了各个模块基因根据TOM拓扑重叠矩阵聚类的结果。只有高度差(近似’1-拓扑重叠指数’得到的距离)在高度阈值(默认0.2)以内的基因，才有可能被归为一个模块。grey是无法归类到任何模块的基因的集合。

图 11.1.3 模块基因聚类树
图 11.1.3 模块基因聚类树

  基于给定模块中的连通性,查找前10个枢纽基因作为核心基因，（云流程中可以调整作为核心基因的数量，如需查看各模块全部结果可将阈值设置为蛋白质总数），如表11.1所示。并通过网络图来展示它们的连接关系,如图11.1.4所示。默认只会显示连接度在100以上的蛋白，如需要调整可使用云流程分步骤进行。

图 11.1.4 Hub核心基因（蛋白）网络图

11.2 模块关联表型

  我们将老师提供的分组变量作为表型因素（如年龄0到20一组，20到40一组，40到60一组，且要求分组数量大于等于3，否则图11.2.3将不会展示各分组名称），我们可以将模块的Module eigengene E（同上表11.1）与表型数据做相关分析，从而推断模块影响的表型。注意：该分析部分使用是属于某模块的全部基因（蛋白），而不是Hub基因（蛋白），请注意区分。   图11.2.1用柱状图的形式展示了模块本征基因（属于该模块的基因，在本分析中为属于该模块的蛋白）在不同分组中的表达量。

图 11.2.1 模块征基因柱状图
图 11.2.1 模块征基因柱状图
注：
    1. 图例: 分组情况
    2. 纵坐标: 模块本征基因（蛋白质）表达量的平均值（经过标准化处理）
    3. 横坐标: 分组名

  图11.2.2用热图的形式展示了模块本征基因（属于该模块的基因，在本分析中为属于该模块的蛋白）在不同分组中的表达量。相较于柱状图，可以看到模块本征基因（蛋白）在各组间的表达模式。理想状态下同一分组下同一模块的本征基因（蛋白）的表达模式（即热图的颜色分布情况）应该是相似的，同一模块的本征基因在不同分组的表达模式相对而言是具有较大差异的（和实验设计，组间差异程度，以及模块的选择有关）。

图 11.2.2 模块征基因热图
图 11.2.2 模块征基因热图
注：
    1. 图例: 分组情况
    2. 纵坐标: 模块本征基因（蛋白质）表达量的平均值（经过标准化处理）
    3. 横坐标: 样本
    4. 热图中颜色越红代表该基因（蛋白）的表达量越高，越黑表达量越低。

  图11.2.3展示的是模型-表型相关性热图。用于推断模块对表型的影响。要查看各样本与模块的原始相关性矩阵可查看文件夹中“样本-模块相关性矩阵.csv”文件。

图 11.2.3 模型-表型相关性热图
图 11.2.3 模型-表型相关性热图
注：
    1. 图例: 相关性范围值，蓝色表示负相关，红色表示正相关
    2. 纵坐标: 分组名（分组小于3组是会显示category名称）
    3. 横坐标: 模块名称,括号中的数字表示该模块包含的基因（蛋白）数量
    4. 热图中数字表示相关性系数（pearson），括号中的数字表示对应的P值。

十二、分析所需软件版本

  分析用到的软件和数据库版本如下表所示：
软件 版本
DESeq2 1.40.0
edgeR 3.42.1
clusterProfiler 4.8.1
stringdb 2.12.0
StringDB v11.5
r-base 4.3.1
python 3.8
annotationhub 3.8.0
emapper 2.0.1
pyvis 0.3.1
PyWGCNA 1.20.4

参考文献

  1. Wu J, An Y, Pu H, Shan Y, Ren X, An M, Wang Q, Wei S, Ji J. Enrichment ofserum low-molecular-weight proteins using C18 absorbent under urea/dithiothreitoldenatured environment. Anal Biochem. 2010 Mar 1;398(1):34-44.↩︎

  2. Wu J, Xie X, Liu Y, He J, Benitez R, Buckanovich RJ, Lubman DM.Identification and confirmation of differentially expressed fucosylated glycoproteinsin the serum of ovarian cancer patients using a lectin array and LC-MS/MS. JProteome Res. 2012 Sep 7;11(9):4541-52.↩︎

  3. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP,Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A,Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G. Geneontology: tool for the unification of biology. The Gene Ontology Consortium. NatGenet. 2000 May;25(1):25-9.↩︎

  4. Love, M. I. , Huber, W. , & Anders, S. . (2014). Moderated estimation of fold change and dispersion for rna-seq data with deseq2. Genome Biology, 15(12), 550.↩︎

  5. Smyth, G. K. . (2010). edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 26(1), 139.↩︎

  6. Kanehisa M, Goto S, Kawashima S, Okuno Y, Hattori M. The KEGG resourcefor deciphering the genome. Nucleic Acids Res. 2004 Jan 1;32(Databaseissue):D277-80.↩︎

  7. Damian, S. , Andrea, F. , Stefan, W. , Kristoffer, F. , Davide, H. , & Jaime, H. C. , et al. (2015). String v10: protein–protein interaction networks, integrated over the tree of life. Nucleic Acids Research, 43(D1). Nucleic Acids Res. 2014 Jan;42(Database issue):D222-30.↩︎

  8. Narges Rezaie, Farilie Reese, Ali Mortazavi, PyWGCNA: a Python package for weighted gene co-expression network analysis, Bioinformatics, Volume 39, Issue 7, July 2023, btad415, https://doi.org/10.1093/bioinformatics/btad415.↩︎