蛋白质组与蛋白质分析

Author

Guangchuang Yu
School of Basic Medical Sciences, Southern Medical University

Published

October 21, 2026

实验基于R语言环境，安装R语言，可参考：

R语言安装

本实验以一组与DNA损伤应答相关的蛋白为例，用R完成一条”从蛋白序列 → 基本属性 → 相互作用网络 → 功能富集 → 质谱定量（拓展）”的分析链。涉及的蛋白：p53 (TP53)、MDM2、BRCA1、BRCA2、ATM、CHK2、CDKN2A。

一、蛋白序列下载与基本属性

用 rentrez 从 NCBI 下载序列

本节直接从 NCBI 下载这 7 个蛋白的 RefSeq 蛋白序列。给定蛋白的 RefSeq accession，用 rentrez::entrez_fetch() 获取 FASTA；NCBI 页面：https://www.ncbi.nlm.nih.gov/protein。id 接受向量，一次请求就能取回全部记录（返回一个拼好的多记录 FASTA 字符串），不必逐条循环。

pload(rentrez)

── loading the package: rentrez ──

acc <- c(TP53  = "NP_000537.3",
         MDM2  = "NP_002383.2",
         BRCA1 = "NP_009225.1",
         BRCA2 = "NP_001419006.1",
         ATM   = "NP_000042.3",
         CHEK2 = "NP_001005735.1",
         CDKN2A= "NP_000068.1")

if (file.exists("proteins_ncbi.fasta")) {
  message("proteins_ncbi.fasta 已存在，跳过 NCBI 下载")
} else {
  fa <- entrez_fetch(db = "protein", id = acc, rettype = "fasta")
  cat(fa, file = "proteins_ncbi.fasta")
}

proteins_ncbi.fasta 已存在，跳过 NCBI 下载

项目目录里另附两份离线文件：proteins_ncbi.fasta 是按上面 7 个 RefSeq accession 下载好的序列（与在线路径完全一致，NCBI 暂时不可用时直接用它）；proteins.fasta 来自 UniProt，注意它选的是各自数据库的代表序列，MDM2（491 aa）与 CHEK2（543 aa，isoform a）与 RefSeq（497 aa / 586 aa，isoform c）不同，长度、分子量、等电点都会差一点。

pload(Biostrings)

── loading the package: Biostrings ──

Loading required package: BiocGenerics

Loading required package: generics


Attaching package: 'generics'

The following objects are masked from 'package:base':

    as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
    setequal, union


Attaching package: 'BiocGenerics'

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    anyDuplicated, aperm, append, as.data.frame, basename, cbind,
    colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
    get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
    rbind, Reduce, rownames, sapply, saveRDS, table, tapply, unique,
    unsplit, which.max, which.min

Loading required package: S4Vectors

Loading required package: stats4


Attaching package: 'S4Vectors'

The following object is masked from 'package:utils':

    findMatches

The following objects are masked from 'package:base':

    expand.grid, I, unname

Loading required package: IRanges

Loading required package: XVector

Loading required package: Seqinfo


Attaching package: 'Biostrings'

The following object is masked from 'package:base':

    strsplit

seqs <- readAAStringSet("proteins_ncbi.fasta")
names(seqs) <- names(acc)   # 按下载顺序对应到基因名
seqs

AAStringSet object of length 7:
    width seq                                               names               
[1]   393 MEEPQSDPSVEPPLSQETFSDLW...SKKGQSTSRHKKLMFKTEGPDSD TP53
[2]   497 MVRSRQMCNTNMSVPTDGAVTTS...KKRNKPCPVCRQPIQMIVLTYFP MDM2
[3]  1863 MDLSALRVEEVQNVINAMQKILE...SVALYQCQELDTYLIPQIPHSHY BRCA1
[4]  3418 MPIGSKERPTFFEIFKTRCNKAD...SSQASTEECEKNKQDTITTKKYI BRCA2
[5]  3056 MSLVLNDLLICCRQLEHDRATER...LLIQQAIDPKNLSRLFPGWKAWV ATM
[6]   586 MSRESDVEAQQSHGSSACSQPHG...RPREGEAEGAETTKRPAVCAAVL CHEK2
[7]   156 MEPAAGSSMEPSADWLATAAARG...AGGTRGSNHARIDAAEGPSDIPD CDKN2A

