這篇Nature子刊文章的蛋白組學(xué)數(shù)據(jù)PCA分析竟花費(fèi)了我兩天時(shí)間來(lái)重現(xiàn)|附全過(guò)程代碼

2020年4月14日,Sanger研究團(tuán)隊(duì)于nature communication在線發(fā)表了題為Single-cell transcriptomics identifies an effectorness gradient shaping the response of CD4+ T cells to cytokines的研究?jī)?nèi)容维哈,作者使用蛋白質(zhì)組學(xué)、bulk RNA-seq和單細(xì)胞轉(zhuǎn)錄組測(cè)序?qū)θ梭w40,000個(gè)以上的na?ve and memory CD4+ T cells進(jìn)行分析放仗,發(fā)現(xiàn)細(xì)胞類型之間的細(xì)胞因子反應(yīng)差異很大秉撇。memory T細(xì)胞不能分化為Th2表型灵迫,但可以響應(yīng)iTreg極化獲得類似Th17的表型照雁。單細(xì)胞分析表明蚕愤,T細(xì)胞構(gòu)成了一個(gè)轉(zhuǎn)錄連續(xù)體(transcriptional continuum)答恶,從幼稚到中樞和效應(yīng)記憶T細(xì)胞饺蚊,形成了一種效應(yīng)梯度,并伴隨著趨化因子和細(xì)胞因子表達(dá)的增加悬嗓。最后污呼,作者表明,T細(xì)胞活化和細(xì)胞因子反應(yīng)受效應(yīng)梯度的影響包竹。

該文獻(xiàn)通過(guò)蛋白質(zhì)組學(xué)((液相色譜-串聯(lián)質(zhì)譜法燕酷,LC-MS/MS)進(jìn)行了探索性分析,樣品對(duì)應(yīng)于從健康個(gè)體的外周血中分離的幼稚和記憶T細(xì)胞周瞎,并用多種細(xì)胞因子刺激5天苗缩,每個(gè)條件平均3個(gè)生物學(xué)重復(fù)。

這次復(fù)現(xiàn)Fig1cPCA圖和Fig2aPCA圖的另一部分声诸,這次作者是通過(guò)蛋白組學(xué)數(shù)據(jù)進(jìn)行PCA的展現(xiàn):

以上是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 ”

我們需要復(fù)現(xiàn)該圖之前彼乌,先需要下載數(shù)據(jù)泻肯,可以點(diǎn)擊https://www.opentargets.org/projects/effectorness對(duì)proteomics的abundances數(shù)據(jù)和metadata數(shù)據(jù)進(jìn)行下載渊迁,然后進(jìn)行以下步驟:

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

加載數(shù)據(jù)

加載標(biāo)準(zhǔn)化后的豐度:

MassSpec_data <- read.table("NCOMMS-19-7936188_MassSpec_scaled_abundances.txt", header = T, stringsAsFactors = F)
View(MassSpec_data)#從以上可以看出,每列除了代表每個(gè)樣本外灶挟,前三列分別為Protein_id琉朽,Gene_id和Gene_name,每行代表一個(gè)蛋白

建立SummarizedExperiment object

創(chuàng)建帶有蛋白質(zhì)注釋的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#構(gòu)成一個(gè)由"Protein_id","Gene_id","Gene_name"的數(shù)據(jù)框MassSpec_data <- MassSpec_data[,-c(1:3)]

創(chuàng)建帶有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)

創(chuàng)建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")

數(shù)據(jù)可視化

將NA值設(shè)置為零
注意:此操作僅出于可視化目的稚铣。執(zhí)行統(tǒng)計(jì)測(cè)試時(shí)箱叁,NA不會(huì)設(shè)置為零。

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

定義函數(shù):

\bullet 提取蛋白質(zhì)表達(dá)值;

\bullet 進(jìn)行主成分分析;

\bullet 返回一個(gè)矩陣惕医,其中包含每個(gè)樣品和樣品注釋的PC坐標(biāo);

\bullet 返回每個(gè)主要成分解釋的方差百分比蝌蹂。

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)}

對(duì)所有樣本執(zhí)行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())

去掉個(gè)體間變異性:

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

重新計(jì)算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())

原圖

細(xì)胞類型特異性分析

將naive和memory T細(xì)胞樣本分為僅包含受刺激細(xì)胞的兩個(gè)不同數(shù)據(jù)集。

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

對(duì) 5 days-stimulated naive T cells進(jìn)行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"))

去掉個(gè)體間變異性:

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標(biāo)識(shí)的異常樣本:

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)))

再次計(jì)算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")

原圖

基本分布還是差不多的沛简,齐鲤,,椒楣,

快去試一試呀给郊!

