Contents

1 Background

Duguet et al. 2017 compared the proteomes of mouse regulatory T cells (Treg) and conventional T cells (Tconv) in order to discover differentially regulated proteins between these two cell populations. For each biological repeat the proteomes were extracted for both Treg and Tconv cell pools, which were purified by flow cytometry. The data in data/quantification/mouseTcell on the pdaData repository are a subset of the data PXD004436 on PRIDE.

We will use a subset of the data with a randomized complete block (RCB) design, i.e. the dataset consists of four mice for which the proteome of both conventional and regulatory T cells are assessed.

We will analyse the data in this file as a completely randomized design, to illustrate the impact of a wrong analyssi that is not accounting for the blocking factor for mouse.

2 Data

We first import the peptides.txt file. This is the file that contains your peptide-level intensities. For a MaxQuant search [6], this peptides.txt file can be found by default in the “path_to_raw_files/combined/txt/” folder from the MaxQuant output, with “path_to_raw_files” the folder where raw files were saved. In this tutorial, we will use a MaxQuant peptides file from MaxQuant that can be found on the pdaData repository. We will use the MSnbase package to import the data.

We generate the object peptideFile with the path to the peptides.txt file. In this file we will fetch the data from the github repository linked to the pda course: https://statomics.github.io/pda. You can also replace the peptideFile with a string that points to the path of a file on your local hard drive.

With the grepEcols function we find the columns that are containing the expression data of the peptides in the peptides.txt file.

library(MSqRob)
library(MSnbase)
library(tidyverse)
library(limma)
peptidesFile <- "https://raw.githubusercontent.com/statOmics/pda/data/quantification/mouseTcell/peptidesRCB.txt"
proteinGroupsFile <-"https://raw.githubusercontent.com/statOmics/pda/data/quantification/mouseTcell/proteinGroups.txt"
ecols<-grepEcols(peptidesFile, "Intensity ", split = "\t")

Next we import the peptide intensities

pepData<-readMSnSet2(peptidesFile,ecol=ecols,fnames="Sequence",sep="\t")
pepData
## MSnSet (storageMode: lockedEnvironment)
## assayData: 55814 features, 8 samples 
##   element names: exprs 
## protocolData: none
## phenoData: none
## featureData
##   featureNames: AAAAAAAAAAGAAGGR AAAAAAAAAAGDSDSWDADTFSMEDPVRK ...
##     YYYDKNIIHK (55814 total)
##   fvarLabels: Sequence N.term.cleavage.window ... MS.MS.Count (74
##     total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:  
## - - - Processing information - - -
##  MSnbase version: 2.10.1

The pepData object is an MSnSet, a container for the data, features information and experimental annotation. They can be accessed using the accessor functions ‘exprs’ (matrix of intensities, features in rows, samples in columns), ‘fData’ (properties for each feature, peptide or protein, on the rows) and ‘pData’ (properties for the samples on the columns).

We will make use from data wrangling functionalities from the tidyverse package. The %>% operator allows us to pipe the output of one function to the next function.