蛋白长度与氨基酸组成

len <- width(seqs)
data.frame(gene = names(seqs), length = len)

    gene length
1   TP53    393
2   MDM2    497
3  BRCA1   1863
4  BRCA2   3418
5    ATM   3056
6  CHEK2    586
7 CDKN2A    156

aa <- alphabetFrequency(seqs, baseOnly = TRUE)
rownames(aa) <- names(seqs)
round(aa[ , apply(aa, 2, function(x) sd(x) > 0)], 2)

         A   R   N   D  C   Q   E   G  H   I   L   K  M   F   P   S   T  W  Y
TP53    24  26  14  20 10  15  30  23 12   8  32  20 12  11  45  38  22  4  9
MDM2    15  23  18  36 19  28  56  21  9  20  36  31  9  11  23  67  25  4 14
BRCA1   84  76 121  85 44  97 198  87 49  77 156 137 30  49  96 224 111 10 31
BRCA2  174 110 230 171 76 154 293 122 83 187 281 322 45 136 149 381 220 20 79
ATM    157 152 137 156 89 154 233 119 72 179 389 214 79 123  98 234 158 38 96
CHEK2   36  29  20  29 14  24  51  33 14  27  56  40  7  21  30  63  32  6 18
CDKN2A  29  17   4  11  1   1  11  14  5   3  17   0  5   1  11   7   5  2  2
         V
TP53    18
MDM2    32
BRCA1  101
BRCA2  185
ATM    179
CHEK2   36
CDKN2A  10

分子量与等电点（`Peptides`）

用 Peptides 计算分子量与等电点。pI() 默认 pKscale = "EMBOSS"；与 ExPASy Compute pI/Mw（Bjellqvist 标度）对照时会差 0.2–0.5。mw()、pI() 都直接吃整个字符向量、返回等长向量，不用逐条循环；但必须先 as.character()——把 AAStringSet 整个传进去会被当成一条序列，只算出水的质量。

pload(Peptides)

── loading the package: Peptides ──

res <- data.frame(
  gene = names(seqs),
  length = len,
  mw  = mw(as.character(seqs)),
  pI  = pI(as.character(seqs))
)
res[, c("mw", "pI")] <- round(res[, c("mw", "pI")], 1)
res

         gene length       mw  pI
TP53     TP53    393  43653.2 6.8
MDM2     MDM2    497  55990.8 4.4
BRCA1   BRCA1   1863 207720.8 5.1
BRCA2   BRCA2   3418 384230.3 6.7
ATM       ATM   3056 350687.5 6.8
CHEK2   CHEK2    586  65418.9 6.1
CDKN2A CDKN2A    156  16532.6 5.6

思考：BRCA2、ATM 为何”又长又重”？pI 与蛋白在细胞中的定位/功能有何关系？

二、蛋白相互作用网络（mechgraph + ggtangle）

mechgraph 从数据库下载并构建”机制证据图”（mechgraph 对象，含节点/边表）；ggtangle 直接接收 mechgraph 对象——ggplot(mg) 自动把节点/边转成网络图，边用 geom_edge()，节点 size 自动为度（degree）。

数据文件

文件	说明
`proteins_ncbi.fasta`	上表 7 个 RefSeq accession 的序列
`proteins.fasta`	同 7 个蛋白的 UniProt 代表序列；MDM2 491 aa、CHEK2 543 aa（isoform a），与 RefSeq（497 / 586 aa，isoform c）不同
`BIOGRID-MV-Physical-LATEST.tab3.zip`	BioGRID Multi-Validated physical 互作数据（35 MB 压缩 / 263 MB 解压），供 `mg_from_biogrid(file = ...)` 使用

