Bohrium
robot
新建

空间站广场

论文
Notebooks
比赛
课程
Apps
我的主页
我的Notebooks
我的论文库
我的足迹

我的工作空间

任务
节点
文件
数据集
镜像
项目
数据库
公开
蛋白质组学入门教程(1)数据预处理流程 · 基础蛋白定量
pyOpenms
python
中文
组学
蛋白质组学
OpenMS
pyOpenmspython中文组学蛋白质组学OpenMS
guolj@dp.tech
发布于 2024-04-21
推荐镜像 :openms:0.3
推荐机型 :c2_m16_cpu
赞 2
蛋白质谱定量
读取OpenMS结果文件
数据预处理
数据筛选
零值填充
蛋白质定量
TOP3
iBAQ
TOP3和iBAQ结果对比
RSD/CV 对比
标准品定量对比

©️ Copyright 2024 @ Authors
作者: guolj@dp.tech 📨
日期:2024-04-21
共享协议:本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始:点击上方的 开始连接 按钮,选择 bohrium-notebook:2023-05-31镜像 和任意配置机型即可开始。


本文将以 OpenMS 的结果为基础,主要进行简单的蛋白无标定量方法介绍,包含常用的TOP3和iBAQ。
关于 OpenMS 的流程,你可以查看上一篇 Notebook
如果你已经具备一定的基础,且主要关注于应用侧,你应该考虑开始接触 NextFlow,例如 QuantMS

代码
文本

自己动手,丰衣足食。
——延安生产动员大会,1939年2月


代码
文本

蛋白质谱定量

代码
文本

蛋白质组学定量方法主要可以分为两大类:

  1. 标记定量

    • 标记定量包括化学标记和代谢标记,前者主流方法有ICATiTRAQTMT等,后者主流方法有SILAC等。

    • 在未来,我们可能会结合实测数据对基于这些标记定量方法的实验数据处理进行更深的解释。

  2. 非标记定量

    • 非标记定量方法的一个显著优势是它不需要昂贵的标记物作为内标准,而样本的制备过程也相对简单。

    • 由于涉及较少的操作步骤,可以有效减少人为和系统误差。

    • 此外,通过在样本中加入已知浓度的标准蛋白,这种方法还可以用于蛋白质的绝对定量。

    • 然而,相对于标记定量方法,非标记定量的数据分析更加复杂且难度更大,可重复性通常也较低,对仪器的稳定性和分辨率有较高的要求。

    • 本篇Notebook将对常用的TOP3和iBAQ两个非标记定量方法进行阐述。

代码
文本

读取OpenMS结果文件

我在这里为你准备好了结果文件,它的产生可以见上一篇Notebook

代码
文本
[40]
# 导入库
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib_venn import venn3
代码
文本
[4]
# 导入数据
consensus_df_with_prot = pd.read_csv("/root/proteomics/example_data/consensus_map.csv")
consensus_df_with_prot.head(2)
, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
sequencescoreprotein_accessionprotein_countchargeRTmzqualityZY_UPS2_C3_Mouse_0_6ugZY_UPS2_D2_Mouse_0_6ugZY_UPS2_C2_Mouse_0_6ugZY_UPS2_C1_Mouse_0_6ugZY_UPS2_B3_Mouse_0_6ugZY_UPS2_A1_Mouse_0_6ugZY_UPS2_D3_Mouse_0_6ugZY_UPS2_B1_Mouse_0_6ugZY_UPS2_A2_Mouse_0_6ugZY_UPS2_A3_Mouse_0_6ugZY_UPS2_B2_Mouse_0_6ugZY_UPS2_D1_Mouse_0_6ug
0IDLAVGDVVK11.477765Q8K1B8|URP2_MOUSE1.023254.861012514.8060140.915005269256000.0417071300.0405560700.0398971200.0285516400.0244054700.0360338700.0361681600.0260860600.0240820200.0319740100.0465238500.0
1VMSQNFTNC(Carbamidomethyl)HTK15.736818P97855|G3BP1_MOUSE1.021492.125672733.8263880.961677166461400.0163396900.0182838000.0181990100.0124401600.0170870100.0170593700.0158083400.0148514500.078017240.0142878300.0153835500.0
,
代码
文本

