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，每行代表一個蛋白
創建帶有蛋白質注釋的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)
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"
)
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
定義函數：
提取蛋白質表達值;
進行主成分分析;
返回一個矩陣，其中包含每個樣品和樣品注釋的PC坐標;
返回每個主要成分解釋的方差百分比。
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)
}
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")]
對 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")
原圖
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")
原圖
基本分布還是差不多的，，，，
快去試一試呀！
你可能還想看往期精品(點擊圖片直達文字對應教程)後臺回復「生信寶典福利第一波」或點擊閱讀原文獲取教程合集