head(exprs(pepData))
##                               Intensity.Tconv.M12_2 Intensity.Tconv.M12_3
## AAAAAAAAAAGAAGGR                          111470000             101350000
## AAAAAAAAAAGDSDSWDADTFSMEDPVRK              13303000              10160000
## AAAAAAAGGAAIAVSTGIETATIQK                         0                     0
## AAAAAAGAASGIPGPVAQGIK                             0               7197400
## AAAAAAGPEMVR                              104440000             101950000
## AAAAAAPGGGGGEPR                            11910000              12683000
##                               Intensity.Tconv.M5_inj1
## AAAAAAAAAAGAAGGR                               965070
## AAAAAAAAAAGDSDSWDADTFSMEDPVRK                  235500
## AAAAAAAGGAAIAVSTGIETATIQK                       48225
## AAAAAAGAASGIPGPVAQGIK                          266620
## AAAAAAGPEMVR                                   947110
## AAAAAAPGGGGGEPR                                190550
##                               Intensity.Tconv.M6_inj1 Intensity.Treg.M12_2
## AAAAAAAAAAGAAGGR                                    0             80606000
## AAAAAAAAAAGDSDSWDADTFSMEDPVRK                 6744500                    0
## AAAAAAAGGAAIAVSTGIETATIQK                      314120                    0
## AAAAAAGAASGIPGPVAQGIK                         2904800             12602000
## AAAAAAGPEMVR                                 12747000             44735000
## AAAAAAPGGGGGEPR                               2166100             12555000
##                               Intensity.Treg.M12_3 Intensity.Treg.M5_inj1
## AAAAAAAAAAGAAGGR                         120130000                1668800
## AAAAAAAAAAGDSDSWDADTFSMEDPVRK             14499000                 437460
## AAAAAAAGGAAIAVSTGIETATIQK                        0                 173070
## AAAAAAGAASGIPGPVAQGIK                     14162000                 520550
## AAAAAAGPEMVR                              70512000                1604600
## AAAAAAPGGGGGEPR                           13873000                 361240
##                               Intensity.Treg.M6_inj1
## AAAAAAAAAAGAAGGR                                   0
## AAAAAAAAAAGDSDSWDADTFSMEDPVRK                3881000
## AAAAAAAGGAAIAVSTGIETATIQK                     317830
## AAAAAAGAASGIPGPVAQGIK                        1608500
## AAAAAAGPEMVR                                 4689900
## AAAAAAPGGGGGEPR                              1168900
pepData %>% exprs %>% head
##                               Intensity.Tconv.M12_2 Intensity.Tconv.M12_3
## AAAAAAAAAAGAAGGR                          111470000             101350000
## AAAAAAAAAAGDSDSWDADTFSMEDPVRK              13303000              10160000
## AAAAAAAGGAAIAVSTGIETATIQK                         0                     0
## AAAAAAGAASGIPGPVAQGIK                             0               7197400
## AAAAAAGPEMVR                              104440000             101950000
## AAAAAAPGGGGGEPR                            11910000              12683000
##                               Intensity.Tconv.M5_inj1
## AAAAAAAAAAGAAGGR                               965070
## AAAAAAAAAAGDSDSWDADTFSMEDPVRK                  235500
## AAAAAAAGGAAIAVSTGIETATIQK                       48225
## AAAAAAGAASGIPGPVAQGIK                          266620
## AAAAAAGPEMVR                                   947110
## AAAAAAPGGGGGEPR                                190550
##                               Intensity.Tconv.M6_inj1 Intensity.Treg.M12_2
## AAAAAAAAAAGAAGGR                                    0             80606000
## AAAAAAAAAAGDSDSWDADTFSMEDPVRK                 6744500                    0
## AAAAAAAGGAAIAVSTGIETATIQK                      314120                    0
## AAAAAAGAASGIPGPVAQGIK                         2904800             12602000
## AAAAAAGPEMVR                                 12747000             44735000
## AAAAAAPGGGGGEPR                               2166100             12555000
##                               Intensity.Treg.M12_3 Intensity.Treg.M5_inj1
## AAAAAAAAAAGAAGGR                         120130000                1668800
## AAAAAAAAAAGDSDSWDADTFSMEDPVRK             14499000                 437460
## AAAAAAAGGAAIAVSTGIETATIQK                        0                 173070
## AAAAAAGAASGIPGPVAQGIK                     14162000                 520550
## AAAAAAGPEMVR                              70512000                1604600
## AAAAAAPGGGGGEPR                           13873000                 361240
##                               Intensity.Treg.M6_inj1
## AAAAAAAAAAGAAGGR                                   0
## AAAAAAAAAAGDSDSWDADTFSMEDPVRK                3881000
## AAAAAAAGGAAIAVSTGIETATIQK                     317830
## AAAAAAGAASGIPGPVAQGIK                        1608500
## AAAAAAGPEMVR                                 4689900
## AAAAAAPGGGGGEPR                              1168900
pepData %>% sampleNames
## [1] "Intensity.Tconv.M12_2"   "Intensity.Tconv.M12_3"  
## [3] "Intensity.Tconv.M5_inj1" "Intensity.Tconv.M6_inj1"
## [5] "Intensity.Treg.M12_2"    "Intensity.Treg.M12_3"   
## [7] "Intensity.Treg.M5_inj1"  "Intensity.Treg.M6_inj1"

The sample names are rather long and contain information on the spike-in concentration and the repeat. We will remove “Intensity.” from the string