上面表格中每个表头的含义如下:

  1. sequence:肽段的氨基酸序列,包含PTM注释。例子中的VMSQNFTNC(Carbamidomethyl)HTK的序列为VMSQNFTNCHTK,其中第九号位上的C为Carbamidomethyl修饰。

  2. score:搜索分数,表示肽段鉴定的置信度。分数越高,鉴定的可信度越高。目前OpenMS并没有给出一个推荐的阈值,但你还是可以查看官方文档了解更多。

  3. protein_accession:蛋白质访问号,通常指的是蛋白质数据库中的一个唯一标识符,用于标记特定的蛋白质。这与我们当初给定的搜索数据库中的蛋白质Description相关。

  4. protein_count:表示鉴定到的该序列肽段所对应的蛋白质数量。通常用于反映一个肽段可能属于多个蛋白质的情况。与protein_accession中的数量一致。

    • 你可以利用consensus_df_with_prot.query("protein_count > 1")查看protein_count > 1的情况。
  5. charge:肽段的电荷状态。在质谱分析中,肽段离子的电荷影响其质量/电荷比(m/z),进而影响检测和分析。你可以查看这里获得更多相关信息。

  6. RT:肽段的LC部分洗脱峰最大值对应的时间点。

  7. mz:质量/电荷比,是质谱分析中的一个基本参数,用于描述肽段或蛋白质离子的质量与其电荷的比值。

  8. quality:质量分数,通常用来表示肽段信号的质量,包括信噪比、峰的对称性等因素。

  9. ZY_***:样品名,其数值即为该样品中对应肽段的定量结果(通常为面积)

代码
文本

数据预处理

代码
文本

数据筛选

代码
文本

接下来,我们看一下quality和score的数值分布情况

代码
文本
[23]
fig, axes = plt.subplots(1, 2, figsize=(6, 1), sharey=True)
consensus_df_with_prot["score"].plot.hist(bins=50, ax=axes[0])
axes[0].set_title("score")
consensus_df_with_prot["quality"].plot.hist(bins=50, ax=axes[1])
axes[1].set_title("quality")
Text(0.5, 1.0, 'quality')
代码
文本

一般地,其他软件会提供E-Value供筛选时需要,但OpenMS目前还未实现这个功能

我们将拍个脑袋,根据quality和score对数据进行过滤,实际项目进行中,我们可以根据内标蛋白/QC样品的RSD(Relative Standard Deviation,或称CV,Coefficient of Variation)或经验来决定阈值。

另外,当一条多肽与多个蛋白质相匹配时,这种情形被称作“多肽共享”或“多肽歧义”。

在利用质谱技术进行蛋白质的非标记定量分析时,这会带来特定的挑战,主要是因为无法确切知道这些多肽具体属于哪个蛋白质。

对于这种情况的处理方法主要有以下几种:

  1. 最保守分配:只有当一个多肽唯一匹配到一个蛋白质时,才将其用于定量。

  2. 分配比例:这种方法根据特定的规则将多肽的信号按照一定比例分配给各个匹配的蛋白质。比如,可以依据每个蛋白质的独有多肽丰度来计算共享多肽的贡献量。

  3. 最小化共享多肽的使用:在定量分析中优先使用独特多肽,只在必要时使用共享多肽。这种策略旨在最大程度上减少因共享多肽带来的计算复杂性和潜在误差。

  4. 统计模型:使用复杂的统计模型来估计每个蛋白质的丰度,这些模型可以考虑多肽共享的影响。例如,软件如MaxQuant中的LFQ (Label-Free Quantification) 算法,可以处理多肽共享问题,通过算法内部的优化处理来估计蛋白质的丰度。

