復(fù)現(xiàn)PCA原圖之蛋白組學(xué)數(shù)據(jù)

NGS系列文章包括NGS基礎(chǔ)、轉(zhuǎn)錄組分析?（Nature重磅綜述|關(guān)于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細(xì)胞測(cè)序分析?(重磅綜述：三萬(wàn)字長(zhǎng)文讀懂單細(xì)胞RNA測(cè)序分析的最佳實(shí)踐教程 （原理、代碼和評(píng)述）)、DNA甲基化分析、重測(cè)序分析、GEO數(shù)據(jù)挖掘（典型醫(yī)學(xué)設(shè)計(jì)實(shí)驗(yàn)GEO數(shù)據(jù)分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內(nèi)容。

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

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

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

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

返回每個(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ò)的！！！

• 什么？你做的差異基因方法不合適？

• NBT：單細(xì)胞轉(zhuǎn)錄組新降維可視化方法PHATE

往期精品(點(diǎn)擊圖片直達(dá)文字對(duì)應(yīng)教程)

后臺(tái)回復(fù)“生信寶典福利第一波”或點(diǎn)擊閱讀原文獲取教程合集

總結(jié)

以上是生活随笔為你收集整理的这篇Nature子刊文章的蛋白组学数据PCA分析竟花费了我两天时间来重现|附全过程代码...的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。