\换行符

简介

ProteoMM是一个同时对单个或多个蛋白质组学数据集进行肽水平差异表达分析的平台。ProteoMM为蛋白质丰度的差异提供了单一的p值和单一的效应量估计。测试统计量是为每个数据集生成的f统计量的和。然后通过排列检验估计p值，因为f统计量的和的分布没有一个封闭形式的解。同时利用多个数据集中的蛋白质中所有可用的多肽，增加了检测条件或处理之间差异的统计能力。此外，在ProteoMM软件包中，我们在之前研究的基础上，提供了规范化、基于模型的缺失肽丰度归一化和肽水平差异蛋白表达分析的功能[1,2]。

目前，多个数据集的组合分析仅限于使用多数据集t-test[3]。由于在蛋白质组学中，蛋白质丰度是根据构成多肽来测量的，因此t检验需要在使用多数据集t检验分析之前，将多肽丰度平均或“汇总”为蛋白质丰度。我们之前已经证明，这种观察数量的减少会导致统计能力的降低和检测差异表达蛋白[1]的能力的降低。ProteoMM为自下而上的基于质谱的蛋白质组学研究提供了从原始肽丰度到蛋白质定量的多个或单个数据集的灵活管道。

本教程将引导读者通过两个模拟数据集的示例分析。函数定义和描述请使用" ?命令。

\换行符

安装

ProteoMM可以从Bioconductor安装:

来源(“//www.andersvercelli.com/biocLite.R”)biocLite(“ProteoMM”)库(ProteoMM)

ProteoMM也可以从GitHub安装:

devtools: install_github(“YuliyaLab / ProteoMM”)库(ProteoMM)

ProteoMM分析管道

ProteoMM Pipeline包括六个步骤，我们建议按照以下顺序执行:

\换行符

负载数据->特征ms归一化->基于模型的Imputation ->基于模型的差分表达式分析和存在/缺席分析->可视化和表输出。

\换行符

个别步骤，如归一化、imputation或存在/不存在分析可以跳过，但必须注意确保传递到基于模型的差异表达分析步骤的肽包含足够数量的观察结果。

我们在本教程中提供的示例遵循图1中建议的ProteoMM分析大纲，并提供了我们认为在蛋白质组学数据分析中有用的额外数据可视化。

\换行符

EigenMS正常化

本例中使用的数据是蛋白质组学实验的一个子集，其中肽id(序列)被打乱，蛋白质和基因id被假的“Prot_#”名称所取代。本文档提供了运行多矩阵分析所需的代码和数据结构示例，包括特征ms归一化、基于模型的imputation和多矩阵统计分析。

对于非蛋白质组学数据，如代谢组学数据，可以提供2列相同的信息。

首先加载数据并定义参数prot.info，两列数据框架，其中代谢物或多肽的id，如果是代谢物，两列是相同的。对于肽，第一列必须包含唯一的肽ID(通常是肽序列)，第二列可以包含蛋白质ID，(在EigenMS中不使用)和任何其他将通过分析管道传播的元数据列。

人类

人类数据集包含695个肽，共13列，其中6列包含强度，其余是描述蛋白质/肽的元数据。CG组和mCG组各有6个样本，各有3个样本。

我们将0替换为NA，并对强度进行log2变换，因为0不应该用于替换缺失的观测值。这种替换会严重扭曲强度分布，产生无效的差异表达结果。更多信息请参见Karpievitch等人，2009[1,2]。

#加载人类数据，然后是鼠标数据(“hs_peptides”)#加载变量hs_peptides dim(hs_peptides) # 695 x 13

## [1] 695 13

intsors = 8:13 #包含强度的列索引m_logInts = make_intents (hs_peptides, intsors) #用NA替换0,NA更适合分析& log2转换m_logInts = convert_log2(m_logInts) metaCols = 1:7 #元数据的列索引，如蛋白质id和序列m_prots .info = make_meta(hs_peptides, metaCols) # m_prots .info - 2+列数据帧与肽id和蛋白质id头(m_prots .info)

# # 1 # #序列MatchedID蛋白质GeneID ProtName CLLAASPENEAGGLKLDGR 3 Prot3 Gene3 Prot3名字# # 2 HNIEGIFTFVDHR 3 Prot3 Gene3 Prot3名字# # 501 Prot501 RLFSGTQISTIAESEDSQESVDSVTDSQKR Gene501 Prot501名字# # 501 LREQYGLGPYEAVTPLTK Prot501 Gene501 Prot501名字# # 502 LINNNPEIFGPLK Prot502 Gene502 Prot502名字# # 6 ENMELEEKEK 14 Prot14 Gene14 Prot14名字# # ProtIDLong GeneIDLong # # 1 Prot3长Gene3 # # 2 Prot3长Gene3 # # 3 Prot501长Gene501 # # 4 Prot501 Gene501长## 5 prote14 long Gene502 long

(m_logInts) # 695 x 6

6 .

GRPS = as。factor(c('CG'，'CG'，'CG'， 'mCG'，'mCG'，'mCG')) # CG和mCG的3个样本#检查缺失值m_nummiss = sum(is.na(m_logInts)) m_nummiss

## [1] 1597

m_numtot = dim(m_logInts)[1] * dim(m_logInts)[2] #观测总数# m_percmiss = m_nummiss/m_numtot # %缺失观测值m_percmiss # 38.29%缺失值，代表真正的较大数据集

## [1] 0.3829736

# plot每个样本中缺失值的数量par(mfcol=c(1,1)) barplot(colsum (is.na(m_logInts))， main="人体样本中缺失值的数量(组顺序)")

注意，mCG治疗组有更多的缺失值。接下来识别偏见趋势eig_norm1 ()．

hs_m_ints_eig1 = eig_norm1(m=m_logInts,treatment=grps,prot.info=m_prot.info)

## eig_norm1(m = m_logInts, treatment = grps, prot.info = m_prot.info)中的警告:该函数使用随机命名器生成器。为了重现性，使用## set.seed(12345)和你选择的参数

##数据维度:6956

##治疗组:CGCGCGmCGmCGmCG

选择完整的肽

获得2+治疗grps

计算SVD，估计特征…

##处理次数:2次

SVD中使用的完整肽(和样品)的数量

# # 1966

##治疗组数(单位:svd.id): 2

##启动引导.....

第50次迭代

##迭代100

##迭代150

第200次迭代

##迭代250

##迭代300

##迭代350

##迭代400

##迭代450

##迭代500

##自动检测到的偏差趋势数量

##准备密谋…

名(hs_m_ints_eig1)

##[1]“m”“treatment”“my”。圣言”“总统”#”[5]n.treatment”“n.u.treatment”“h.c”“礼物“# #”[9]prot.info”“完成”“toplot1”“Tk“# #”[13]ncompl grp”

我们的模拟数据集很小，在没有缺失值的多肽中仅识别出1个偏倚趋势。但从视觉上看，似乎至少有两个。

hs_m_ints_eig1 h.c美元

## [1]

运行EigenMS归一化，消除1个偏差趋势

Hs_m_ints_norm_1bt = eig_norm2(rv=hs_m_ints_eig1)

##独特的治疗组合数量:2

# #正常化……

##处理100肽

##处理200肽

##处理300肽

##处理肽400

##正常化完成!!

正如右上角的百分比所示，趋势解释的百分比方差增加了15%。但是下一个(中间)趋势解释了18%的变化，所以这一趋势的偏差效应可能需要去除。

名(hs_m_ints_eig1)

##[1]“m”“treatment”“my”。圣言”“总统”#”[5]n.treatment”“n.u.treatment”“h.c”“礼物“# #”[9]prot.info”“完成”“toplot1”“Tk“# #”[13]ncompl grp”

数据中有多少肽没有缺失值(完整)?Dim (hs_m_ints_eig1$complete)#偏差趋势识别基于196个肽

## [1] 196 6

我们的模拟数据集很小，只有196个肽没有缺失值，用于识别偏差趋势。只确定了一种偏差趋势，但从视觉上看，似乎至少有两种。这里我们手动将h。c设为2个将被消除的支架。

Hs_m_ints_eig1 $h.c = 2 #视觉上有超过1个偏差趋势，设置为2

Hs_m_ints_norm = eig_norm2(rv=hs_m_ints_eig1)

##独特的治疗组合数量:2

# #正常化……

##处理100肽

##处理200肽

##处理300肽

##处理肽400

##正常化完成!!

图4显示了由趋势解释的28%的方差百分比增加，其中组间的差异解释了数据中71%的总变化，如右上角的百分比所示。下一个(中间)趋势解释了16%的变化，但去除更多趋势的影响可能会过度正态化，因此我们将使用消除两个偏差趋势的规范化数据。

\ newpage

鼠标

小鼠数据集包含1102个肽，共13列，其中6列包含强度，其余为描述蛋白质/肽的元数据。

CG组和mCG组各有6个样本，各有3个样本。数据准备类似于我们对人类数据所做的工作。

Data ("mm_peptides") #加载变量mm_peptides

## [1] 1102 13

Dim (mm_peptides) # 1102 x 13

## [1] 1102 13

头(mm_peptides)

# # 1 # #序列MatchedID蛋白质GeneID ProtName GFAYVQFEDVRDAEDALYNLNRK 64 Prot64 Gene64 Prot64名字# # 61 SKCEELSSLHGQLKEAR Prot61 Gene61 Prot61名字# # 61 Prot61 QDAGSEPVTPASLAALQSDVQPVGHDYVEEVR Gene61 Prot61名字# # 4 TGDQEERQDYINLDESEAAAFDDEWRR 1 Prot1 Gene1 Prot1名字# # 5 IPAYFITVHDPAVPPGEDPDGR 60 Prot60 Gene60 Prot60名字# # 6 GGTPGSGAAAAAGSKPPPSSSASASSSSSSFAQQR 60 Prot60 Gene60 Prot60名字# # ProtIDLong GeneIDLong CG1 CG2 CG3 mCG1 mCG2 mCG3 # # 1 Prot64长Gene64 372590011642000 4872400 0 12850000 3751700 ## 2 Prot61长Gene61长19699000 38055000 30661000 15896000 55187000 20356000 ## 3 Prot61长Gene61长000 5277500 0 38698000 ## 4 Prot1长Gene1长0000 0000 0000 0000 0000 0000 00 5 Prot60长Gene60长9391200 00 4689800 8305300 0 ## 6 Prot60长Gene60长00 20406000 5809800 000

intsles = 8:13 #可能不同于每个数据集，用户需要调整m_logInts = make_intentices (mm_peptides, intsors) #重用名称m_logInts m_logInts = convert_log2(m_logInts) metaCols = 1:7 m_prots .info = make_meta(mm_peptides, metaCols) head(m_prots .info)

# # 1 # #序列MatchedID蛋白质GeneID ProtName GFAYVQFEDVRDAEDALYNLNRK 64 Prot64 Gene64 Prot64名字# # 61 SKCEELSSLHGQLKEAR Prot61 Gene61 Prot61名字# # 61 Prot61 QDAGSEPVTPASLAALQSDVQPVGHDYVEEVR Gene61 Prot61名字# # 4 TGDQEERQDYINLDESEAAAFDDEWRR 1 Prot1 Gene1 Prot1名字# # 5 IPAYFITVHDPAVPPGEDPDGR 60 Prot60 Gene60 Prot60名字# # 6 GGTPGSGAAAAAGSKPPPSSSASASSSSSSFAQQR 60 Prot60 Gene60 Prot60名字# # ProtIDLong GeneIDLong # # 1 Prot64长Gene64 # # 2 Prot61长Gene61 # # 3Prot61长Gene61长## 4 Prot1长Gene1长## 5 Prot60长Gene60长## 6 Prot60长Gene60长

(m_logInts)# 1102 x 6

## [1] 1102

#检查鼠标样本中缺失值的数量m_nummiss = sum(is.na(m_logInts)

## [1] 2698

m_numtot = dim(m_logInts)[1] * dim(m_logInts)[2] #总观测值m_percmiss = m_nummiss/m_numtot # %缺失观测值m_percmiss # 40.8%缺失值，代表真正的较大数据集

## [1] 0.408046

# plot每个样本中缺失值的数量par(mfcol=c(1,1)) barplot(colsum (is.na(m_logInts))， main="鼠标样本中缺失值的数量(组顺序)")

mm_m_ints_eig1 = eig_norm1(m=m_logInts,treatment=grps, prott .info= m_prott .info)

## eig_norm1(m = m_logInts, treatment = grps, prot.info = m_prot.info)中的警告:该函数使用随机命名器生成器。为了重现性，使用## set.seed(12345)和你选择的参数

数据尺寸:11026

##治疗组:CGCGCGmCGmCGmCG

选择完整的肽

获得2+治疗grps

下面的对象从TREATS (pos = 3)中屏蔽:## ## TREATS

计算SVD，估计特征…

##处理次数:2次

SVD中使用的完整肽(和样品)的数量

# # 2706

##治疗组数(单位:svd.id): 2

##启动引导.....

第50次迭代

##迭代100

##迭代150

第200次迭代

##迭代250

##迭代300

##迭代350

##迭代400

##迭代450

##迭代500

##自动检测到的偏差趋势数量

##准备密谋…

解释小鼠数据中大部分变化(45%)的特征趋势不能代表治疗组差异(图5)。原始数据中的第二个趋势仅解释了22%的总变化，类似于治疗组差异，需要归一化。数据的变化以右上角的百分比表示。

mm_m_ints_eig1 h.c美元

## [1]

Mm_m_ints_norm_1bt = eig_norm2(rv=mm_m_ints_eig1)

##独特的治疗组合数量:2

# #正常化……

##处理100肽

##处理200肽

##处理300肽

##处理肽400

##处理500肽

##处理600肽

##处理700肽

##正常化完成!!

解释归一化小鼠数据中大部分变异(43%)的特征是治疗组差异的代表性。原始数据中的第二个趋势只解释了总变化的27%，应该被认为是偏差。

Mm_m_ints_eig1 $h.c = 2

Mm_m_ints_norm = eig_norm2(rv=mm_m_ints_eig1)

##独特的治疗组合数量:2

# #正常化……

##处理100肽

##处理200肽

##处理300肽

##处理肽400

##处理500肽

##处理600肽

##处理700肽

##正常化完成!!

# 190个没有缺失值的pedes用于标识bais趋势($complete)

解释标准化小鼠数据中大多数变异的特征代表了治疗组差异，现在解释了58%的变异。规范化数据中的第二个趋势解释了比图6(24%)中更少的变化，这仍然有点高，但我们将使用这些数据进行分析，以避免过拟合。

length(mm_m_ints_eig1$prot.info$ matchdid) # 1102 -正确

## [1] 1102

length(hs_m_ints_eig1$prot.info$ matchdid) # 695 -可以规范化所有

## [1] 695

length(unique(mm_m_ints_eig1$prot.info$ matchdid)) # 69

69

length(unique(hs_m_ints_eig1$prot.info$ matchdid)) # 69

69

# 787肽被归一化，其余因观察数量低而被消除。

## [1] 787 6

Dim (hs_m_ints_norm$norm_m) # 480肽归一化

## [1] 480 6

\ newpage

\换行符

基于模型的归责

基于模型的imputation使用统计模型来解释肽强度中的信息缺失，并提供无偏倚的，基于模型的，在肽水平[1]进行的蛋白水平差异表达分析。

基于模型的归责方法模拟了两种缺失机制，一种是完全随机缺失，另一种是丰度依赖。完全随机缺失发生在一个样本中没有观察到的肽与它的丰度或任何其他肽的丰度无关。这通常会影响分析中考虑的一小部分多肽。从我们过去的经验来看，它接近5%或所有观测值。丰度依赖缺失是由于左检，其中肽要么不存在，要么存在于仪器检测不到的浓度过低。在这种情况下，我们对肽强度有部分信息，因为我们知道它一定小于其余观察到的肽强度。

人类

我们需要设置元数据和强度以用于imputation。我们将根据ProtID -在矩阵中的位置为蛋白质标识符进行赋值。在这个示例数据集中，ProtID和matchdid可以互换使用。

hs_prots .info = hs_m_ints_norm$normalized[，metaCols] hs_norm_m = hs_m_ints_norm$normalized[， intsols] head(hs_prots .info)

序列匹配ProtID ## HNIEGIFTFVDHR HNIEGIFTFVDHR 3 Prot3 ## RLFSGTQISTIAESEDSQESVDSVTDSQKR RLFSGTQISTIAESEDSQESVDSVTDSQKR 501 Prot501 ## linnpeifgplk LINNNPEIFGPLK 502 Prot502 # ENMELEEKEK ENMELEEKEK 14 Prot14 ## GeneID ProtName ProtIDLong GeneIDLong ## HNIEGIFTFVDHR Gene3 Prot3 Name Prot3 long Gene3 long ##RLFSGTQISTIAESEDSQESVDSVTDSQKR Gene501 Prot501名称Prot501长Gene501长LINNNPEIFGPLK Gene502 Prot502名称Prot502长Gene502长ENMELEEKEK Gene14 Prot14名称Prot14长Gene14长GHEFYNPQKK Gene14 Prot14名称Prot14长Gene14长

头(hs_norm_m)

## CLLAASPENEAGGLKLDGR 24.16344 25.11800 25.39066 24.73530 24.47494 ## HNIEGIFTFVDHR 21.81538 NA 21.42956 21.90027 21.74596 ## RLFSGTQISTIAESEDSQESVDSVTDSQKR 23.52846 22.73723 23.53173 23.03903 23.51463 ## linnpeifgplk NA 22.34531 21.88714 NA 21.09684 ## ENMELEEKEK 27.31511 26.85826 27.39201 27.89371 28.18741 ## ghfynpqkk 24.69609 24.27661 24.96221 24.42590 24.74535 ## # HNIEGIFTFVDHR NA ## # RLFSGTQISTIAESEDSQESVDSVTDSQKR 22.95478## linnnpeifgplk 21.24429 ## enmeleekek 27.83388 ## ghefynpqkk 24.34182

Dim (hs_norm_m) # 480 x 6, raw: 695, 215肽因缺乏而被淘汰

## [1] 480 6

#观察长度(唯一(hs_prot.info$ matchdid)

59

长度(唯一(hs_prot.info$ProtID)) # 59

59

imp_hs = MBimpute(hs_norm_m, grps, prot.info=hs_prot.info, pr_ppos=3, my.pi=0.05, compute_pi=FALSE) #使用默认pi #历史pi=。05代表了%缺失#完全随机缺失的观测值

#检查imputation长度之后的一些数字(唯一的(imp_hs$imp_prot.info$MatchedID)) # 59 - MatchedID id

59

长度(unique(imp_hs$imp_prot.info$ProtID)) # 59 -蛋白质id

59

长度(唯一(imp_hs$imp_prot.info$GeneID)) # 59

59

Dim (imp_hs$imp_prot.info) # 480 x 7估算肽

## [1] 480 7

6 . Dim (imp_hs$y_imputed) # 480 x 6

## [1] 480 6

mylabs = c('CG'，'CG'，'CG'， 'mCG'，'mCG'，'mCG') #与grps相同，这是一个字符串prot_to_plot = 'Prot32' # 43 gene_to_plot = 'Gene32' plot_3_pep_trends_NOfile(as.matrix(hs_m_ints_eig1$m)， hs_ints_eig1 $ prot_info, as.matrix(hs_norm_m)， hs_prot_info, imp_hs$y_imputed, imp_hs$ imp_prot_info, prot_to_plot, 3, gene_to_plot, 4, mylabs)

鼠标

Mm_prot.info = mm_m_ints_norm$normalized[，1:7] mm_norm_m = mm_m_ints_norm$normalized[，8:13] head(Mm_prot.info)

# # # # # #序列GFAYVQFEDVRDAEDALYNLNRK GFAYVQFEDVRDAEDALYNLNRK SKCEELSSLHGQLKEAR SKCEELSSLHGQLKEAR # # IPAYFITVHDPAVPPGEDPDGR IPAYFITVHDPAVPPGEDPDGR # # GGTPGSGAAAAAGSKPPPSSSASASSSSSSFAQQR GGTPGSGAAAAAGSKPPPSSSASASSSSSSFAQQR # # NLGGNYPEK NLGGNYPEK # # ISCAGPQTYKEHLEGQKHK ISCAGPQTYKEHLEGQKHK # # MatchedID蛋白质GeneID ProtName # # GFAYVQFEDVRDAEDALYNLNRK 64 Prot64 Gene64 Prot64名字# # SKCEELSSLHGQLKEAR 61 Prot61 Gene61 Prot61名字# # IPAYFITVHDPAVPPGEDPDGR 60 Prot60 Gene60 Prot60名字# #GGTPGSGAAAAAGSKPPPSSSASASSSSSSFAQQR 60 Prot60 Gene60 Prot60名字# # NLGGNYPEK 28 Prot28 Gene28 Prot28名字# # ISCAGPQTYKEHLEGQKHK 53 Prot53 Gene53 Prot53名字# # ProtIDLong GeneIDLong # # GFAYVQFEDVRDAEDALYNLNRK Prot64长Gene64 # # SKCEELSSLHGQLKEAR Prot61长Gene61 # # IPAYFITVHDPAVPPGEDPDGR Prot60长Gene60 # # GGTPGSGAAAAAGSKPPPSSSASASSSSSSFAQQR Prot60长Gene60 # # NLGGNYPEK Prot28长Gene28 # # ISCAGPQTYKEHLEGQKHK Prot53长Gene53长

头(mm_norm_m)

## GFAYVQFEDVRDAEDALYNLNRK 21.99076 22.78591 22.56999 ## SKCEELSSLHGQLKEAR 24.24259 24.78175 25.25876 22.56999 ## ggtpgsgaaaaagskpppsssasassssfaqqr NA NA 24.28249 22.47006 ## NLGGNYPEK 24.19505 24.89556 24.52888 NA ## ISCAGPQTYKEHLEGQKHK NA 22.50866 23.56617 23.18408 ## mCG2 mCG3 ## GFAYVQFEDVRDAEDALYNLNRK 22.68752 22.83741 ## SKCEELSSLHGQLKEAR 24.76317 24.58578 ## IPAYFITVHDPAVPPGEDPDGR 22.63033 NA ##ggtpgsgaaaaagskpppsssasassssfaqqr na na ## nlggnypek na 24.74703 ## iscagpqtykehlegqkhk na 22.84175

Dim (mm_norm_m) # 787 x 6，原始有:1102肽/行

## [1] 787 6

length(unique(mm_prot.info$ matchdid)) # 56

56

length(unique(mm_prot.info$ProtID)) # 56

56

set.seed(12131) #为重现性设置随机数生成器种子，#否则将获得重复尝试的各种估算值#对于Human，基于ProtID -在蛋白质标识符矩阵中的位置进行impute imp_mm = MBimpute(mm_norm_m, grps, protd .info= mm_prott .info, pr_ppos=3, my.pi=0.05, compute_pi=FALSE)

# #的改动……

# PI =。05通常是一个很好的估计

检查返回的行数是否与规范化数据中的行数相同。

Dim (imp_mm$imp_prot.info) # 787 x 7 -估算肽和787归一化

## [1] 7

6 . Dim (imp_mm$y_imputed) # 787 x

## [1] 787 6

\ newpage

\换行符

基于模型的微分表达式分析

我们将对在小鼠和人类数据集中检测到的蛋白质进行基于模型的组合差异表达分析。对于仅在两个数据集中的一个中识别的蛋白质，将分别对该特定物种进行分析。由于观察数量的增加，多个数据集的组合分析在检测差异表达蛋白方面具有更高的敏感性。

基于模型的组合微分表达式分析

多矩阵分析可推广到2个或更多数据集，因此并行列表用于存储强度、元数据和治疗组信息。第二列元数据帧必须是两个数据集中都存在的蛋白质标识符。在这个模拟数据集中，protid和matchdid将在人和小鼠之间进行匹配，实际上，蛋白质id将有所不同，因为同一种蛋白质，人和小鼠的蛋白质id是不同的。基因id通常只在大写和小写之间存在差异，由于我们未知的原因，一些基因具有不同的id。因此，当比较不同生物之间的蛋白质丰度时，ProtID不是一个在不同生物之间使用的良好标识符，相反，蛋白质id可以基于Ensembl id进行匹配。

首先创建并行列表，将其作为参数传递给微分表达式函数prot_level_multi_part()。首先，将数据分为两个数据集中共有的蛋白质(可以多于2个)和只存在于其中一个或另一个中的蛋白质(唯一的一个或另一个)。这里我们将分析仅在其中一个数据集中观察到的蛋白质，注意“grps”变量对于这里的两个模拟数据集是相同的，但对于用户需要检查的样本的数量和顺序以及grps变量设置为每个数据集的适当因素。还要注意，所有数据集中的治疗组顺序应该相同。不设置组为

控制，控制，控制，治疗，治疗

在一个样本中

治疗治疗治疗控制控制控制

在另一个。

#使并行列表作为参数传递mms = list() treats = list() protinfos = list() mms[[1]] = imp_mm$y_imputed mms[[2]] = imp_hs$y_imputed treats[[1]] = grps treats[[2]] = grps protinfos[[1]] = imp_mm$ imp_protinfo protinfos[[2]] = imp_hs$ imp_protinfo subset_data = subset_proteins(mm_list=mms, prot_info =protinfos， 'MatchedID') names(subset_data)

##[4]“sub_unique_prot.info”“common_list”“sub_unique_mm_list”

mm_dd_only = subset_data$sub_unique_prot.info[[1]] hs_dd_only = subset_data$sub_unique_prot.info[[2]] ugene_mm_dd = unique(mm_dd_only$ matche_did) ugene_hs_dd = unique(hs_dd_only$ matche_did) length(ugene_mm_dd) # 24 -只在鼠标中

## [1]

length(ugene_hs_dd) # 27 -仅限人类使用

27

Nsets = length(mms) nperm = 50 #排列的数量应该是500+发布质量

ptm = proc.time() set.seed=(12357) comb_MBDE = prot_level_multi_part(mm_list=mms, treat=treats, prot_info =protinfos, prot_col_name='ProtID'， nperm=nperm, dataset_suffix=c('MM'， 'HS'))

##计算统计

##进行排列测试

数据集1

数据集2

Proc.time () - PTM #显示运行测试所需的时间

mybreaks = seq(0,1, by=.05) #调整排列测试是通过拉伸值在#区间[0 1]作为理论p值分布par(mfcol=c(1,3)) #总是检查p值#聚在区间[0 .5]hist(comb_MBDE$P_val, breaks=mybreaks, xlab='未调整p值'，main= ") #调整p值看起来不错hist(comb_MBDE$BH_P_val, breaks=mybreaks, xlab='调整p值'，main= ") #聚在区间[0 .5]hist(p。调整(comb_MBDE$P_val, method='BH')， breaks=mybreaks, xlab='BH调整p值'，main= ")

#水平条纹对应于排列测试产生0或#非常小的值，这些重置以提高可视化par(mfcol=c(1,1)) #火山通常看起来更好的大型数据集…plot_volcano_wLab(comb_MBDE$FC, comb_MBDE$BH_P_val, comb_MBDE$GeneID, FC_cutoff=1.2, PV_cutoff=。05、“CG vs mCG”)

仅在人类中观察到的蛋白质的基于模型的差异表达分析

有人类(HS)特定的蛋白质可以用基于模型的差异表达分析进行分析，所以没有对这个子集进行分析。

仅在小鼠中观察到的蛋白质的基于模型的差异表达分析

# subset_data包含“sub_unique_mm_list”“sub_unique_prot.info”列表#为每个数据集，按提供给子集函数mms_mm_dd = subset_data$sub_unique_mm_list[[1]] #鼠标dim(mms_mm_dd) # 258 x 6，

## [1] 258 6

protinfos_mm_dd = subset_data$ sub_unique_prod .info[[1]] length(unique(protinfos_mm_dd$ProtID)) # 24

## [1]

长度(unique(protinfos_mm_dd$GeneID)) # 24

## [1]

长度(unique(protinfos_mm_dd$ matchdid)) # 24

## [1]

DE_mCG_CG_mm_dd = peptideLevel_DE(mms_mm_dd, grps, protinfo =protinfos_mm_dd, pr_ppos=2) #火山图FCval = 1.2 #更改此值为替代折叠变化截断plot_volcano_wLab(DE_mCG_CG_mm_dd$FC, DE_mCG_CG_mm_dd$BH_P_val, DE_mCG_CG_mm_dd$GeneID, FC_cutoff=FCval, PV_cutoff=。05，“鼠标专用- CG vs mCG”)

\ newpage

\换行符

Presence-Absence分析

结合老鼠和人的分析

在存在-缺失分析中，我们只使用不在标准化数据中的蛋白质。例如，由于许多缺失值，对于某些蛋白质，一些肽可能已经被消除，但如果一些肽仍然存在于基于模型的差异表达分析中，我们不会在存在-缺失分析中分析肽的子集，因为我们将获得2个p值。我们坚信，基于模型的差异表达分析是一种更敏感的方法，因此它是在两个治疗组中都有足够数量观察到的蛋白质的首选分析方法。

#让数据结构适合get_presAbs_prots()函数raw_list =列表()norm_imp_prot.info_list = () raw_list [[1]] = mm_m_ints_eig1 $ m raw_list [[2]] = hs_m_ints_eig1 $ m norm_imp_prot.info_list [[1]] = mm_m_ints_eig1美元prot.info norm_imp_prot.info_list [[2]] = hs_m_ints_eig1美元prot.info protnames_norm_list =列表()protnames_norm_list[[1]] =独特(mm_m_ints_norm规范化MatchedID美元)# 56/69 protnames_norm_list[[2]] =独特(hs_m_ints_norm规范化MatchedID美元)# 59 presAbs_dd =get_presAbs_prots(mm_list=raw_list, prot.info=norm_imp_prot.info_list, protnames_norm=protnames_norm_list, prot_col_name=2)

规范化的肽数:1072

肽Pres/Abs数:30

多肽归一化数量:663

肽Pres/Abs数:32

ints_presAbs = list() protmeta_presAbs = list() ints_presAbs[[1]] = presAbs_dd[[1]][[1]] #鼠标ints_presAbs[[2]] = presAbs_dd[[1]][[2]] # HS protmeta_presAbs[[1]] = presAbs_dd[[2]][[1]] protmeta_presAbs[[2]] = presAbs_dd[[2]][[2]] dim(protmeta_presAbs[[2]]) # 32 x 7肽

## [1] 32 7

length(unique(protmeta_presAbs[[2]]$ matchdid)) # 10 -蛋白质

## [1]

dim(protmeta_presAbs[[1]]) # 30 x 7肽

## [1] 30 7

length(unique(protmeta_presAbs[[1]]$ matchdid)) # 13 -蛋白质

## [1]

subset_presAbs = subset_proteins(mm_list=ints_presAbs, prott .info=protmeta_presAbs，' matchdid ') names(subset_presAbs)

##[4]“sub_unique_prot.info”“common_list”“sub_unique_mm_list”

暗(subset_presAbs sub_unique_prot.info美元[[1]])

## [1] 17 7

暗(subset_presAbs sub_unique_prot.info美元[[2]])

## [1] 14 7

暗(subset_presAbs sub_prot.info美元[[1]])

## [1] 13 7

暗(subset_presAbs sub_prot.info美元[[2]])

## [1] 18 7

nperm= 50 #设置为500+用于发布ptm = proc.time() set.seed=(123372) presAbs_comb=prot_level_multiMat_PresAbs(mm_list=subset_presAbs$sub_mm_list, treat=treats, prot_info =subset_presAbs$ sub_prot_info, prot_col_name=' matchdid '， nperm=nperm, dataset_suffix=c('MM'， 'HS'))

在prot_level_multiMat_PresAbs(mm_list = subset_presAbs$sub_mm_list，:该函数使用随机命名器生成器。为了重现性，使用## set.seed(12345)和你选择的参数

##计算统计

##进行排列测试

数据集1

数据集2

Proc.time () - PTM

plot_volcano_wLab(presAbs_comb$FC, presAbs_comb$BH_P_val, presAbs_comb$GeneID, FC_cutoff=.)5, PV_cutoff =。05、“综合Pres/Abs CG vs mCG”)

#只是检查这里的数字dim(subset_presAbs$sub_unique_mm_list[[1]])

## [1] 17 6

暗(subset_presAbs sub_unique_mm_list美元[[2]])

## [1] 14 6

独特的(subset_presAbs sub_unique_prot.info美元[[1]]蛋白质美元)# 8

级别:Prot1 Prot10 Prot11 Prot12 Prot13 Prot14 Prot15 Prot16…Prot9

独特的(subset_presAbs sub_unique_prot.info美元[[2]]蛋白质美元)# 5

级别:Prot1 Prot10 Prot11 Prot12 Prot13 Prot14 Prot15 Prot16…Prot9

仅在小鼠中发现的蛋白质的存在/缺失分析

mm_presAbs = peptideLevel_PresAbsDE(subset_presAbs$sub_unique_mm_list[[1]]， treats[[1]]， subset_presAbs$sub_unique_prot.info[[1]]， pr_ppos=3) plot_volcano_wLab(mm_presAbs$FC, mm_presAbs$BH_P_val, mm_presAbs$GeneID, FC_cutoff=.)5, PV_cutoff =。05、“MM Pres/Abs CG vs mCG”)