mg_from_string() 可改用 STRING（keytype = "symbol", score_threshold = 700），该数据需在线下载。

用 mechgraph 下载网络 + 聚焦到我们的基因

pload(mechgraph)

── loading the package: mechgraph ──

mg <- mg_from_biogrid(file = "BIOGRID-MV-Physical-LATEST.tab3.zip",
                      taxon_id = 9606, dataset = "mv_physical")
mg

<mechgraph>
  nodes: 12041
  edges: 101036
  node types: protein
  edge types: ppi

genes <- c("TP53", "MDM2", "BRCA1", "BRCA2", "ATM", "CHEK2", "CDKN2A")
nd <- mg_nodes(mg)
ids <- nd$id[nd$label %in% genes]   # 节点 id 是 BIOGRID:xxxxx，基因名在 label 列
mg2 <- mg_induced_subgraph(mg, ids)
mg_nodes(mg2)

                 id    type  label taxon_id
18   BIOGRID:107140 protein  BRCA1     9606
203  BIOGRID:107142 protein  BRCA2     9606
237  BIOGRID:106962 protein    ATM     9606
245  BIOGRID:113010 protein   TP53     9606
2671 BIOGRID:107463 protein CDKN2A     9606
3142 BIOGRID:116369 protein  CHEK2     9606
3248 BIOGRID:110358 protein   MDM2     9606

mg_edges(mg2)

                from             to type  source
5039  BIOGRID:106962 BIOGRID:106962  ppi BioGRID
5043  BIOGRID:106962 BIOGRID:107140  ppi BioGRID
5065  BIOGRID:106962 BIOGRID:110358  ppi BioGRID
5089  BIOGRID:106962 BIOGRID:113010  ppi BioGRID
5116  BIOGRID:106962 BIOGRID:116369  ppi BioGRID
6826  BIOGRID:107140 BIOGRID:107140  ppi BioGRID
6827  BIOGRID:107142 BIOGRID:107140  ppi BioGRID
6835  BIOGRID:107463 BIOGRID:107140  ppi BioGRID
6978  BIOGRID:107140 BIOGRID:113010  ppi BioGRID
7058  BIOGRID:116369 BIOGRID:107140  ppi BioGRID
7166  BIOGRID:107142 BIOGRID:107142  ppi BioGRID
7186  BIOGRID:107142 BIOGRID:113010  ppi BioGRID
7203  BIOGRID:116369 BIOGRID:107142  ppi BioGRID
11932 BIOGRID:107463 BIOGRID:107463  ppi BioGRID
11946 BIOGRID:110358 BIOGRID:107463  ppi BioGRID
11960 BIOGRID:107463 BIOGRID:113010  ppi BioGRID
39314 BIOGRID:110358 BIOGRID:110358  ppi BioGRID
39385 BIOGRID:110358 BIOGRID:113010  ppi BioGRID
39452 BIOGRID:116369 BIOGRID:110358  ppi BioGRID
68186 BIOGRID:113010 BIOGRID:113010  ppi BioGRID
68377 BIOGRID:116369 BIOGRID:113010  ppi BioGRID
90490 BIOGRID:116369 BIOGRID:116369  ppi BioGRID
                                                                                                                                                                                                                                                                                            evidence
