这篇Nature子刊文章的蛋白组学数据PCA分析竟花费了我两天时间来重现|附全过程代码

2020年4月14日，Sanger研究团队于nature communication在线发表了题为Single-cell transcriptomics identifies an effectorness gradient shaping the response of CD4+ T cells to cytokines的研究内容，作者使用蛋白质组学、bulk RNA-seq和单细胞转录组测序对人体40,000个以上的naïve and memory CD4+ T cells进行分析，发现细胞类型之间的细胞因子反应差异很大。memory T细胞不能分化为Th2表型，但可以响应iTreg极化获得类似Th17的表型。单细胞分析表明，T细胞构成了一个转录连续体（transcriptional continuum），从幼稚到中枢和效应记忆T细胞，形成了一种效应梯度，并伴随着趋化因子和细胞因子表达的增加。最后，作者表明，T细胞活化和细胞因子反应受效应梯度的影响。

该文献通过蛋白质组学（(液相色谱-串联质谱法，LC-MS/MS）进行了探索性分析，样品对应于从健康个体的外周血中分离的幼稚和记忆T细胞，并用多种细胞因子刺激5天，每个条件平均3个生物学重复。

这次复现Fig1cPCA图和Fig2aPCA图的另一部分，这次作者是通过蛋白组学数据进行PCA的展现：

以上是Fig1c原图，图注为“PCA plots from the whole transcriptome of TN and TM cells. Different colors correspond to cell types and different shades to stimulation time points. PCA plots were derived using 21 naive and 19 memory T cell samples for proteomics”

以上为Fig 2a原图，图注为“PCA plot from the full transcriptome of TN and TM cells following five days of cytokine stimulations. Only stimulated cells were included in this analysis. PCA plots were derived using 18 naive and 17 memory T cells samples ”

我们需要复现该图之前，先需要下载数据，可以点击https://www.opentargets.org/projects/effectorness对proteomics的abundances数据和metadata数据进行下载，然后进行以下步骤：

library(SummarizedExperiment)
library(annotables)
library(rafalib)
library(ggplot2)
library(ggrepel)library(limma)

加载数据

加载标准化后的丰度：

MassSpec_data <- read.table("NCOMMS-19-7936188_MassSpec_scaled_abundances.txt", header = T, stringsAsFactors = F)
View(MassSpec_data)#从以上可以看出，每列除了代表每个样本外，前三列分别为Protein_id，Gene_id和Gene_name，每行代表一个蛋白

建立SummarizedExperiment object

创建带有蛋白质注释的dataframe

protein_annotations <- data.frame(MassSpec_data[,c("Protein_id","Gene_id","Gene_name")], row.names = MassSpec_data$Gene_name)
rownames(MassSpec_data) <- MassSpec_data$Gene_name#构成一个由"Protein_id","Gene_id","Gene_name"的数据框MassSpec_data <- MassSpec_data[,-c(1:3)]

创建带有sample注释的dataframe

sample_ids <- colnames(MassSpec_data)
sample_annotations <- data.frame(row.names = sample_ids,
donor_id = sapply(sample_ids, function(x){strsplit(x, split = "_")[[1]][1]}),
cell_type = paste("CD4",
sapply(sample_ids, function(x){strsplit(x, split = "_")[[1]][3]}),
sep="_"),
cytokine_condition = sapply(sample_ids, function(x){strsplit(x, split = "_")[[1]][4]}),
stringsAsFactors = T)
sample_annotations$activation_status <- ifelse(sample_annotations$cytokine_condition == "resting", "Resting", "Activated")View(sample_annotations)

创建relevant metadata的变量

meta <- list(
Study="Mapping cytokine induced gene expression changes in human CD4+ T cells",
Experiment="Quantitative proteomics (LC-MS/MS) panel of cytokine induced T cell polarisations",
Laboratory="Trynka Group, Wellcome Sanger Institute",
Experimenter=c("Eddie Cano-Gamez",
"Blagoje Soskic",
"Deborah Plowman"),
Description="To study cytokine-induced cell polarisation, we isolated human naive and memory CD4+ T cells in triplicate from peripheral blood of healthy individuals. Next, we polarised the cells with different cytokine combinations linked to autoimmunity and performed LC-MS/MS.",
Methdology="LC-MS/MS with isobaric labelling",
Characteristics="Data type: Normalised, scaled protein abundances",
Date="September, 2019",
URL="https://doi.org/10.1101/753731")

建立SummarizedExperiment object

proteomics_data <- SummarizedExperiment(assays=list(counts=as.matrix(MassSpec_data)),
colData=sample_annotations,
rowData=protein_annotations,
metadata=meta)saveRDS(proteomics_data, file="proteinAbundances_summarizedExperiment.rds")

数据可视化

将NA值设置为零
注意：此操作仅出于可视化目的。执行统计测试时，NA不会设置为零。

assay(proteomics_data)[is.na(assay(proteomics_data))] <- 0

定义函数：

\bullet 提取蛋白质表达值;

\bullet 进行主成分分析;

\bullet 返回一个矩阵，其中包含每个样品和样品注释的PC坐标;

\bullet 返回每个主要成分解释的方差百分比。

getPCs <- function(exp){
pcs <- prcomp(t(assay(exp)))
pVar <- pcs$sdev^2/sum(pcs$sdev^2)
pca.mat <- data.frame(pcs$x)
pca.mat$donor_id <- colData(exp)$donor_id
pca.mat$cell_type <- colData(exp)$cell_type
pca.mat$cytokine_condition <- colData(exp)$cytokine_condition
pca.mat$activation_status <- colData(exp)$activation_status
res <- list(pcs = pca.mat, pVar=pVar)
return(res)}

对所有样本执行PCA

pcs <- getPCs(proteomics_data)
ggplot(data=pcs$pcs, aes(x=PC1, y=PC2, color=cell_type, shape=activation_status)) +
geom_point(size = 8) +
xlab(paste0("PC1:", round(pcs$pVar[1]*100), "% variance")) +
ylab(paste0("PC2: ", round(pcs$pVar[2]*100), "% variance")) +
scale_colour_manual(values = c("#5AB4AC","#AF8DC3")) +
scale_alpha_discrete(range = c(0.5,1)) +
coord_fixed() + theme_bw() +theme(panel.grid = element_blank())

去掉个体间变异性：

proteomics_data_regressed <- proteomics_data
assay(proteomics_data_regressed) <- removeBatchEffect(assay(proteomics_data_regressed),
batch = factor(as.vector(colData(proteomics_data_regressed)$donor_id)))

重新计算PCA：

pcs <- getPCs(proteomics_data_regressed)
ggplot(data=pcs$pcs, aes(x=PC1, y=PC2, color=cell_type, shape=activation_status)) +
geom_point(size = 8) +
xlab(paste0("PC1:", round(pcs$pVar[1]*100), "% variance")) +
ylab(paste0("PC2: ", round(pcs$pVar[2]*100), "% variance")) +
scale_colour_manual(values = c("#5AB4AC","#AF8DC3")) +
scale_alpha_discrete(range = c(0.5,1)) +
coord_fixed() + theme_bw() +theme(panel.grid = element_blank())

原图

细胞类型特异性分析

将naive和memory T细胞样本分为仅包含受刺激细胞的两个不同数据集。

proteomics_data_naive <- proteomics_data[,(proteomics_data$cell_type=="CD4_naive") & (proteomics_data$activation_status=="Activated")]proteomics_data_memory <- proteomics_data[,(proteomics_data$cell_type=="CD4_memory") & (proteomics_data$activation_status=="Activated")]

Naive T cells

对 5 days-stimulated naive T cells进行PCA:

pcs_naive <- getPCs(proteomics_data_naive)
ggplot(data=pcs_naive$pcs, aes(x=PC1, y=PC2)) + geom_point(aes(color=donor_id), size=8) +
xlab(paste0("PC1:", round(pcs_naive$pVar[1]*100), "% variance")) +
ylab(paste0("PC2: ", round(pcs_naive$pVar[2]*100), "% variance")) +
coord_fixed() + theme_bw() +theme(plot.title=element_text(size=20, hjust=0.5), axis.title=element_text(size=14), panel.grid = element_blank(), axis.text=element_text(size=12),legend.text=element_text(size=12), legend.title=element_text(size=12), legend.key.size = unit(1.5,"lines"))

去掉个体间变异性：

assay(proteomics_data_naive) <- removeBatchEffect(assay(proteomics_data_naive),
batch = factor(as.vector(colData(proteomics_data_naive)$donor_id))
)
pcs_naive <- getPCs(proteomics_data_naive)
ggplot(data=pcs_naive$pcs, aes(x=PC1, y=PC2, color=cytokine_condition)) +
geom_point(size = 8) + geom_label_repel(aes(label=cytokine_condition, color=cytokine_condition)) +
xlab(paste0("PC1: ", round(pcs_naive$pVar[1]*100), "% variance")) +
ylab(paste0("PC2: ", round(pcs_naive$pVar[2]*100), "% variance")) +
scale_colour_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
coord_fixed() + theme_bw() +theme(panel.grid = element_blank(), legend.position = "none")

删除由PCA标识的异常样本:

proteomics_data_naive <- proteomics_data_naive[, colnames(proteomics_data_naive) != "D257_CD4_naive_Th1"]
pcs_naive <- getPCs(proteomics_data_naive)
ggplot(data=pcs_naive$pcs, aes(x=PC1, y=PC2, color=cytokine_condition)) +
geom_point(size = 8) + geom_label_repel(aes(label=cytokine_condition, color=cytokine_condition)) +
xlab(paste0("PC1: ", round(pcs_naive$pVar[1]*100), "% variance")) +
ylab(paste0("PC2: ", round(pcs_naive$pVar[2]*100), "% variance")) +
scale_colour_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
coord_fixed() + theme_bw() +theme(panel.grid = element_blank(), legend.position = "none")

原图

Memory T cells

again。。。

Performing PCA on 5 days-stimulated memory T cells only.
```{r compute_pca_naive, message=FALSE, warning=FALSE}pcs_memory <- getPCs(proteomics_data_memory)
ggplot(data=pcs_memory$pcs, aes(x=PC1, y=PC2)) + geom_point(aes(color=donor_id), size=8) +
xlab(paste0("PC1:", round(pcs_memory$pVar[1]*100), "% variance")) +
ylab(paste0("PC2: ", round(pcs_memory$pVar[2]*100), "% variance")) +
coord_fixed() + theme_bw() +theme(plot.title=element_text(size=20, hjust=0.5), axis.title=element_text(size=14), panel.grid = element_blank(), axis.text=element_text(size=12),legend.text=element_text(size=12), legend.title=element_text(size=12), legend.key.size = unit(1.5,"lines"))

Regressing out inter-individual variability

assay(proteomics_data_memory) <- removeBatchEffect(assay(proteomics_data_memory),
batch = factor(as.vector(colData(proteomics_data_memory)$donor_id)))

再次计算PCs

pcs_memory <- getPCs(proteomics_data_memory)
ggplot(data=pcs_memory$pcs, aes(x=PC1, y=PC2, color=cytokine_condition)) +
geom_point(size = 8) + geom_label_repel(aes(label=cytokine_condition, color=cytokine_condition)) +
xlab(paste0("PC1: ", round(pcs_memory$pVar[1]*100), "% variance")) +
ylab(paste0("PC2: ", round(pcs_memory$pVar[2]*100), "% variance")) +
scale_colour_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
coord_fixed() + theme_bw() +theme(panel.grid = element_blank(), legend.position = "none")

原图

基本分布还是差不多的，，，，

快去试一试呀！

你可能还想看

  • 蛋白质组学研究概述

  • PCA主成分分析实战和可视化 附R代码和测试数据

  • 用了这么多年的PCA可视化竟然是错的！！！

  • 什么？你做的差异基因方法不合适？

  • NBT：单细胞转录组新降维可视化方法PHATE


    • ©著作权归作者所有,转载或内容合作请联系作者
    • 人面猴
      序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
      沈念sama阅读 204,590评论 6赞 478
    • 死咒
      序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
      沈念sama阅读 86,808评论 2赞 381
    • 救了他两次的神仙让他今天三更去死
      文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
      开封第一讲书人阅读 151,151评论 0赞 337
    • 道士缉凶录：失踪的卖姜人
      文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
      开封第一讲书人阅读 54,779评论 1赞 277
    • 港岛之恋（遗憾婚礼）
      正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
      茶点故事阅读 63,773评论 5赞 367
    • 恶毒庶女顶嫁案：这布局不是一般人想出来的
      文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
      开封第一讲书人阅读 48,656评论 1赞 281
    • 城市分裂传说
      那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
      沈念sama阅读 38,022评论 3赞 398
    • 双鸳鸯连环套：你想象不到人心有多黑
      文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
      开封第一讲书人阅读 36,678评论 0赞 258
    • 万荣杀人案实录
      序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
      沈念sama阅读 41,038评论 1赞 299
    • 护林员之死
      正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
      茶点故事阅读 35,659评论 2赞 321
    • 白月光启示录
      正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
      茶点故事阅读 37,756评论 1赞 330
    • 活死人
      序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
      沈念sama阅读 33,411评论 4赞 321
    • 日本核电站爆炸内幕
      正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
      茶点故事阅读 39,005评论 3赞 307
    • 男人毒药：我在死后第九天来索命
      文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
      开封第一讲书人阅读 29,973评论 0赞 19
    • 一桩弑父案，背后竟有这般阴谋
      文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
      开封第一讲书人阅读 31,203评论 1赞 260
    • 情欲美人皮
      我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
      沈念sama阅读 45,053评论 2赞 350
    • 代替公主和亲
      正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
      茶点故事阅读 42,495评论 2赞 343