仅在人类中发现的蛋白质的存在/缺失分析

hs_presAbs = peptideLevel_PresAbsDE(subset_presAbs$sub_unique_mm_list[[2]]， treats[[2]]， subset_presAbs$sub_unique_prot.info[[2]]， pr_ppos=3) plot_volcano_wLab(hs_presAbs$FC, hs_presAbs$BH_P_val, hs_presAbs$GeneID, FC_cutoff=.)5, PV_cutoff =。05、“HS Pres/Abs CG vs mCG”)

\换行符

参考文献

1. Karpievitch, y.v.，等人，自下而上的基于ms的蛋白质组学中蛋白质定量的统计框架。生物信息学,2009年。25(16): p. 2028-34。

2. Karpievitch, y.v.，等，自下而上的基于ms的蛋白质组学中使用奇异值分解的峰值强度归一化。生物信息学,2009年。25(19): p. 2573-80。

3. Taylor, s.l.，等，多元两部分统计分析相关的质谱数据从多个生物标本。生物信息学,2017年。33(1): p. 17-25。

\换行符

R会话信息

sessionInfo ()

## R版本4.2.0 RC (2022-04-21 r82226) ##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 20.04.4 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas。/home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack。所以## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# # [3] LC_TIME=en_GB LC_COLLATE= c# # [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# # [9] LC_ADDRESS=C lc_phone = c# # [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:## [1]stats graphics grDevices utils datasets methods基础## ##其他附加包:## [1]ProteoMM_1.15.0 ## ##通过命名空间加载(且未附加):[1] Rcpp_1.0.8.3 pillar_1.7.0 compiler_4.2.0 highr_0.9 ## [5] tools_4.2.0 digest_0.6.29 evaluate_0.15 lifecycle_1.0.1 ## [9] tibble_3.1.6 gtable_0.3.0 pkgconfig_2.0.3 rlang_1.0.2 ## [13] cli_3.3.0 DBI_1.1.2 ggrepel_0.9.1 xfun_0.30 ## [17] dplyr_1.0.8 string_1 .4.0 knitr_1.38 generics_0.1.2 ## [21] vctrs_0.4.1 gtools_3.9.2 grid_4.2.0 tidyselect_1.1.2 ## [25] glue_1.6.2 R6_2.5.1 fansi_1.0.3 ggplot2_3.3.5 ## [29] purrr_0.3.4 farver_2.1.0 magrittr_2.0.3 scales_1.2.0 ## [33] ellipsis_0.3.2matrixStats_0.62.0 assertthat_0.2.1 colorspace_2.0-3 ## [37] labeling_0.4.2 utf8_1.2.2 stringi_1.7.6 munsell_0.5.0 ## [41] crayon_1.5.1