5039                                                                                                                                                                        Protein-peptide; Affinity Capture-MS; Biochemical Activity; Co-purification; Reconstituted Complex; Co-crystal Structure
5043                                                                                                                                                                     Protein-peptide; Affinity Capture-Western; Reconstituted Complex; Co-purification; Biochemical Activity; Proximity Label-MS
5065                                                                                                                                                                                                                                                           Biochemical Activity; Co-localization
5089                                                                                                                                                                                         Reconstituted Complex; Protein-peptide; Biochemical Activity; Co-localization; Affinity Capture-Western
5116                                                                                                                                                                                                                                                     Biochemical Activity; Reconstituted Complex
6826                                                                                                                         Protein-peptide; Co-crystal Structure; Reconstituted Complex; Biochemical Activity; Affinity Capture-Western; Affinity Capture-MS; Proximity Label-MS; Co-fractionation
6827                                                                                                                                                                           Affinity Capture-Western; Two-hybrid; Co-purification; Reconstituted Complex; Affinity Capture-MS; Proximity Label-MS
6835                                                                                                                                                                                                                                                       Affinity Capture-Western; Co-localization
6978                                                                                                                                                                                                           Affinity Capture-Western; Affinity Capture-MS; Reconstituted Complex; Co-purification
7058                                                                                                                                                                                                          Reconstituted Complex; Affinity Capture-Western; Biochemical Activity; Co-localization
7166                                                                                                                                                                                                           Co-localization; Affinity Capture-MS; Reconstituted Complex; Affinity Capture-Western
7186                                                                                                                                                                                       Affinity Capture-Western; Affinity Capture-MS; Co-purification; Reconstituted Complex; Proximity Label-MS
7203                                                                                                                                                                                                                                 Affinity Capture-Western; Protein-peptide; Biochemical Activity
11932                                                                                                                                                                                                                                                                    Reconstituted Complex; FRET
11946                                                                                                                                                                                        Affinity Capture-MS; Affinity Capture-Western; Two-hybrid; Reconstituted Complex; Protein-peptide; FRET
11960                                                                                                                                                                                                       Affinity Capture-Western; Co-localization; FRET; Affinity Capture-MS; Proximity Label-MS
39314                                                                                                          Biochemical Activity; Affinity Capture-Western; Co-crystal Structure; Reconstituted Complex; Two-hybrid; Co-fractionation; Protein-peptide; Affinity Capture-MS; Co-localization; PCA
39385 Biochemical Activity; Affinity Capture-Western; Reconstituted Complex; Co-fractionation; Affinity Capture-MS; Co-crystal Structure; PCA; Co-localization; Affinity Capture-RNA; Two-hybrid; Protein-peptide; Protein-RNA; Far Western; FRET; Affinity Capture-Luminescence; Proximity Label-MS
39452                                                                                                                                                                                                                          Biochemical Activity; Affinity Capture-Western; Reconstituted Complex
68186                                                                                                                                                Reconstituted Complex; Co-purification; Affinity Capture-Western; Two-hybrid; FRET; Co-fractionation; Co-crystal Structure; Affinity Capture-MS
68377                                                                                                                                                                                                                          Affinity Capture-Western; Biochemical Activity; Reconstituted Complex
90490                                                                                                                                                                  Affinity Capture-Western; Reconstituted Complex; Biochemical Activity; Co-crystal Structure; Co-purification; Protein-peptide
      evidence_type pmid n_evidence
5039       physical              13
5043       physical              12
5065       physical               3
5089       physical              34
5116       physical               7
6826       physical              27
6827       physical              15
6835       physical               3
6978       physical              23
7058       physical              10
7166       physical               5
7186       physical               9
7203       physical               3
11932      physical               2
11946      physical              48
11960      physical              17
39314      physical              66
39385      physical             666
39452      physical               7
68186      physical              34
68377      physical               9
90490      physical              19

BioGRID 同一对蛋白之间常有多条”证据边”，还可能有反向重复对与自环，画图前可先化简，观察会更清晰： g <- igraph::simplify(igraph::graph_from_data_frame(mg_edges(mg2)[, c("from","to")], directed = FALSE))

用 ggtangle 直接可视化 mechgraph

pload(ggtangle)

── loading the package: ggtangle ──

ggtangle v0.1.4 Learn more at https://yulab-smu.top/

pload(ggplot2)

── loading the package: ggplot2 ──

pload(shadowtext)

── loading the package: shadowtext ──

ggplot(mg2, layout = "circular") +
    geom_edge(alpha = .4, color = "#9AA5B1") +
    geom_point(aes(size = size), color = "#0B4F6C", show.legend = FALSE) +
    geom_shadowtext(aes(label = name), color = "black",
                    bg.color = "white", bg.r = 0.18, size = 3.5) +
    theme_void() +
    theme(plot.background = element_rect(fill = "white", color = NA))