你可能還想看

  • 蛋白質(zhì)組學(xué)研究概述

  • PCA主成分分析實(shí)戰(zhàn)和可視化 附R代碼和測(cè)試數(shù)據(jù)

  • 用了這么多年的PCA可視化竟然是錯(cuò)的!E趸摇淆九!

  • 什么?你做的差異基因方法不合適毛俏?

  • NBT:?jiǎn)渭?xì)胞轉(zhuǎn)錄組新降維可視化方法PHATE


  • ?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
    • 序言:七十年代末炭庙,一起剝皮案震驚了整個(gè)濱河市,隨后出現(xiàn)的幾起案子煌寇,更是在濱河造成了極大的恐慌焕蹄,老刑警劉巖,帶你破解...
      沈念sama閱讀 206,013評(píng)論 6 481
    • 序言:濱河連續(xù)發(fā)生了三起死亡事件阀溶,死亡現(xiàn)場(chǎng)離奇詭異腻脏,居然都是意外死亡,警方通過(guò)查閱死者的電腦和手機(jī)银锻,發(fā)現(xiàn)死者居然都...
      沈念sama閱讀 88,205評(píng)論 2 382
    • 文/潘曉璐 我一進(jìn)店門永品,熙熙樓的掌柜王于貴愁眉苦臉地迎上來(lái),“玉大人击纬,你說(shuō)我怎么就攤上這事鼎姐。” “怎么了?”我有些...
      開(kāi)封第一講書人閱讀 152,370評(píng)論 0 342
    • 文/不壞的土叔 我叫張陵症见,是天一觀的道長(zhǎng)喂走。 經(jīng)常有香客問(wèn)我,道長(zhǎng)谋作,這世上最難降的妖魔是什么芋肠? 我笑而不...
      開(kāi)封第一講書人閱讀 55,168評(píng)論 1 278
    • 正文 為了忘掉前任,我火速辦了婚禮遵蚜,結(jié)果婚禮上帖池,老公的妹妹穿的比我還像新娘。我一直安慰自己吭净,他們只是感情好睡汹,可當(dāng)我...
      茶點(diǎn)故事閱讀 64,153評(píng)論 5 371
    • 文/花漫 我一把揭開(kāi)白布。 她就那樣靜靜地躺著寂殉,像睡著了一般囚巴。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上友扰,一...
      開(kāi)封第一講書人閱讀 48,954評(píng)論 1 283
    • 那天彤叉,我揣著相機(jī)與錄音,去河邊找鬼村怪。 笑死秽浇,一個(gè)胖子當(dāng)著我的面吹牛,可吹牛的內(nèi)容都是我干的甚负。 我是一名探鬼主播柬焕,決...
      沈念sama閱讀 38,271評(píng)論 3 399
    • 文/蒼蘭香墨 我猛地睜開(kāi)眼,長(zhǎng)吁一口氣:“原來(lái)是場(chǎng)噩夢(mèng)啊……” “哼梭域!你這毒婦竟也來(lái)了斑举?” 一聲冷哼從身側(cè)響起,我...
      開(kāi)封第一講書人閱讀 36,916評(píng)論 0 259
    • 序言:老撾萬(wàn)榮一對(duì)情侶失蹤碰辅,失蹤者是張志新(化名)和其女友劉穎懂昂,沒(méi)想到半個(gè)月后介时,有當(dāng)?shù)厝嗽跇?shù)林里發(fā)現(xiàn)了一具尸體没宾,經(jīng)...
      沈念sama閱讀 43,382評(píng)論 1 300
    • 正文 獨(dú)居荒郊野嶺守林人離奇死亡,尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
      茶點(diǎn)故事閱讀 35,877評(píng)論 2 323
    • 正文 我和宋清朗相戀三年驰坊,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了践惑。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片蟋座。...
      茶點(diǎn)故事閱讀 37,989評(píng)論 1 333
    • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡,死狀恐怖会钝,靈堂內(nèi)的尸體忽然破棺而出,到底是詐尸還是另有隱情,我是刑警寧澤迁酸,帶...
      沈念sama閱讀 33,624評(píng)論 4 322
    • 正文 年R本政府宣布先鱼,位于F島的核電站,受9級(jí)特大地震影響奸鬓,放射性物質(zhì)發(fā)生泄漏焙畔。R本人自食惡果不足惜,卻給世界環(huán)境...
      茶點(diǎn)故事閱讀 39,209評(píng)論 3 307
    • 文/蒙蒙 一串远、第九天 我趴在偏房一處隱蔽的房頂上張望宏多。 院中可真熱鬧,春花似錦澡罚、人聲如沸伸但。這莊子的主人今日做“春日...
      開(kāi)封第一講書人閱讀 30,199評(píng)論 0 19
    • 文/蒼蘭香墨 我抬頭看了看天上的太陽(yáng)更胖。三九已至,卻和暖如春隔显,著一層夾襖步出監(jiān)牢的瞬間函喉,已是汗流浹背。 一陣腳步聲響...
      開(kāi)封第一講書人閱讀 31,418評(píng)論 1 260
    • 我被黑心中介騙來(lái)泰國(guó)打工荣月, 沒(méi)想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留管呵,地道東北人。 一個(gè)月前我還...
      沈念sama閱讀 45,401評(píng)論 2 352
    • 正文 我出身青樓哺窄,卻偏偏與公主長(zhǎng)得像捐下,于是被迫代替她去往敵國(guó)和親。 傳聞我的和親對(duì)象是個(gè)殘疾皇子萌业,可洞房花燭夜當(dāng)晚...
      茶點(diǎn)故事閱讀 42,700評(píng)論 2 345