作为入门教学,在本篇Notebook中,我们将采用最保守分配法进行数据预处理。

代码
文本
[46]
# 按quality > 0.8、score > 10和protein_count = 1过滤数据
min_quality = 0.8
min_score = 10
filtered_df = consensus_df_with_prot.query("quality > @min_quality and score > @min_score and protein_count == 1")
filtered_df.shape, consensus_df_with_prot.shape
((6836, 20), (33161, 20))
代码
文本

可以看到,这里筛选的过程使得,我们丢掉了几乎80%的数据,他们的Venn图如下

代码
文本
[47]
meet_quality = set(consensus_df_with_prot.query("quality > @min_quality").index)
meet_score = set(consensus_df_with_prot.query("score > @min_score").index)
meet_protein_count = set(consensus_df_with_prot.query("protein_count == 1").index)
venn3([meet_quality, meet_score, meet_protein_count], ["quality", "score", "protein_count"])
<matplotlib_venn._common.VennDiagram at 0x7fed763a04f0>
代码
文本

另外,更加严格地,我们只保留鉴定出的独特肽段 > 1的蛋白质

代码
文本
[69]
protein_peptide_count = filtered_df["protein_accession"].value_counts()
meet_protein_peptide_count = protein_peptide_count[protein_peptide_count > 1].index
filtered_df_more_than_one_peptide = filtered_df[filtered_df["protein_accession"].isin(meet_protein_peptide_count)]
filtered_df_more_than_one_peptide
, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
sequencescoreprotein_accessionprotein_countchargeRTmzqualityZY_UPS2_C3_Mouse_0_6ugZY_UPS2_D2_Mouse_0_6ugZY_UPS2_C2_Mouse_0_6ugZY_UPS2_C1_Mouse_0_6ugZY_UPS2_B3_Mouse_0_6ugZY_UPS2_A1_Mouse_0_6ugZY_UPS2_D3_Mouse_0_6ugZY_UPS2_B1_Mouse_0_6ugZY_UPS2_A2_Mouse_0_6ugZY_UPS2_A3_Mouse_0_6ugZY_UPS2_B2_Mouse_0_6ugZY_UPS2_D1_Mouse_0_6ug
0IDLAVGDVVK11.477765Q8K1B8|URP2_MOUSE1.023254.861012514.8060140.915005269256000.0417071300.0405560700.0398971200.0285516400.0244054700.0360338700.0361681600.0260860600.0240820200.0319740100.0465238500.0
1VMSQNFTNC(Carbamidomethyl)HTK15.736818P97855|G3BP1_MOUSE1.021492.125672733.8263880.961677166461400.0163396900.0182838000.0181990100.0124401600.0170870100.0170593700.0158083400.0148514500.078017240.0142878300.0153835500.0
2FVDEEDGGDGQAGPDEGEVDSC(Carbamidomethyl)LR40.109421Q9CPW4|ARPC5_MOUSE1.022891.4788781277.0211520.96023561894790.047759900.064873340.061048470.060498050.041455210.048406750.063682410.043125160.048718520.060900750.056396670.0
3EKPQALVTSPATPLPAGSGIK38.518372O88271|CFDP1_MOUSE1.032909.676320688.0570990.94373782818720.069078570.080557970.080324010.064225850.064915200.074057170.074447260.055183950.059126850.068997690.074162200.0
5LDLEAWFPGSGAFR18.670906P26638|SYSC_MOUSE1.024483.495941783.3935570.97626543836660.039589420.040812830.035405420.033558470.045752620.035412970.046204930.042939580.037579030.037180490.052378800.0
...............................................................
14238SSEMNVLIPTEGGDFNEFPVPEQFK26.028290P40124|CAP1_MOUSE1.034354.917646937.7780260.802763168539700.0159314500.0190462600.0167633400.0165573700.097061050.0157201000.0188045500.091055600.00.0195871800.00.0
14247LAGVTALSC(Carbamidomethyl)WLPLR21.785582P97823|LYPA1_MOUSE1.024460.107468778.9365590.80735859122810.047262840.053875060.045753870.037983950.040344040.044332540.048201210.044635930.00.047577000.00.0
14248MFIGGLSWDTTK19.285593Q60668|HNRPD_MOUSE1.024071.956379678.3379350.802992216008800.0226267900.0210280400.0227036000.0191345200.0175056200.00.0221198200.0175835800.0165210800.0207849600.00.0
14335QAFTDVATGSLGQGLGAAC(Carbamidomethyl)GMAYTGK48.753551P40142|TKT_MOUSE1.024007.7535991266.5965950.800172129932500.0122962900.0118945200.00.0108757800.0130031000.0113474400.0117755100.0131244700.0135655600.0121221200.00.0
14338ITGEAFVQFASQELAEK18.519884Q9Z2X1|HNRPF_MOUSE1.034331.446865623.3196500.80739090952630.00.084579110.044385430.045189660.061535430.044376070.0108048800.026259210.034415850.043703790.00.0
,

