蛋白质组学数据分析实践（一）

2021-04-22 本文已影响0人克里克的钟

蛋白质组学分析实践（一）

文章标题：The Primary Effect on the Proteome of ARID1A-mutated Ovarian Clear Cell Carcinoma is Downregulation of the Mevalonate Pathway at the Post-transcriptional Level
期刊：Molecular & cellular proteomics
年份 : 2016
DOI: 10.1074/mcp.M116.062539.
数据来源：http://proteomecentral.proteomexchange.org， ID: PXD004570.

说明：本文是学习了生信技能树公众号蛋白质组学相关推文学习整理的。

MaxQuant搜库

MaxQuant搜库的部分主要是电脑在跑程序，设置几个参数，注意experiment中的设置，应该为每一个重复试验文件为一个名字，相互之间各不相同，才能得到6个LFQ值进行统计分析。（第一次搜库就是没有设置这个，导致只有一个LFQ的值）。

顺便说一下，MaxQuant搜库，OVCA429细胞敲除ARID1A和对照的这六个样本数据，电脑配置很菜：锐龙2600+8G DDR5 3000，跑了一天的时间才结束（设置的使用核心数是6个，电脑卡的几乎不能做其他事情）。

后面的分析直接使用了搜库结果的txt数据。

Perseus部分

导入数据

导入combined文件夹 --> txt文件夹 --> ProteinGroups.txt，选择LFQ intensity [分组]的6列数据为Main数据，其他的基本自动填充的。

数据筛选质控

Filter rows --> Filter rows based ong category column，里面有3项，分3次除去之后，得到5681个蛋白质；然后剔除只匹配到一个肽段的蛋白（single peptide hits），这里使用的是Filter rows --> Filter rows based ong numeric/main column，选择Razer + unique peptides，Relation 1 写入x>1， OK。

Filter rows

得到4973行，也就是4973个蛋白质，与文章中的描述一致。

4973个蛋白质

聚类分析

文章结果的第二节（Label-free Comparisons Between ARID1A Knockout and Control Proteomes），就讲到使用的是LFQ来计算ARID1A敲除蛋白组和对照蛋白组中蛋白质组的相对丰度。

归一化和缺失值的处理

首先，点击Annot.Row，将样本分为两组，然后再进行数据转化操作。

Annot.Row

然后点击Basci --> transform，转换为log2的数据。这里会引入一些空值（NaN），因为数据本身有很多0，也就是没有匹配到蛋白信息。这些可能跟转录组有些差异，不能直接log(2+1)，因为这样会导致数据偏差太大。文中的描述是：缺失值被假定为偏向于低于质谱检测限的低丰度蛋白质，称为：“missing not at random";(这是蛋白质组学研究中经常作出的假定); 缺失值被替换为中位数下移高斯分布中的随机值，以模拟低丰度LFQ值；每个样本分别从宽度为0.3,downshift为1.8的分布中进行估算。第二步，进行缺失值的处理，在imputation --> Replace missing values from normal distribution，默认参数，确定。

imputation

第三步，normalization，选择normalization --> Z-score，确定。2.3 聚类分析选择归一化后的数据，点击clustering图标，稍经调整，就可以得到文章中类似的图了

HeatMap

原文的聚类图

可以看到，趋势基本上是一致的。

差异分析-火山图

火山图选择log2后填充缺失值的数据，然后点击火山图的图标。使用的是双侧t检验，FDR在这里为T检验的p值，而S0是方差，当S0设置为0时表示仅p值起作用。number of randomization不知道是什么意思。结果如图：

ScatterPlot

目前只能达到这个效果，不知道设置。t检验显著的蛋白有2896个（奇怪的是每次计算都有所差异，是随机化的问题吗？）【用R做t.test的差异蛋白是2613，这其中的差别在哪里？】，然后计算|log2FC| > 1的有422个。文献讲到的是430和2606个，有些差别。在图的调整上，自由度也是非常有限。所以，如果要画火山图，还是将数据导出来，在R中绘制比较好。

通过R来绘制火山图：

火山图jpeg

原文的图：

原文的火山图

是相差不多的。火山图的代码：

draw_volcano_plot <- function(need_DEG,logFC_cutoff){
  if(! logFC_cutoff){
    logFC_cutoff <- with(need_DEG,mean(abs(log2FoldChange)) + 2*sd(abs( log2FoldChange)) )
  }
  #logFC_cutoff=1
  
  need_DEG$change = as.factor(ifelse(need_DEG$Pvalue < 0.05 & abs(need_DEG$log2FoldChange) > logFC_cutoff,
                                     ifelse(need_DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT')
  )
  
  library(ggplot2)
  g = ggplot(data=need_DEG,
             aes(x=log2FoldChange, y=-log10(Pvalue))) +
    geom_point(aes(color=change)) +
    xlab("log2 fold-change") + ylab("-log(p-value)") +
    scale_x_continuous(limits = c(-10, 10))+
    scale_y_continuous(limits = c(0, 8))+
    scale_colour_manual(values = c("#00AAAA",'darkgray','#00AAAA')) + ## corresponding to the levels(res$change)
    geom_hline(yintercept = -log10(0.05), linetype =2,size = 1, color = 'red')+
    geom_vline(xintercept = logFC_cutoff, linetype =2,size = 1, color = 'red')+
    geom_vline(xintercept = -logFC_cutoff, linetype =2, size =1,color = 'red')+
    geom_text(aes(x = -6,y = 7),label = "179 Proteins\n downregulated", color = "black", size = 4)+
    geom_text(aes(x = 6,y = 7),label = "91 Proteins\n upregulated",color = "black", size = 4)+
    theme_prism(border = TRUE) +
    coord_cartesian(clip = "off")
  print(g)
}
draw_volcano_plot(dat, 1)

原文的上调蛋白是95，下调169。本次分析上调91，下调179。【可能有存在差异的地方】但是，后面的倍数变化的表格，也是对得上的。

倍数变化的表格

如果数据分析仅仅是这样的话，完全可以通过R来操作，归一化，缺失值的处理，剔除不符合条件的列，这些都可以通过R来分析。热图也没有问题。可能Perseus还有一些其他重要的内容吧，有时间把B站的视频学一遍。

Bio Zhong
本人博客：www.clockstudy.com