sampleNames(pepData) <- pepData %>% sampleNames %>% str_replace(., pattern="Intensity.", replacement="")

Next we will add information on the proteins to the feature data. The feature data of the object pepData contains information for each peptide in the experiment. This info on the proteins can be found in the proteingroups.txt file of maxquant. The peptides.txt file contains many data on each feature (peptide). We will only retain specified columns Protein, Sequence, Reverse and Potential.Contaminant, indicating the protein, peptide sequence, if the peptide is a decoy and if it is a contaminant, respectively.

In older and more recent versions of maxquant the names of the fields might vary. We can therefore explore the names we want to select with the shiny app that can be launched from the selectFeatureData from the MSnbase package, i.e. by running the function without an fcol argument that specifies the columns you would like to select.

The code can be found below, but is not executed by this Rmd script (note, that the flag eval=FALSE has been used as option for the R chunk).

pepData<-selectFeatureData(pepData)

If we know the columns, we can specify them as a vector in the fcol slot.

pepData<-selectFeatureData(pepData,fcol=c("Proteins","Protein.names","Gene.names", "Potential.contaminant","Reverse","Sequence"))

Calculate how many non zero intensities we have per peptide. This will be useful for filtering.

fData(pepData)$nNonZero<- rowSums(exprs(pepData) > 0)

Peptides with zero intensities are missing peptides and should be represent with a NA value instead of 0.

exprs(pepData)[exprs(pepData)==0]<-NA

Next we create the data on the experimental layout. We can do this based on the samplenames. We will add an effect for cell type and also for the animal (mouse ID).

pData(pepData)$cellType <- pepData %>% sampleNames %>% substr(1,4) %>% as.factor
pData(pepData)$mouseID <- factor(paste0("m",rep(1:4,2)))

We will also replace the sample names to make them shorter

sampleNames(pepData)<-pepData %>% sampleNames %>% substr(1,4) %>% paste0(.,"_m",rep(1:4,2))

2.1 Data exploration

We can inspect the missingness in our data with the plotNA() function provided with MSnbase. 38% of all peptide intensities are missing and for some peptides we don’t even measure a signal in any sample. The missingness is similar across samples. Note, that we plot the peptide data, so the label protein in the plot refers to peptides.

plotNA(pepData)

When we look at density plots of the intensities, we see that most intensities are low and that there is a long tail to the right. It is more useful to explore the data upon log transformation.

plotDensities(exprs(pepData))

nrow(pepData)
## [1] 55814

There are 55814 peptides before preprocessing.

3 Preprocessing

We will log transform, normalize, filter and summarize the data.

3.1 Log transform the data

pepData <- log(pepData, base = 2)
plotDensities(exprs(pepData))

When we look at the densities we see small shifts in location and shape of the density distributions. So normalisation will be required. But, we will first perform filtering.

3.2 Filtering

3.2.1 Handling overlapping protein groups

In our approach a peptide can map to multiple proteins, as long as there is none of these proteins present in a smaller subgroup.

pepData<-pepData[fData(pepData)$Proteins %in% smallestUniqueGroups(fData(pepData)$Proteins),]

3.2.2 Remove reverse sequences (decoys) and contaminants

We now remove the contaminants, peptides that map to decoy sequences and proteins, which were only identified by peptides with modifications.

pepData <- pepData[fData(pepData)$Reverse!="+",]
pepData <- pepData[fData(pepData)$Potential.contaminant!="+",]

3.2.3 Remove peptides of proteins that were only identified with modified peptides

Proteins for which all peptides are carrying modifications (PTMs) can be considered as unreliable. We will filter out these proteins. This information is included in the Only.identified.by.site column of proteinGroups.txt file of maxquant. The variable is containing a “+” character if the protein is only identified by modified peptides and an empty string “” if this is not the case. Sometimes an NA value is also present. We will replace these by the empty character string. The information of the peptides.txt file can be linked to the proteinGroups.txt file by using the Proteins column from the peptides.txt file and the Protein.IDs column in the proteinGroups.txt file. If NAs are included in the

proteinGroups<-read.table(proteinGroupsFile, sep = "\t", header = TRUE, quote = "", comment.char = "")
only_site <- proteinGroups$Only.identified.by.site
only_site[is.na(only_site)] <- ""
proteinsOnlyIdentifiedWithPTMs <- proteinGroups$Protein.IDs[only_site=="+"]
pepData<-pepData[!(fData(pepData)$Proteins %in% proteinsOnlyIdentifiedWithPTMs),]