6205 rows × 20 columns

,
代码
文本

零值填充

对于检出为0的结果,我们需要对其数值进行填充。

这里我们采用Scikit-Learn中的KNNImputer

代码
文本
[150]
from sklearn.impute import KNNImputer

# 丢掉不需要的列
# 当然,你也可以用filtered_df,这会包含独特肽段仅为1的蛋白质
data = filtered_df_more_than_one_peptide.drop(columns=["score", "protein_count", "charge", "RT", "mz", "quality"]).reset_index(drop=True)

imputer = KNNImputer(
missing_values=0,
n_neighbors=5,
weights="uniform",
metric="nan_euclidean",
keep_empty_features=True,
)

imputed_data = imputer.fit_transform(data.iloc[:, 2:])
imputed_data = pd.concat([data.iloc[:, :2], pd.DataFrame(imputed_data, columns=data.columns[2:])], axis=1)

# 由于包含多个组,我们先将数据“融化”开
imputed_data = imputed_data.melt(id_vars=["protein_accession", "sequence"], var_name="sample", value_name="intensity")
代码
文本

蛋白质定量

代码
文本

由于蛋白质在蛋白酶的催化下降解本质上是一个化学过程,这个反应如下:

在我们鉴定出的结果中,一般而言,(当然,你可以自行再验证一下)

于是,直觉上,蛋白质的浓度应当于其消化产物中肽段的浓度基本一致(虽然实际上不同肽段会因为性质差异表现出质谱检出强度的偏差,但在入门教学中我们暂时不引入过多考虑)

代码
文本

TOP3

接下来,我们采用TOP3法对蛋白质进行定量。

上述提到,不同肽段会因为性质差异,表现出质谱检出强度的偏差,即:

而对于检出强度较高的肽段,这个偏差对强度本身的影响较小

代码
文本
[161]
def top3(in_df: pd.DataFrame):
return in_df.sort_values("intensity", ascending=False)["intensity"].head(3).sum()

top3_result = imputed_data.groupby(["sample", "protein_accession"]).apply(top3).reset_index().rename(columns={0: "top3"})
代码
文本

简单运行上面代码,我们便获得了TOP3法的蛋白非标记定量结果,之后,我们运行完iBAQ之后再对比二者的结果

代码
文本

iBAQ

接下来,我们采用iBAQ法对蛋白质进行定量。

iBAQ(intensity-Based Absolute Quantification)是一种用于蛋白质组学定量的技术,主要用于估算蛋白质在样本中的绝对丰度。

iBAQ的原理基于蛋白质的肽段信号强度总和来估计蛋白质的绝对量。