思考：图中最大的节点（度最高）是哪个蛋白？为何它在 DNA 损伤应答中如此核心？（这组基因在 BioGRID 中高度连通，TP53 是最核心的 hub。）

需要具体的度数值时，可用 igraph::as.igraph(mg2) 把 mechgraph 转成 igraph（节点/边信息会作为 igraph 属性保留）。注意 BioGRID 的 physical 数据里自相互作用（同源二聚体）也是一条边，igraph 会把自环算作 2 度；算度数前先用 g <- igraph::simplify(igraph::as.igraph(mg2), remove.multiple = TRUE, remove.loops = TRUE)，再 igraph::degree(g)。若改用 STRING，只需把下载语句换成：mg <- mg_from_string(taxon_id = 9606, keytype = "symbol", score_threshold = 700)。

三、功能富集分析

这组蛋白”一起做什么”？用它们的基因做 GO 富集。

entrez <- c(TP53 = 7157, MDM2 = 4193, BRCA1 = 672, BRCA2 = 675,
            ATM = 472, CHEK2 = 11200, CDKN2A = 1029)

pload(clusterProfiler)

── loading the package: clusterProfiler ──

clusterProfiler v4.21.3 Learn more at https://yulab-smu.top/contribution-knowledge-mining/

Please cite:

G Yu. Thirteen years of clusterProfiler. The Innovation. 2024,
5(6):100722


Attaching package: 'clusterProfiler'

The following object is masked from 'package:XVector':

    slice

The following object is masked from 'package:IRanges':

    slice

The following object is masked from 'package:S4Vectors':

    rename

The following object is masked from 'package:stats':

    filter

pload(org.Hs.eg.db)

── loading the package: org.Hs.eg.db ──

Loading required package: AnnotationDbi

Loading required package: Biobase

Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.


Attaching package: 'AnnotationDbi'

The following object is masked from 'package:clusterProfiler':

    select

x <- enrichGO(entrez, OrgDb = org.Hs.eg.db, ont = "BP",
              pvalueCutoff = 0.2, qvalueCutoff = 0.2, readable = TRUE)
x

#
# over-representation test
#
#...@organism    Homo sapiens 
#...@ontology    BP 
#...@keytype     ENTREZID 
#...@gene    chr [1:7] "7157" "4193" "672" "675" "472" "11200" "1029"
#...pvalues adjusted by 'BH' with cutoff < 0.2
#...629 enriched terms found
'data.frame':   629 obs. of  13 variables:
 $ ID            : chr  "GO:0042770" "GO:0090399" "GO:0071479" "GO:0030330" ...
 $ Description   : chr  "signal transduction in response to DNA damage" "replicative senescence" "cellular response to ionizing radiation" "DNA damage response, signal transduction by p53 class mediator" ...
 $ GeneRatio     : chr  "7/7" "4/7" "5/7" "5/7" ...
 $ BgRatio       : chr  "185/18842" "14/18842" "63/18842" "72/18842" ...
 $ RichFactor    : num  0.0378 0.2857 0.0794 0.0694 0.0459 ...
 $ FoldEnrichment: num  102 769 214 187 123 ...
 $ oddsRatio     : num  Inf 2510 809 700 450 ...
 $ zScore        : num  26.6 55.4 32.6 30.5 24.7 ...
 $ pvalue        : num  7.85e-15 6.66e-12 7.42e-12 1.48e-11 1.23e-10 ...
 $ p.adjust      : num  4.94e-12 1.56e-09 1.56e-09 2.32e-09 1.55e-08 ...
 $ qvalue        : num  1.09e-12 3.44e-10 3.44e-10 5.13e-10 3.42e-09 ...
 $ geneID        : chr  "TP53/MDM2/BRCA1/BRCA2/ATM/CHEK2/CDKN2A" "TP53/ATM/CHEK2/CDKN2A" "TP53/MDM2/BRCA1/BRCA2/ATM" "TP53/MDM2/ATM/CHEK2/CDKN2A" ...
 $ Count         : int  7 4 5 5 5 5 5 6 5 5 ...
