## ----Loadpackage----------------------------------------------------------- library(DEqMS) ## ----DownloadProteinTable-------------------------------------------------- url <- "ftp://ftp.pride.ebi.ac.uk/pride/data/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") ## ----Extractcolumn--------------------------------------------------------- # filter at 1% protein FDR and extract TMT quantifications TMT_columns = seq(15,33,2) dat = df.prot[df.prot$miR.FASP_q.value<0.01,TMT_columns] 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) ## ----log2transform1-------------------------------------------------------- dat.log = log2(dat) #remove rows with NAs dat.log = na.omit(dat.log) ## ----boxplot1-------------------------------------------------------------- boxplot(dat.log,las=2,main="TMT10plex data PXD004163") # Here the data is already median centered, we skip the following step. # dat.log = equalMedianNormalization(dat.log) ## ----design---------------------------------------------------------------- # 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 colnames(design) = gsub("cond","",colnames(design)) ## ----limma----------------------------------------------------------------- # 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) fit1 <- lmFit(dat.log, design) fit2 <- contrasts.fit(fit1,contrasts = contrast) fit3 <- eBayes(fit2) ## ---- DEqMS1--------------------------------------------------------------- # 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) psm.count.table = data.frame(count = rowMins( as.matrix(df.prot[,count_columns])), row.names = df.prot$Protein.accession) fit3$count = psm.count.table[rownames(fit3$coefficients),"count"] fit4 = spectraCounteBayes(fit3) ## ---- plot----------------------------------------------------------------- # n=30 limits the boxplot to show only proteins quantified by <= 30 PSMs. VarianceBoxplot(fit4,n=30,main="TMT10plex dataset PXD004163",xlab="PSM count") VarianceScatterplot(fit4,main="TMT10plex dataset PXD004163") ## ---- result--------------------------------------------------------------- DEqMS.results = outputResult(fit4,coef_col = 1) #if you are not sure which coef_col refers to the specific contrast,type head(fit4$coefficients) # 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) ## ----volcanoplot1---------------------------------------------------------- library(ggrepel) # Use ggplot2 allows more flexibility in plotting DEqMS.results$log.sca.pval = -log10(DEqMS.results$sca.P.Value) ggplot(DEqMS.results, aes(x = logFC, y =log.sca.pval )) + geom_point(size=0.5 )+ theme_bw(base_size = 16) + # change theme xlab(expression("log2(miR372/ctrl)")) + # x-axis label ylab(expression(" -log10(P-value)")) + # y-axis label geom_vline(xintercept = c(-1,1), colour = "red") + # Add fold change cutoffs geom_hline(yintercept = 3, colour = "red") + # Add significance cutoffs geom_vline(xintercept = 0, colour = "black") + # Add 0 lines scale_colour_gradient(low = "black", high = "black", guide = FALSE)+ geom_text_repel(data=subset(DEqMS.results, abs(logFC)>1&log.sca.pval > 3), aes( logFC, log.sca.pval ,label=gene)) # add gene label ## ----volcanoplot2---------------------------------------------------------- fit4$p.value = fit4$sca.p # volcanoplot highlight top 20 proteins ranked by p-value here volcanoplot(fit4,coef=1, style = "p-value", highlight = 20, names=rownames(fit4$coefficients)) ## ----DownloadLabelfreeData------------------------------------------------- url2 <- "ftp://ftp.pride.ebi.ac.uk/pride/data/archive/2014/09/PXD000279/proteomebenchmark.zip" download.file(url2, destfile = "./PXD000279.zip",method = "auto") unzip("PXD000279.zip") ## ----LFQprotein------------------------------------------------------------ df.prot = read.table("proteinGroups.txt",header=T,sep="\t",stringsAsFactors = F, comment.char = "",quote ="") # remove decoy matches and matches to contaminant df.prot = df.prot[!df.prot$Reverse=="+",] df.prot = df.prot[!df.prot$Contaminant=="+",] # Extract columns of LFQ intensites df.LFQ = df.prot[,89:94] df.LFQ[df.LFQ==0] <- NA rownames(df.LFQ) = df.prot$Majority.protein.IDs df.LFQ$na_count_H = apply(df.LFQ,1,function(x) sum(is.na(x[1:3]))) df.LFQ$na_count_L = apply(df.LFQ,1,function(x) sum(is.na(x[4:6]))) # Filter protein table. DEqMS require minimum two values for each group. df.LFQ.filter = df.LFQ[df.LFQ$na_count_H<2 & df.LFQ$na_count_L<2,1:6] ## ----pepCountTable--------------------------------------------------------- library(matrixStats) # we use minimum peptide count among six samples # count unique+razor peptides used for quantification pep.count.table = data.frame(count = rowMins(as.matrix(df.prot[,19:24])), row.names = df.prot$Majority.protein.IDs) # Minimum peptide count of some proteins can be 0 # add pseudocount 1 to all proteins pep.count.table$count = pep.count.table$count+1 ## ----labelfreeDEqMS-------------------------------------------------------- protein.matrix = log2(as.matrix(df.LFQ.filter)) class = as.factor(c("H","H","H","L","L","L")) design = model.matrix(~0+class) # fitting without intercept fit1 = lmFit(protein.matrix,design = design) cont <- makeContrasts(classH-classL, levels = design) fit2 = contrasts.fit(fit1,contrasts = cont) fit3 <- eBayes(fit2) fit3$count = pep.count.table[rownames(fit3$coefficients),"count"] #check the values in the vector fit3$count #if min(fit3$count) return NA or 0, you should troubleshoot the error first min(fit3$count) fit4 = spectraCounteBayes(fit3) ## ----LFQboxplot------------------------------------------------------------ VarianceBoxplot(fit4, n=20, main = "Label-free dataset PXD000279", xlab="peptide count + 1") ## ----LFQresult------------------------------------------------------------- DEqMS.results = outputResult(fit4,coef_col = 1) # Add Gene names to the data frame rownames(df.prot) = df.prot$Majority.protein.IDs DEqMS.results$Gene.name = df.prot[DEqMS.results$gene,]$Gene.names head(DEqMS.results) write.table(DEqMS.results,"H-L.DEqMS.results.txt",sep = "\t", row.names = F,quote=F) ## ----retrieveExampleData, message=FALSE------------------------------------ ### retrieve example PSM dataset from ExperimentHub library(ExperimentHub) eh = ExperimentHub() query(eh, "DEqMS") dat.psm = eh[["EH1663"]] ## ----log2transform2-------------------------------------------------------- dat.psm.log = dat.psm dat.psm.log[,3:12] = log2(dat.psm[,3:12]) head(dat.psm.log) ## ----boxplot2-------------------------------------------------------------- dat.gene.nm = medianSweeping(dat.psm.log,group_col = 2) boxplot(dat.gene.nm,las=2,ylab="log2 ratio",main="TMT10plex dataset PXD004163") ## ----DEqMS2---------------------------------------------------------------- gene.matrix = as.matrix(dat.gene.nm) # make design table cond = as.factor(c("ctrl","miR191","miR372","miR519","ctrl", "miR372","miR519","ctrl","miR191","miR372")) design = model.matrix(~0+cond) colnames(design) = gsub("cond","",colnames(design)) #limma part analysis fit1 <- lmFit(gene.matrix,design) x <- c("miR372-ctrl","miR519-ctrl","miR191-ctrl") contrast = makeContrasts(contrasts=x,levels=design) fit2 <- eBayes(contrasts.fit(fit1,contrasts = contrast)) #DEqMS part analysis psm.count.table = as.data.frame(table(dat.psm$gene)) rownames(psm.count.table) = psm.count.table$Var1 fit2$count = psm.count.table[rownames(fit2$coefficients),2] fit3 = spectraCounteBayes(fit2) # extract DEqMS results DEqMS.results = outputResult(fit3,coef_col = 1) head(DEqMS.results) write.table(DEqMS.results,"DEqMS.results.miR372-ctrl.fromPSMtable.txt", sep = "\t",row.names = F,quote=F) ## ----PriorVarianceTrend---------------------------------------------------- VarianceBoxplot(fit3,n=20, xlab="PSM count",main="TMT10plex dataset PXD004163") ## ----PSMintensity---------------------------------------------------------- peptideProfilePlot(dat=dat.psm.log,col=2,gene="TGFBR2") # col=2 is tell in which column of dat.psm.log to look for the gene ## ----PriorVar-------------------------------------------------------------- VarianceScatterplot(fit3, xlab="log2(PSM count)") limma.prior = fit3$s2.prior abline(h = log(limma.prior),col="green",lwd=3 ) legend("topright",legend=c("DEqMS prior variance","Limma prior variance"), col=c("red","green"),lwd=3) ## ----Residualplot---------------------------------------------------------- op <- par(mfrow=c(1,2), mar=c(4,4,4,1), oma=c(0.5,0.5,0.5,0)) Residualplot(fit3, xlab="log2(PSM count)",main="DEqMS") x = fit3$count y = log(limma.prior) - log(fit3$sigma^2) plot(log2(x),y,ylim=c(-6,2),ylab="Variance(estimated-observed)", pch=20, cex=0.5, xlab = "log2(PSMcount)",main="Limma") ## ----PostVar, echo=TRUE, fig.height=5, fig.width=10------------------------ library(LSD) op <- par(mfrow=c(1,2), mar=c(4,4,4,1), oma=c(0.5,0.5,0.5,0)) x = fit3$count y = fit3$s2.post heatscatter(log2(x),log(y),pch=20, xlab = "log2(PSMcount)", ylab="log(Variance)", main="Posterior Variance in Limma") y = fit3$sca.postvar heatscatter(log2(x),log(y),pch=20, xlab = "log2(PSMcount)", ylab="log(Variance)", main="Posterior Variance in DEqMS") ## ----t-test---------------------------------------------------------------- pval.372 = apply(dat.gene.nm, 1, function(x) t.test(as.numeric(x[c(1,5,8)]), as.numeric(x[c(3,6,10)]))$p.value) logFC.372 = rowMeans(dat.gene.nm[,c(3,6,10)])-rowMeans(dat.gene.nm[,c(1,5,8)]) ## ----echo=TRUE------------------------------------------------------------- ttest.results = data.frame(gene=rownames(dat.gene.nm), logFC=logFC.372,P.Value = pval.372, adj.pval = p.adjust(pval.372,method = "BH")) ttest.results$PSMcount = psm.count.table[ttest.results$gene,"count"] ttest.results = ttest.results[with(ttest.results, order(P.Value)), ] head(ttest.results) ## ----Anova----------------------------------------------------------------- ord.t = fit1$coefficients[, 1]/fit1$sigma/fit1$stdev.unscaled[, 1] ord.p = 2*pt(abs(ord.t), fit1$df.residual, lower.tail = FALSE) ord.q = p.adjust(ord.p,method = "BH") anova.results = data.frame(gene=names(fit1$sigma), logFC=fit1$coefficients[,1], t=ord.t, P.Value=ord.p, adj.P.Val = ord.q) anova.results$PSMcount = psm.count.table[anova.results$gene,"count"] anova.results = anova.results[with(anova.results,order(P.Value)),] head(anova.results) ## ----echo=TRUE------------------------------------------------------------- limma.results = topTable(fit2,coef = 1,n= Inf) limma.results$gene = rownames(limma.results) #Add PSM count values in the data frame limma.results$PSMcount = psm.count.table[limma.results$gene,"count"] head(limma.results) ## ----pvalueall------------------------------------------------------------- plot(sort(-log10(limma.results$P.Value),decreasing = TRUE), type="l",lty=2,lwd=2, ylab="-log10(p-value)",ylim = c(0,10), xlab="Proteins ranked by p-values", col="purple") lines(sort(-log10(DEqMS.results$sca.P.Value),decreasing = TRUE), lty=1,lwd=2,col="red") lines(sort(-log10(anova.results$P.Value),decreasing = TRUE), lty=2,lwd=2,col="blue") lines(sort(-log10(ttest.results$P.Value),decreasing = TRUE), lty=2,lwd=2,col="orange") legend("topright",legend = c("Limma","DEqMS","Anova","t.test"), col = c("purple","red","blue","orange"),lty=c(2,1,2,2),lwd=2) ## ----pvalue500------------------------------------------------------------- plot(sort(-log10(limma.results$P.Value),decreasing = TRUE)[1:500], type="l",lty=2,lwd=2, ylab="-log10(p-value)", ylim = c(2,10), xlab="Proteins ranked by p-values", col="purple") lines(sort(-log10(DEqMS.results$sca.P.Value),decreasing = TRUE)[1:500], lty=1,lwd=2,col="red") lines(sort(-log10(anova.results$P.Value),decreasing = TRUE)[1:500], lty=2,lwd=2,col="blue") lines(sort(-log10(ttest.results$P.Value),decreasing = TRUE)[1:500], lty=2,lwd=2,col="orange") legend("topright",legend = c("Limma","DEqMS","Anova","t.test"), col = c("purple","red","blue","orange"),lty=c(2,1,2,2),lwd=2)