其包括下列步骤:

  1. 对一个蛋白质的所有鉴定到的肽段的信号强度进行求和,得到该蛋白质的总信号强度。

  2. 将这个总信号强度除以蛋白质中可观测到的肽段数量(理论上可以被鉴定的肽段数量),得到iBAQ值。

  3. iBAQ值 = 总信号强度 / 可观测肽段数量

通过比较不同蛋白质的iBAQ值,可以估算它们在样本中的相对丰度。如果有标准蛋白或已知浓度的蛋白质参与实验,还可以通过iBAQ值估算蛋白质的绝对浓度。

代码
文本
[175]
from pyopenms import AASequence, ProteaseDigestion
from Bio import SeqIO
from tqdm import tqdm

# 读入搜库时使用的数据库文件
database = SeqIO.parse("/root/proteomics/sp_mouse_ups2.fasta", "fasta")
acc_to_seq = {i.id: str(i.seq) for i in database}

# 利用OpenMS的方法获取蛋白质理论上可鉴定出的肽段数量
def unique_digest_products_num(
sequence: str,
enzyme: str = "Trypsin",
min_length: int = 6,
max_length: int = 40,
allow_missed_cleavage: int = 2
):
digester = ProteaseDigestion()
digester.setEnzyme(enzyme)
digester.setMissedCleavages(allow_missed_cleavage)
# 下面我们使用了一个黑魔法:=,这是Python 3.8引入的海象运算符
# 因为它很可爱,所以我们在这里用了它,并进行了强调
# 但是需要注意,它会降低代码的可读性,所以在实际工作中,你可能不会经常看到它
digester.digest(AASequence.fromString(sequence), digest_products:= [], min_length, max_length)
return len(set(digest_products))

acc_to_product_num = {
acc: unique_digest_products_num(acc_to_seq[acc]) for acc in tqdm(imputed_data["protein_accession"].unique())
}

def ibaq(in_df: pd.DataFrame):
return in_df["intensity"].sum() / acc_to_product_num[in_df["protein_accession"].iloc[0]]

ibaq_result = imputed_data.groupby(["sample", "protein_accession"]).apply(ibaq).reset_index().rename(columns={0: "ibaq"})
100%|██████████| 1086/1086 [00:00<00:00, 2485.33it/s]
代码
文本

运行上面的代码,我们便获得了iBAQ的定量结果,接下来,我们简短地欣赏一下海象后,来对比一下两个方法的定量结果

海象

:=海象运算符

顺便推荐一个Python黑魔法手册,虽然在工程代码中我会希望你少用这些怪东西

代码
文本

TOP3和iBAQ结果对比

代码
文本
[183]
import numpy as np

results_df = top3_result.merge(ibaq_result, on=["sample", "protein_accession"])
# group为组别,repeat为重复实验的编号
results_df[["group", "repeat"]] = np.vstack(results_df["sample"].apply(lambda x: tuple(x.split("_")[2])).values)
# 当数据过多的时候(比如上亿),你可以考虑提前生成一个字典,然后用map
# sample_group_repeat = {
# i: tuple(i.split("_")[2]) for i in results_df["sample"].unique()
# }
# results_df[["group", "repeat"]] = np.vstack(results_df["sample"].map(sample_group_repeat).values)
代码
文本

RSD/CV 对比

代码
文本
[190]
def calculate_rsd(in_df: pd.DataFrame):
v = in_df[["top3", "ibaq"]]
return v.std() / v.mean()

rsd_df = results_df.groupby(["group", "protein_accession"]).apply(calculate_rsd)
代码
文本
[ ]
fig, ax = plt.subplots()
rsd_df["top3"].plot.hist(bins=50, ax=ax, alpha=0.5, label="top3")
rsd_df["ibaq"].plot.hist(bins=50, ax=ax, alpha=0.5, label="ibaq")
plt.legend()
<matplotlib.legend.Legend at 0x7fecbaa5ea10>
代码
文本
[206]
rsd_df.mean()
top3    0.121571
,ibaq    0.118506
,dtype: float64
代码
文本