#...Citation

dotplot(x)

cnetplot(x, showCategory = 5)

思考：富集到的最显著生物学过程是什么？与你对这组蛋白功能的直觉一致吗？

四、拓展：质谱定量分析（MSstats）

如果拿到质谱（如 MaxQuant）的定量数据，可用 MSstats 做蛋白水平的差异丰度分析。

用 MSstats 做蛋白差异丰度分析

pload(MSstats)

── loading the package: MSstats ──


Attaching package: 'MSstats'

The following object is masked from 'package:grDevices':

    savePlot

data(DDARawData)   # 内置示例：label-free DDA 掺入（spike-in）实验，6 蛋白 / 6 条件（Mueller et al. 2007，latin square 设计）
head(DDARawData)

  ProteinName PeptideSequence PrecursorCharge FragmentIon ProductCharge
1      bovine     S.PVDIDTK_5               5          NA            NA
2      bovine     S.PVDIDTK_5               5          NA            NA
3      bovine     S.PVDIDTK_5               5          NA            NA
4      bovine     S.PVDIDTK_5               5          NA            NA
5      bovine     S.PVDIDTK_5               5          NA            NA
6      bovine     S.PVDIDTK_5               5          NA            NA
  IsotopeLabelType Condition BioReplicate Run Intensity
1                L        C1            1   1   2636792
2                L        C1            1   2   1992418
3                L        C1            1   3   1982146
4                L        C2            1   4   5019594
5                L        C2            1   5   4560468
6                L        C2            1   6   3627849

dim(DDARawData)

[1] 2070   10

table(DDARawData$ProteinName)


   bovine   chicken cyc_horse myg_horse    rabbit     yeast 
      252       198       576       216       558       270

processed <- dataProcess(DDARawData)   # 肽段→蛋白水平汇总

INFO  [2026-10-21 11:37:34] ** Features with one or two measurements across runs are removed.
INFO  [2026-10-21 11:37:34] ** Fractionation handled.
INFO  [2026-10-21 11:37:34] ** Updated quantification data to make balanced design. Missing values are marked by NA
INFO  [2026-10-21 11:37:34] ** Log2 intensities under cutoff = 13.456  were considered as censored missing values.
INFO  [2026-10-21 11:37:34] ** Log2 intensities = NA were considered as censored missing values.
INFO  [2026-10-21 11:37:34] ** Use top100 features that have highest average of log2(intensity) across runs.
INFO  [2026-10-21 11:37:34] 
 # proteins: 6
 # peptides per protein: 11-32
 # features per peptide: 1-1
INFO  [2026-10-21 11:37:34] 
                    C1 C2 C3 C4 C5 C6
             # runs  3  3  3  3  3  3
    # bioreplicates  1  1  1  1  1  1
 # tech. replicates  3  3  3  3  3  3
INFO  [2026-10-21 11:37:34] Some features are completely missing in at least one condition:  
 -.PHSHPALTPEQK_35_35_NA_NA,
 A.DIAGHGQEVLIR_321_321_NA_NA,
 A.PGFTYTDANK_129_129_NA_NA,
 A.PGFTYTDANK_130_130_NA_NA,
 A.PGFTYTDANK_131_131_NA_NA ...
INFO  [2026-10-21 11:37:34]  == Start the summarization per subplot...

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |======================================================================| 100%
INFO  [2026-10-21 11:37:34]  == Summarization is done.

names(processed)

[1] "FeatureLevelData" "ProteinLevelData" "SummaryMethod"

构造对比矩阵（比较 C6 与 C4 两个条件；latin square 设计下各蛋白浓度走向不同，log2FC 有正有负）：