We will now remove the proteinGroups object from the R session to free memory.

rm(proteinGroups,only_site)
gc()
##           used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells 4158132 222.1    8603538 459.5         NA  5452430 291.2
## Vcells 8237828  62.9   24434645 186.5      16384 24430794 186.4

3.2.4 Drop peptides that were only identified in one sample

We want to keep peptide that were at least observed twice.

pepData<-pepData[fData(pepData)$nNonZero>=2,]
nrow(pepData)
## [1] 44396

We keep 44396 peptides upon filtering.

3.3 Quantile normalize the data

pepData <- normalise(pepData, "quantiles")

3.4 Explore quantile normalized data

Upon normalisation the density curves for all samples coincide.

plotDensities(exprs(pepData))

We can visualize our data using a Multi Dimensional Scaling plot, eg. as provided by the limma package.

plotMDS(exprs(pepData),col=as.double(pData(pepData)$cellType))

The first axis in the plot is showing the leading log fold changes (differences on the log scale) between the samples. We notice that the leading differences (log FC) in the peptide data seems to be driven by technical variability. Indeed the samples do not seem to be clearly separated according to the spike in condition. Because there are missing values for the peptide level data we use the na.rm=TRUE argument to summarize the data based on the observed peptides for every protein.

3.5 Summarization to protein level

protData<-combineFeatures(pepData,fcol="Proteins",method="robust",na.rm=TRUE)
## Your data contains missing values. Please read the relevant section
## in the combineFeatures manual page for details the effects of
## missing values on data aggregation.

If you summarize using robust summarization please refer to Sticker et al. (2019) https://www.biorxiv.org/content/10.1101/668863v1

We filter proteins for which we have at least two protein expression values in both conditions:

protData<-protData[
protData %>% exprs %>% is.na %>% `!` %>% rowSums>4,]
plotMDS(exprs(protData),col=as.double(pData(pepData)$cellType))

4 Data Analysis

4.1 Estimation

MSqRob is currently working with a format where we have one dataframe for each protein. This will be changed in the next release. Therefore we first have to reorganise the data.

Next the models are fitted. This is done using the fit.model function. We only have to model the data using the factor condition from the pheno data of the protein level MSnSet. The name of the factor variable is specified in the fixed argument (if multiple predictors have to be incorporated in the model, a vector of variable names has to be provided in this argument.). The argument shrinkage is used to specify if ridge regression has to be adopted. For the sake of speed we do not do this in the tutorial. The shrinkage has to be specified for each variable in the fixed effects. We also have to indicate this for the intercept (which we never shrink). So we specify it at c(0,0) to indicate that the intercept (first 0) and the parameters for the factor condition (second 0) are not penalized. We set the robust_var function equal to FALSE, this functionality will be removed from the package in the next release.

We analysis the data as a completely randomized design and only include cellType in the model.

protMSqRob <- MSnSet2protdata(protData, "Proteins",annotations=c("Protein.names","Gene.names"))
models <- fit.model(protdata=protMSqRob, response="quant_value", fixed=c("cellType"),shrinkage.fixed=c(0,0),robust_var=FALSE)

4.2 Inference

Many biologists have problems with the reference coding. In MSqRob we have opted to formulate contrasts using all levels of a factor. Internally, the contrasts are than recasted according to the factor level that is the reference class.