可以看到,iBAQ的RSD略小于TOP3(越小越好),在更复杂的样品和更粗的筛选条件中,iBAQ的领先一般会更多。

你可以调整数据筛选阶段的阈值来试试

代码
文本

标准品定量对比

实验在样品中加入了UPS2,这是一个包含了一些的已知相对浓度蛋白质的标准品,常用来作为样品间绝对定量的内标。

在这里,我们可以根据它们的结果评估定量方法。

代码
文本
[234]
ups_amount = pd.read_csv("/root/proteomics/StandardProteins.txt", sep="\t", index_col=0)["Amount"].to_dict()

resuls_df_ups = results_df.query("protein_accession.str.lower().str.endswith('_ups')").copy()
resuls_df_ups["real_amount"] = resuls_df_ups["protein_accession"].apply(lambda x: ups_amount[x.split("|")[0]])
代码
文本

我们一般会采用非线性拟合(特定的高线性区段也可以采用线性拟合)根据内标的定量结果来定量其他蛋白

因此我们也可以简单地观察内标定量结果与真实值的Spearman相关系数来考察不同方法定量的可靠性

代码
文本
[246]
def get_corr(in_df: pd.DataFrame):
agg_df = in_df.groupby("protein_accession")[["top3", "ibaq", "real_amount"]].agg("mean")
# 当正态性得到验证后(一般没问题),你也可以将spearman换成pearson,从而考察线性相关性
return agg_df.corr("spearman").loc["real_amount", ["top3", "ibaq"]]

resuls_df_ups.groupby("group").apply(get_corr)
, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
real_amounttop3ibaq
group
A0.7004640.700464
B0.8731810.850792
C0.7372460.794819
D0.5357430.698865
,
代码
文本

其中A - B - C - D四个group为逐步稀释后的样品

可以看到,在样品浓度较低的情况下,iBAQ方法显著优于top3

代码
文本

本Notebook对常用的TOP3和iBAQ进行了代码实现和阐述。随着硬件和算法的进步,非标记定量领域也一直在发展,一些新的算法也被提出,例如LFAQ、FlashLFQ等。

未来,我们或许会对各个算法进行复现和解读。 :sweat_drops:

另外,根据情况,或许也会写一些Notebook来展示MaxQuant等工具的使用。

当然,这非常取决于作者的心情,所以,你不妨通过点赞的形式催更。 :stuck_out_tongue_winking_eye:

代码
文本
pyOpenms
python
中文
组学
蛋白质组学
OpenMS
pyOpenmspython中文组学蛋白质组学OpenMS
已赞2
本文被以下合集收录
蛋白质组学系列
guolj@dp.tech
更新于 2024-05-06
2 篇1 人关注
推荐阅读
公开
PyTorch Lightning官网教程搬运整合-02:设置最佳保存点和早停(Early Stopping)
pythonPyTorch LightningTutorial中文Deep LearningPyTorch
pythonPyTorch LightningTutorial中文Deep LearningPyTorch
guolj@dp.tech
发布于 2023-11-08
3 赞1 转存文件
公开
PyTorch Lightning官网教程搬运整合-01:训练模型初体验
pythonDeep LearningPyTorchPyTorch LightningTutorial
pythonDeep LearningPyTorchPyTorch LightningTutorial
guolj@dp.tech
发布于 2023-09-24
3 赞
评论
 由于蛋白质在蛋白酶的催化下降解本质上是一...

Suice。

04-24 05:05
n = a = b = c = ...这是不是写错了

guolj@dp.tech

作者
04-24 20:20
回复 Suice。 一般而言,假设所有检出肽段都是独特的话(即在蛋白全长中无重复),它们确实是相等的

Suice。

04-25 20:17
回复 guolj@dp.tech 了解了~谢谢
评论