蛋白组学分析包——DEqMS学习

url <- "https://ftp.ebi.ac.uk/pride-archive/2016/06/PXD004163/Yan_miR_Protein_table.flatprottable.txt"

download.file(url, destfile = "./miR_Proteintable.txt",method = "auto")

df.prot = read.table("miR_Proteintable.txt",stringsAsFactors = FALSE, header = TRUE, quote = "", comment.char = "",sep = "\t") 

# filter at 1% protein FDR and extract TMT quantifications
TMT_columns = seq(15,33,2)      #选择15-33列中、间隔1列的数据，也就是所有的COUNT数据
dat = df.prot[df.prot$miR.FASP_q.value<0.01,TMT_columns]    #筛选q值小于0.01的数据
rownames(dat) = df.prot[df.prot$miR.FASP_q.value<0.01,]$Protein.accession

# The protein dataframe is a typical protein expression matrix structure
# Samples are in columns, proteins are in rows
# use unique protein IDs for rownames
# to view the whole data frame, use the command View(dat)

View(dat)

dat.log = log2(dat)    #数据取对数
dat.log = na.omit(dat.log)    #remove rows with NAs
View(dat.log)
boxplot(dat.log,las=2,main="TMT10plex data PXD004163")    #Use boxplot to check if the samples have medians centered. if not, do median centering.

# Here the data is already median centered, we skip the following step. 
#由于数据表达量足够一致，所以可以不进行标准化
# dat.log = equalMedianNormalization(dat.log)

# if there is only one factor, such as treatment. You can define a vector with
# the treatment group in the same order as samples in the protein table.

cond = as.factor(c("ctrl","miR191","miR372","miR519","ctrl",
                   "miR372","miR519","ctrl","miR191","miR372"))

# The function model.matrix is used to generate the design matrix
design = model.matrix(~0+cond) # 0 means no intercept for the linear model
#这里是设计一个对比矩阵，其做法原理可以详见：https://treeh.cn/?id=21

colnames(design) = gsub("cond","",colnames(design))    #去除表格中的cond字符
View(design)

# you can define one or multiple contrasts here
#设计哪些组间需要对比
x <- c("miR372-ctrl","miR519-ctrl","miR191-ctrl", "miR372-miR519","miR372-miR191","miR519-miR191")
contrast =  makeContrasts(contrasts=x,levels=design)    #按照x的对比方式，对design的样本分组对样本进行比较
View(contrast)

fit1 <- lmFit(dat.log, design)    #线性拟合模型构建
fit2 <- contrasts.fit(fit1,contrasts = contrast)    #Compute Contrasts from Linear Model Fit

#Given a linear model fit to microarray data, compute estimated coefficients and standard errors for a given set of contrasts.
fit3 <- eBayes(fit2)#Empirical Bayes Statistics for Differential Expression
#Given a linear model fit from lmFit, compute moderated t-statistics, moderated F-statistic, and log-odds of differential expression by empirical Bayes moderation of the standard errors towards a global value.
#具体介绍参考自help(contrasts.fit)、help(eBayes)

# assign a extra variable `count` to fit3 object, telling how many PSMs are quantifed for each protein
library(matrixStats)
count_columns = seq(16,34,2)     #选择16-34列中、间隔1列的数据，也就是所有的PSMs数据
psm.count.table = data.frame(count = rowMins(as.matrix(df.prot[,count_columns])), row.names =  df.prot$Protein.accession)

#rowMins: Calculates the minimum for each row (column) of a matrix-like object

fit3$count = psm.count.table[rownames(fit3$coefficients),"count"]    #数据导入fit3中
fit4 = spectraCounteBayes(fit3)
#Peptide/Spectra Count Based Empirical Bayes Statistics for Differential Expression. This function is to calculate peptide/PSM count adjusted t-statistics, p-values.
View(psm.count.table)

# n=30 limits the boxplot to show only proteins quantified by <= 30 PSMs.

VarianceBoxplot(fit4,n=30,main="TMT10plex dataset PXD004163",xlab="PSM count")

#横坐标为不同PSM count值，纵坐标为log值，但是是什么数据的Log值并不清楚
#查看VarianceBoxplot源码，x轴y轴输入如下
#x <- fit$count
#y <- fit$sigma^2，这个sigma值不知道是什么，应该是涉及上机检测的数据？在公司给的数据里面好像有见到类似值？等待一个解答***

VarianceScatterplot(fit4,main="TMT10plex dataset PXD004163")

#if you are not sure which coef_col refers to the specific contrast,type
#查看我们设定的对照组
head(fit4$coefficients)

#提取第一个对比组
DEqMS.results = outputResult(fit4,coef_col = 1)#miR372-ctrl
#提取第二个对比组
#DEqMS.results = outputResult(fit4,coef_col = 2)#miR519-ctrl
# a quick look on the DEqMS results table
head(DEqMS.results)
# Save it into a tabular text file
write.table(DEqMS.results,"DEqMS.results.miR372-ctrl.txt",sep = "\t", row.names = F,quote=F)