L <- makeContrast("cellTypeTreg - cellTypeTcon",  levels=c("cellTypeTreg","cellTypeTcon"))
tests <- test.contrast_adjust(models, L)
nSig <- sum(tests$signif,na.rm=TRUE)
head(tests,nSig)
##                                                                         Protein.names
## P24452                                                     Macrophage-capping protein
## P21550                                                                   Beta-enolase
## P07356                                                                     Annexin A2
## Q9DC16                 Endoplasmic reticulum-Golgi intermediate compartment protein 1
## Q61503                                                                 5-nucleotidase
## Q3UW53                                                                  Protein Niban
## Q80SU7                                         Interferon-induced very large GTPase 1
## Q60611                                                      DNA-binding protein SATB1
## Q61205                    Platelet-activating factor acetylhydrolase IB subunit gamma
## Q9QXY6                                                 EH domain-containing protein 3
## O09131                                              Glutathione S-transferase omega-1
## P54071                                 Isocitrate dehydrogenase [NADP], mitochondrial
## Q7TPR4                                                                Alpha-actinin-1
## Q3UN02                                              Lysocardiolipin acyltransferase 1
## Q04447                                                         Creatine kinase B-type
## Q03267                                                     DNA-binding protein Ikaros
## P13020                                                                       Gelsolin
## P09055                                                                Integrin beta-1
## Q9CQ62                                       2,4-dienoyl-CoA reductase, mitochondrial
## Q99JB6                                                        Forkhead box protein P3
## P57016                                                                      Ladinin-1
## P81183                                                     Zinc finger protein Helios
## P08207                                                               Protein S100-A10
## Q8BGW0                                                                 Protein THEMIS
## P00329                                                        Alcohol dehydrogenase 1
## Q8CFB4                                                    Guanylate-binding protein 5
## Q7TMB8                                         Cytoplasmic FMR1-interacting protein 1
## Q05186                                                               Reticulocalbin-1
## Q60710                        Deoxynucleoside triphosphate triphosphohydrolase SAMHD1
## P49718                                          DNA replication licensing factor MCM5
## P07742                             Ribonucleoside-diphosphate reductase large subunit
## O70400                                                   PDZ and LIM domain protein 1
## P70271                                                   PDZ and LIM domain protein 4
## Q61107                                                    Guanylate-binding protein 4
## O88508                                          DNA (cytosine-5)-methyltransferase 3A
## Q00417                                                         Transcription factor 7
## Q9QYB5                                                                  Gamma-adducin
## Q8VCT3                                                               Aminopeptidase B
## P70302                                                 Stromal interaction molecule 1
## Q9WTK5   Nuclear factor NF-kappa-B p100 subunit;Nuclear factor NF-kappa-B p52 subunit
## P70677                          Caspase-3;Caspase-3 subunit p17;Caspase-3 subunit p12
## Q8BFP9 [Pyruvate dehydrogenase (acetyl-transferring)] kinase isozyme 1, mitochondrial
## O70370                                                                    Cathepsin S
## P48758                                                   Carbonyl reductase [NADPH] 1
## Q8VCW8                             Acyl-CoA synthetase family member 2, mitochondrial
## O70293                                            G protein-coupled receptor kinase 6
## Q9Z0F6                                    Cell cycle checkpoint control protein RAD9A
## Q9CYL5                          Golgi-associated plant pathogenesis-related protein 1
## Q9WUZ9                               Ectonucleoside triphosphate diphosphohydrolase 5
## P54310                                                       Hormone-sensitive lipase
## P10630                                             Eukaryotic initiation factor 4A-II
## P42230                            Signal transducer and activator of transcription 5A
## P61028                                                     Ras-related protein Rab-8B
## Q8K157                                                             Aldose 1-epimerase
## O88673                                                    Diacylglycerol kinase alpha
## Q6PHZ2              Calcium/calmodulin-dependent protein kinase type II subunit delta
## Q9R1Q7                                                          Proteolipid protein 2
## P47856               Glutamine--fructose-6-phosphate aminotransferase [isomerizing] 1
## Q3UGR5               Haloacid dehalogenase-like hydrolase domain-containing protein 2
## P97311                                          DNA replication licensing factor MCM6
## P25799   Nuclear factor NF-kappa-B p105 subunit;Nuclear factor NF-kappa-B p50 subunit
## Q68FL4                                              Putative adenosylhomocysteinase 3
## P29416                                              Beta-hexosaminidase subunit alpha
## P37913                                                                   DNA ligase 1
## Q9DBX6                                                            Cytochrome P450 2S1
## Q8C0P5                                                                     Coronin-2A
## Q9D8S4                                               Oligoribonuclease, mitochondrial
## Q3UDE2                                       Tubulin--tyrosine ligase-like protein 12
## Q99MN9                            Propionyl-CoA carboxylase beta chain, mitochondrial
##        Gene.names   estimate        se       df      Tval         pval
## P24452       Capg  2.8361630 0.1481032 8.738754 19.149907 1.937725e-08
## P21550       Eno3  2.6045120 0.1872245 8.738754 13.911171 2.915818e-07
## P07356      Anxa2  2.2241620 0.1849286 8.738754 12.027137 9.819877e-07
## Q9DC16     Ergic1  2.1805381 0.2002555 8.738754 10.888779 2.228514e-06
## Q61503       Nt5e  2.5718806 0.2700720 8.738754  9.522944 6.617576e-06
## Q3UW53    Fam129a  3.0426788 0.3267307 8.738754  9.312497 7.917717e-06
## Q80SU7      Gvin1  1.5269534 0.1692601 8.738754  9.021343 1.020549e-05
## Q60611      Satb1 -1.8750336 0.2088208 8.738754 -8.979151 1.059379e-05
## Q61205   Pafah1b3  1.2694962 0.1422006 8.738754  8.927504 1.109148e-05
## Q9QXY6       Ehd3 -1.2786926 0.1448946 8.738754 -8.824985 1.215762e-05
## O09131      Gsto1  3.0611645 0.3214327 7.738754  9.523502 1.527546e-05
## P54071       Idh2 -1.3263895 0.1583619 8.738754 -8.375683 1.837007e-05
## Q7TPR4      Actn1 -1.6049906 0.1922279 8.738754 -8.349414 1.882911e-05
## Q3UN02     Lclat1  2.1579100 0.2554461 7.738754  8.447615 3.589713e-05
## Q04447        Ckb -1.6112296 0.2181731 8.738754 -7.385097 4.875052e-05
## Q03267      Ikzf1 -1.0968334 0.1513320 8.738754 -7.247863 5.624519e-05
## P13020        Gsn -1.6477614 0.2299394 8.738754 -7.166069 6.130837e-05
## P09055      Itgb1  1.4669577 0.1876525 7.738754  7.817418 6.180708e-05
## Q9CQ62      Decr1  1.4413715 0.2048099 8.738754  7.037607 7.029981e-05
## Q99JB6      Foxp3  2.5400776 0.2584505 5.738754  9.828103 8.371882e-05
## P57016       Lad1  2.4836822 0.3393424 7.738754  7.319105 9.742648e-05
## P81183      Ikzf2  2.5985811 0.2774924 5.738754  9.364512 1.089134e-04
## P08207    S100a10  2.6645861 0.4031672 8.738754  6.609134 1.124747e-04
## Q8BGW0     Themis -2.0237546 0.3070787 8.738754 -6.590344 1.148729e-04
## P00329       Adh1 -1.4440678 0.2214092 8.738754 -6.522166 1.240567e-04
## Q8CFB4       Gbp5  1.5306967 0.2376108 8.738754  6.442034 1.358909e-04
## Q7TMB8     Cyfip1  2.3860884 0.3087404 6.738754  7.728461 1.380601e-04
## Q05186       Rcn1  1.7849845 0.2611635 7.738754  6.834740 1.552554e-04
## Q60710     Samhd1  1.0397824 0.1660315 8.738754  6.262561 1.671310e-04
## P49718       Mcm5  0.9113382 0.1482407 8.738754  6.147691 1.912018e-04
## P07742       Rrm1  1.1504309 0.1877895 8.738754  6.126174 1.961186e-04
## O70400     Pdlim1 -1.2622753 0.2082099 8.738754 -6.062514 2.114913e-04
## P70271     Pdlim4 -3.1346604 0.4821311 7.738754 -6.501677 2.170376e-04
## Q61107       Gbp4  1.9450820 0.3220050 8.738754  6.040534 2.171007e-04
## O88508     Dnmt3a  1.1899880 0.1993696 8.738754  5.968753 2.365794e-04
## Q00417       Tcf7 -1.5831026 0.2688529 8.738754 -5.888360 2.606851e-04
## Q9QYB5       Add3 -0.9784364 0.1674141 8.738754 -5.844409 2.749858e-04
## Q8VCT3      Rnpep  0.8356827 0.1436498 8.738754  5.817499 2.841620e-04
## P70302      Stim1 -0.9272306 0.1605023 8.738754 -5.777055 2.985867e-04
## Q9WTK5      Nfkb2  1.0016477 0.1750723 8.738754  5.721338 3.197796e-04
## P70677      Casp3  1.4794873 0.2598880 8.738754  5.692789 3.312689e-04
## Q8BFP9       Pdk1 -1.5355360 0.2708411 8.738754 -5.669509 3.409705e-04
## O70370       Ctss  1.0142563 0.1804486 8.738754  5.620751 3.623059e-04
## P48758       Cbr1 -0.9470080 0.1689381 8.738754 -5.605651 3.692044e-04
## Q8VCW8      Acsf2 -1.0588995 0.1908312 8.738754 -5.548881 3.964466e-04
## O70293       Grk6 -0.8791346 0.1585372 8.738754 -5.545288 3.982427e-04
## Q9Z0F6      Rad9a -2.1001814 0.3253428 6.738754 -6.455288 4.083312e-04
## Q9CYL5     Glipr2  1.1974930 0.2210053 8.738754  5.418392 4.677162e-04
## Q9WUZ9     Entpd5 -1.1587764 0.2006132 7.738754 -5.776172 4.709142e-04
## P54310       Lipe  1.1598976 0.1847450 6.738754  6.278371 4.808805e-04
## P10630     Eif4a2 -0.9820406 0.1833182 8.738754 -5.357027 5.059427e-04
## P42230     Stat5a  0.9179775 0.1713736 8.738754  5.356588 5.062283e-04
## P61028      Rab8b  0.9839867 0.1838247 8.738754  5.352853 5.086631e-04
## Q8K157       Galm  1.3934667 0.2604875 8.738754  5.349457 5.108888e-04
## O88673       Dgka -0.8249853 0.1542696 8.738754 -5.347684 5.120542e-04
## Q6PHZ2     Camk2d  0.8123895 0.1548639 8.738754  5.245827 5.841353e-04
## Q9R1Q7       Plp2  1.6357893 0.3149087 8.738754  5.194487 6.245749e-04
## P47856      Gfpt1  0.6693069 0.1292158 8.738754  5.179759 6.367286e-04
## Q3UGR5      Hdhd2 -0.8213089 0.1585940 8.738754 -5.178689 6.376215e-04
## P97311       Mcm6  0.7436700 0.1444667 8.738754  5.147692 6.640853e-04
## P25799      Nfkb1  0.6876488 0.1348708 8.738754  5.098575 7.084850e-04
## Q68FL4     Ahcyl2  0.9590399 0.1887357 8.738754  5.081391 7.247690e-04
## P29416       Hexa -0.8075849 0.1595596 8.738754 -5.061338 7.442857e-04
## P37913       Lig1  1.4654814 0.2905106 8.738754  5.044503 7.611093e-04
## Q9DBX6     Cyp2s1 -1.8538148 0.3216508 6.738754 -5.763440 7.897876e-04
## Q8C0P5     Coro2a  3.1046432 0.5402204 6.738754  5.746994 8.028159e-04
## Q9D8S4      Rexo2  0.8160401 0.1638234 8.738754  4.981219 8.281250e-04
## Q3UDE2     Ttll12 -0.8302515 0.1675411 8.738754 -4.955509 8.571497e-04
## Q99MN9       Pccb -0.7381137 0.1496424 8.738754 -4.932518 8.840369e-04
##                qval signif
## P24452 7.485433e-05   TRUE
## P21550 5.631903e-04   TRUE
## P07356 1.264473e-03   TRUE
## Q9DC16 2.152188e-03   TRUE
## Q61503 4.696490e-03   TRUE
## Q3UW53 4.696490e-03   TRUE
## Q80SU7 4.696490e-03   TRUE
## Q60611 4.696490e-03   TRUE
## Q61205 4.696490e-03   TRUE
## Q9QXY6 4.696490e-03   TRUE
## O09131 5.364465e-03   TRUE
## P54071 5.595142e-03   TRUE
## Q7TPR4 5.595142e-03   TRUE
## Q3UN02 9.905045e-03   TRUE
## Q04447 1.255489e-02   TRUE
## Q03267 1.326449e-02   TRUE
## P13020 1.326449e-02   TRUE
## P09055 1.326449e-02   TRUE
## Q9CQ62 1.429306e-02   TRUE
## Q99JB6 1.617029e-02   TRUE
## P57016 1.792183e-02   TRUE
## P81183 1.848976e-02   TRUE
## P08207 1.848976e-02   TRUE
## Q8BGW0 1.848976e-02   TRUE
## P00329 1.916924e-02   TRUE
## Q8CFB4 1.975282e-02   TRUE
## Q7TMB8 1.975282e-02   TRUE
## Q05186 2.141970e-02   TRUE
## Q60710 2.226300e-02   TRUE
## P49718 2.443891e-02   TRUE
## P07742 2.443891e-02   TRUE
## O70400 2.466647e-02   TRUE
## P70271 2.466647e-02   TRUE
## Q61107 2.466647e-02   TRUE
## O88508 2.611161e-02   TRUE
## Q00417 2.797296e-02   TRUE
## Q9QYB5 2.871000e-02   TRUE
## Q8VCT3 2.888731e-02   TRUE
## P70302 2.957540e-02   TRUE
## Q9WTK5 3.088272e-02   TRUE
## P70677 3.121200e-02   TRUE
## Q8BFP9 3.136117e-02   TRUE
## O70370 3.241447e-02   TRUE
## P48758 3.241447e-02   TRUE
## Q8VCW8 3.344373e-02   TRUE
## O70293 3.344373e-02   TRUE
## Q9Z0F6 3.356135e-02   TRUE
## Q9CYL5 3.596482e-02   TRUE
## Q9WUZ9 3.596482e-02   TRUE
## P54310 3.596482e-02   TRUE
## P10630 3.596482e-02   TRUE
## P42230 3.596482e-02   TRUE
## P61028 3.596482e-02   TRUE
## Q8K157 3.596482e-02   TRUE
## O88673 3.596482e-02   TRUE
## Q6PHZ2 4.029490e-02   TRUE
## Q9R1Q7 4.174800e-02   TRUE
## P47856 4.174800e-02   TRUE
## Q3UGR5 4.174800e-02   TRUE
## P97311 4.275603e-02   TRUE
## P25799 4.486685e-02   TRUE
## Q68FL4 4.515778e-02   TRUE
## P29416 4.563771e-02   TRUE
## P37913 4.594008e-02   TRUE
## Q9DBX6 4.693769e-02   TRUE
## Q8C0P5 4.698906e-02   TRUE
## Q9D8S4 4.774697e-02   TRUE
## Q3UDE2 4.869367e-02   TRUE
## Q99MN9 4.949325e-02   TRUE