lv <- c("C1", "C2", "C3", "C4", "C5", "C6")
cm <- matrix(0, nrow = 1, ncol = 6, dimnames = list("C6-vs-C4", lv))
cm[1, "C4"] <- -1
cm[1, "C6"] <- 1

cmp <- groupComparison(contrast.matrix = cm, data = processed)

INFO  [2026-10-21 11:37:34]  == Start to test and get inference in whole plot ...

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |======================================================================| 100%
INFO  [2026-10-21 11:37:34]  == Comparisons for all proteins are done.

res <- cmp$ComparisonResult   # 每蛋白的 log2FC、p 值等
res

    Protein    Label    log2FC        SE     Tvalue DF       pvalue
1    bovine C6-vs-C4  3.644656 0.4245943   8.583853 11 3.322944e-06
2   chicken C6-vs-C4 -5.108220 0.2205455 -23.161746 12 2.498091e-11
3 cyc_horse C6-vs-C4  2.691865 0.1955787  13.763586 12 1.035056e-08
4 myg_horse C6-vs-C4  3.175654 0.2807086  11.312993 12 9.300920e-08
5    rabbit C6-vs-C4  2.020085 0.2209612   9.142261 12 9.352836e-07
6     yeast C6-vs-C4 -6.366990 0.3134408 -20.313216 12 1.163274e-10
    adj.pvalue issue MissingPercentage ImputationPercentage
1 3.322944e-06    NA        0.40476190           0.40476190
2 1.498854e-10    NA        0.24242424           0.24242424
3 2.070113e-08    NA        0.49479167           0.49479167
4 1.395138e-07    NA        0.22222222           0.22222222
5 1.122340e-06    NA        0.09677419           0.09677419
6 3.489822e-10    NA        0.42222222           0.42222222

火山图只有 6 个点：DDARawData 本来就只有 6 个蛋白（每个 18–32 条肽、共 2070 行特征）。它是受控掺入实验——把已知浓度的蛋白掺进复杂背景，用来检验定量方法本身，不是用来找差异蛋白的发现型数据集。要看蛋白更多的火山图，需要换成真实的大规模定量数据（如 MaxQuant 的 evidence.txt）。

pload(ggplot2)

── loading the package: ggplot2 ──

ggplot(res, aes(log2FC, -log10(adj.pvalue),
                color = adj.pvalue < 0.05 & abs(log2FC) > 1)) +
    geom_point(size = 3) +
    scale_color_manual(values = c("grey80", "firebrick")) +
    theme_minimal() +
    theme(legend.position = "none")

下面用 ivolcano 出交互式版本，鼠标悬停可看数值：

pload(ivolcano)

── loading the package: ivolcano ──

vdat <- data.frame(gene      = res$Protein,
                   logFC     = res$log2FC,
                   adj.P.Val = res$adj.pvalue,
                   pvalue    = res$pvalue,
                   SE        = res$SE)
ivolcano(vdat, logFC_col = "logFC", pval_col = "adj.P.Val", gene_col = "gene",
         title = "C6 vs C4 (MSstats, DDARawData)",
         interactive = TRUE)

实际的质谱数据可先用 MaxQtoMSstatsFormat(evidence, annotation) 将 MaxQuant 的 evidence.txt 转成 MSstats 格式，再走上面 dataProcess() + groupComparison() 流程。

思考：质谱定量的”差异分析”思路与转录组的 DESeq2/limma 有何异同？

五、小结

R 分析链：序列下载（rentrez）→ 基本属性（Biostrings/Peptides）→ 互作网络（mechgraph + ggtangle）→ 功能富集（clusterProfiler）→ 定量（MSstats）。
进一步可拓展：蛋白修饰位点分析、序列 Motif 检索、用 igraph::as.igraph() 把 mechgraph 交给别的 graph 工具进一步分析等。
注意：第 4 节的 DDARawData 每个条件只有 1 个生物学重复（3 个技术重复），MSstats 用技术重复估计误差，这里的 p 值不能当生物学显著性读。