There are 69 proteins with a significant effect at the 5% FDR level.

4.3 Plots

4.3.1 Volcano-plot

Note, that we find less proteins.

volcano<- ggplot(tests,aes(x=estimate,y=-log10(pval),color=signif)) + geom_point() + scale_color_manual(values=c("black","red"))
volcano

An interactive volcano plot can be obtained by using the plotly library. We alter the object tests to add the protein name to the data frame in ggplot, this will allow us to see the protein name in the interactive plot.

library(plotly)
volcano<- ggplot(tests %>% rownames_to_column("protein"), aes(x=estimate,y=-log10(pval),color=signif,protein=protein,annotation=Protein.names,gene=Gene.names)) + geom_point() + scale_color_manual(values=c("black","red"))
ggplotly(volcano)

4.3.2 Heatmap

There are no significant proteins so we will plot a heatmap of the top 10 proteins, instead. Note, that these results are not statistically significant.

sigNames<- tests %>% rownames_to_column("protein") %>% filter(signif) %>% pull(protein)
heatmap(exprs(protData)[sigNames,])

4.3.3 Detail plot for the top 5 proteins.

We first extract the normalized protein expression values for a particular protein and we create a dataframe that contains the protein expression values and the relevant variables of the experimental design.

for (protName in rownames(tests)[1:5])
{
plotDat<-pepData[fData(pepData)$Proteins==protName] %>% exprs
plotDataStack<-data.frame(quant_value=c(plotDat),sequence=rep(rownames(plotDat),rep=ncol(plotDat)),sample=rep(colnames(plotDat),each=nrow(plotDat)),condition=rep(pData(pepData)$cellType,each=nrow(plotDat)))
plot1 <- ggplot(plotDataStack, aes(x=sample, y=quant_value,fill=condition))
print(plot1 +geom_boxplot(outlier.shape=NA) + geom_point(position=position_jitter(width=.1),aes(shape=sequence)) + scale_shape_manual(values=1:nrow(plotDat)) +labs(title = protName, x="sample", y="Peptide intensity (log2)"))
}