Protein Identification and Initial Quantification
Multiple algorithms and software implementations have been developed for quantitation label-free proteomics experiments, in particular for extracted ion chromatograms (XIC). For background information you may look at Wikipedia labell-free Proteomics. Here, the use of the output for 3 such implementations for extracting peptide/protein quantifications is shown. These 3 software implementations were run individually using equivalent settings, ie identifcation based on the same fasta-database, starting at a single peptide with 1% FDR, MS mass tolerance for ion precursors at 0.7 ppm, oxidation of methionins and N-terminal acetylation as fixed as well as carbamidomethylation of cysteins as variable modifications.
Since in this context it is crucial to recognize all UPS1 proteins as such, the import-functions make use of the specPref argument, allowing to define custom tags. Most additional arguments to the various import-functions have been kept common for conventient use and for generating output structured the same way. Indeed, simply separating proteins by their species origin is not sufficient since common contaminants like human keratin might get considered by error as UPS1.
MaxQuant
MaxQuant is free software provided by the Max-Planck-Institute, see also Tyanova et al 2016. Later in this document data from MaxQuant will by frequently abbreviated as MQ.
Typically MaxQuant exports quantitation data on level of consensus-proteins by default to a folder called txt with a file called “proteinGroups.txt” . So in a standard case (when the file name has not been changed manually) it is sufficient to provide the path to this file. Of course, you can explicitely point to a specific file, as shown below. With the data presented here MaxQuant version 1.6.10 was run. Files compressed as .gz can be read, too (like in the example below).
path1 <- system.file("extdata", package="wrProteo")
fiNaMQ <- "proteinGroups.txt.gz"
specPrefMQ <- list(conta="CON_|LYSC_CHICK", mainSpecies="OS=Saccharomyces cerevisiae", spike=UPS1$ac)
dataMQ <- readMaxQuantFile(path1, file=fiNaMQ, specPref=specPrefMQ, refLi="mainSpe")
#> -> readMaxQuantFile : Suspicious filename, this function was designed for reading tabulated text files produced by MaxQuant
#> -> readMaxQuantFile : Note: Found 11 proteins marked as 'REV_' (reverse peptide identification) - Removing
#> -> readMaxQuantFile : Display 2845(9.5%) initial '0' values as 'NA'
#> -> readMaxQuantFile : Could not find peptide counts columns (argument 'pepCountCol') matching to 'Unique peptides MS\.MS\.count '
#> -> readMaxQuantFile : Note: 6 proteins with unknown species
#> by species : Homo sapiens: 49, Mus musculus: 1, Saccharomyces cerevisiae: 1047, Sus scrofa: 1,
#> -> readMaxQuantFile : Found 75 composite accession-numbers, truncating (eg P00761;CON__P00761)
#> -> readMaxQuantFile : Use column 'Accession' as identifyer (has fewest, ie 0 duplicated entries) as rownames
#> -> readMaxQuantFile : Normalize using subset of 1047
#> -> readMaxQuantFile -> vioplotW -> asSepList : argument 'fxArg' is depreciated, please use argument 'exclElem' (still using since 'exclElem=NULL'
#> -> readMaxQuantFile -> vioplotW -> asSepList : argument 'fxArg' is depreciated, please use argument 'exclElem' (still using since 'exclElem=NULL'
The data were imported and median-normalized, the protein annotation was parsed to automatically extract IDs, protein-names and species information.
## the number of lines and colums
dim(dataMQ$quant)
#> [1] 1104 27
## a summary of the quantitation data
summary(dataMQ$quant[,1:8]) # the first 8 cols
#> 12500amol_1 12500amol_2 12500amol_3 125amol_1
#> Min. :17.52 Min. :15.64 Min. :14.90 Min. :15.19
#> 1st Qu.:22.50 1st Qu.:22.49 1st Qu.:22.50 1st Qu.:22.39
#> Median :23.46 Median :23.46 Median :23.46 Median :23.43
#> Mean :23.69 Mean :23.65 Mean :23.67 Mean :23.60
#> 3rd Qu.:24.82 3rd Qu.:24.76 3rd Qu.:24.78 3rd Qu.:24.82
#> Max. :30.31 Max. :30.29 Max. :30.35 Max. :30.25
#> NA's :98 NA's :100 NA's :104 NA's :114
#> 125amol_2 125amol_3 25000amol_1 25000amol_2
#> Min. :14.86 Min. :14.91 Min. :15.78 Min. :18.42
#> 1st Qu.:22.36 1st Qu.:22.40 1st Qu.:22.52 1st Qu.:22.54
#> Median :23.42 Median :23.44 Median :23.53 Median :23.53
#> Mean :23.59 Mean :23.62 Mean :23.74 Mean :23.73
#> 3rd Qu.:24.81 3rd Qu.:24.79 3rd Qu.:24.95 3rd Qu.:24.89
#> Max. :30.25 Max. :30.30 Max. :30.31 Max. :30.13
#> NA's :106 NA's :114 NA's :109 NA's :115Now we can summarize the presence of UPS1 proteins after treatment by MaxQuant : In sum 48 UPS1 proteins were found, 0 are missing.
ProteomeDiscoverer
ProteomeDiscoverer is commercial software from ThermoFisher (www.thermofisher.com). Later in this document data from ProteomeDiscoverer will by frequently abbreviated as PD.
With the data used here, the identification was performed using the XCalibur module of ProteomeDiscoverer version 2.4 . Quantitation data at the level of consensus-proteins can be exported to tabulated text files, which can be treated by the function shown below. The resultant data were export in tablulated format and the file automatically named ‘_Proteins.txt’ (the option R-headers was checked, however this option is not mandatory). Files compressed as .gz can be read, too (like in the example below).
path1 <- system.file("extdata", package="wrProteo")
fiNaPd <- "pxd001819_PD24_Proteins.txt.gz"
## Note: data exported from ProteomeDiscoverer frequently do not have proper column-names, providing names here
sampNa <- paste(rep(UPSconc, each=3),"amol_",rep(1:3,length(UPSconc)),sep="")
specPrefPD <- list(conta="Bos tauris|Gallus", mainSpecies="OS=Saccharomyces cerevisiae", spike=UPS1$ac)
dataPD <- readProtDiscovFile(file=fiNaPd, path=path1, sampleNames=sampNa, refLi="mainSpe", specPref=specPrefPD)
#> -> readProtDiscovFile : setting 'annotCol' to export of 'R-friendly' colnames
#> -> readProtDiscovFile : Note: 10 (out of 1297) lines with unrecognized species
#> -> readProtDiscovFile : Count by 'specPref' : Gallus gallus: 1 ; Homo sapiens: 47 ; Saccharomyces cerevisiae: 1239 ;
#> -> readProtDiscovFile : Removing 1 lines due to duplicated Accessions (typically due to contaminants)
#> -> readProtDiscovFile : Normalize using (custom) subset of 1239 lines
#> -> readProtDiscovFile -> vioplotW -> asSepList : argument 'fxArg' is depreciated, please use argument 'exclElem' (still using since 'exclElem=NULL'
#> -> readProtDiscovFile -> vioplotW -> asSepList : argument 'fxArg' is depreciated, please use argument 'exclElem' (still using since 'exclElem=NULL'
The data were imported and median-normalized, the protein annotation was parsed to atomatically extract IDs, protein-names and species information.
Note, that the original column-heads in the file exported from ProteomeDiscoverer has simply increasing numbers as names. Thus, much care should be taken on the order when preparing the vector sampleNames with the names to use instead.
## the number of lines and colums
dim(dataPD$quant)
#> [1] 1296 27
## a summary of the quantitation data
summary(dataPD$quant[,1:8]) # the first 8 cols
#> 50amol_1 50amol_2 50amol_3 125amol_1
#> Min. :13.68 Min. :11.60 Min. : 9.609 Min. :10.92
#> 1st Qu.:18.15 1st Qu.:18.19 1st Qu.:18.237 1st Qu.:18.23
#> Median :19.36 Median :19.37 Median :19.380 Median :19.37
#> Mean :19.52 Mean :19.53 Mean :19.539 Mean :19.52
#> 3rd Qu.:20.83 3rd Qu.:20.88 3rd Qu.:20.895 3rd Qu.:20.84
#> Max. :26.33 Max. :26.38 Max. :26.418 Max. :26.41
#> NA's :92 NA's :93 NA's :87 NA's :88
#> 125amol_2 125amol_3 250amol_1 250amol_2
#> Min. :10.51 Min. :11.15 Min. :10.04 Min. :10.21
#> 1st Qu.:18.21 1st Qu.:18.20 1st Qu.:18.15 1st Qu.:18.17
#> Median :19.40 Median :19.40 Median :19.40 Median :19.37
#> Mean :19.51 Mean :19.53 Mean :19.50 Mean :19.47
#> 3rd Qu.:20.82 3rd Qu.:20.83 3rd Qu.:20.78 3rd Qu.:20.78
#> Max. :26.37 Max. :26.47 Max. :26.50 Max. :26.45
#> NA's :99 NA's :86 NA's :87 NA's :77
# there are proteins where the 'OS='-tag won't be visible as Species (orig Fasta-header and Protein-name not accessible) :
which(is.na(dataPD$annot[,"Species"]) & dataPD$annot[,"SpecType"]=="species2") # is NA as Species
#> P02768
#> 25Confirming the presence of UPS1 proteins by ProteomeDiscoverer:
Now we can summarize the presence of UPS1 proteins after treatment by ProteomeDiscoverer : In sum 48 UPS1 proteins were found, 0 are missing.
Proline
Proline is open-source software provided by the Profi-consortium (see also proline-core on github), published by Bouyssié et al 2020. Later in this document data from Proline will by frequently abbreviated as PL.
Protein identification in Proline gets performed by SearchGUI, see also Vaudel et al 2015. In this case X!Tandem (see also Duncan et al 2005) was used as search engine.
Quantitation data at the level of consensus-proteins can be exported from Proline as .xlsx or tabulated text files, both formats can be treated by the import-functions shown below. Here, Proline version 1.6.1 was used.
path1 <- system.file("extdata", package="wrProteo")
fiNaPl <- "pxd001819_PL.xlsx"
specPrefPL <- c(conta="_conta", mainSpecies="OS=Saccharomyces cerevisiae", spike="_ups")
dataPL <- readProlineFile(file.path(path1,fiNaPl), specPref=specPrefPL, normalizeMeth="median", refLi="mainSpe", silent=TRUE)
#> -> readProlineFile -> vioplotW -> asSepList : argument 'fxArg' is depreciated, please use argument 'exclElem' (still using since 'exclElem=NULL'
#> -> readProlineFile -> vioplotW -> asSepList : argument 'fxArg' is depreciated, please use argument 'exclElem' (still using since 'exclElem=NULL'
In addition, we need to correct the quantification column-heads (like ‘Levure2ug+ UPS1-100amol-1’) and bring them to a simpler version :
head(colnames(dataPL$raw),8)
#> [1] "Levure2ug+ UPS1-100amol-1" "Levure2ug+ UPS1-100amol-2"
#> [3] "Levure2ug+ UPS1-100amol-3" "Levure2ug+ UPS1-250amol-1"
#> [5] "Levure2ug+ UPS1-250amol-2" "Levure2ug+ UPS1-250amol-3"
#> [7] "Levure2ug+ UPS1-500amol-1" "Levure2ug+ UPS1-500amol-2"
sub1 <- cbind(ini=paste0("-",c("100a","250a","500a","1f","5f","10f","25f","50f","100f"),"mol-"),
paste0(fin=c("50a","125a","250a","500a","2500a","5000a","12500a","25000a","50000a"),"mol_"))
dataPL <- cleanListCoNames(dataPL, rem=c("abundance_","Levure2ug+ UPS1"), subst=sub1)
## let's check the result
head(colnames(dataPL$raw),8)
#> [1] "50amol_1" "50amol_2" "50amol_3" "125amol_1" "125amol_2" "125amol_3"
#> [7] "250amol_1" "250amol_2"## the number of lines and colums
dim(dataPL$quant)
#> [1] 1186 27
## a summary of the quantitation data
summary(dataPL$quant[,1:8]) # the first 8 cols
#> 50amol_1 50amol_2 50amol_3 125amol_1
#> Min. :14.29 Min. :14.53 Min. :13.77 Min. :14.41
#> 1st Qu.:19.64 1st Qu.:19.66 1st Qu.:19.70 1st Qu.:19.75
#> Median :21.64 Median :21.55 Median :21.65 Median :21.65
#> Mean :21.66 Mean :21.66 Mean :21.71 Mean :21.75
#> 3rd Qu.:23.64 3rd Qu.:23.61 3rd Qu.:23.63 3rd Qu.:23.70
#> Max. :29.52 Max. :29.42 Max. :29.44 Max. :29.50
#> NA's :44 NA's :40 NA's :40 NA's :46
#> 125amol_2 125amol_3 250amol_1 250amol_2
#> Min. :13.53 Min. :14.76 Min. :14.55 Min. :14.74
#> 1st Qu.:19.76 1st Qu.:19.76 1st Qu.:19.74 1st Qu.:19.72
#> Median :21.64 Median :21.65 Median :21.65 Median :21.61
#> Mean :21.73 Mean :21.73 Mean :21.71 Mean :21.69
#> 3rd Qu.:23.69 3rd Qu.:23.65 3rd Qu.:23.62 3rd Qu.:23.63
#> Max. :29.49 Max. :29.43 Max. :29.42 Max. :29.40
#> NA's :46 NA's :48 NA's :54 NA's :51Now we can summarize the presence of UPS1 proteins after treatment by Proline : In sum 48 UPS1 proteins were found, 0 are missing.
Uniform Re-Arranging of Data
For easy and proper comparisons we need to make sure all columns are in the same order.
# get all results (MaxQuant,ProteomeDiscoverer, ...) in same order
grp9 <- paste0(rep(UPSconc,each=3),"amol")
## it is more convenient to re-order columns this way in each project
dataPD <- corColumnOrder(dataPD, sampNames=sampNa) # already in good order
#> -> corColumnOrder : Order already correct
dataMQ <- corColumnOrder(dataMQ, sampNames=sampNa)
dataPL <- corColumnOrder(dataPL, sampNames=sampNa)
#> -> corColumnOrder : Order already correctAt import we made use of the argument specPref (specifying _‘mainSpecies’, ‘conta’ and ‘spike’) which allows to build categories based on searching keywords based on the initial annotation. In turn, we obtain the labels : ‘main Spe’ for yeast (ie matrix), ‘species2’ for the UPS1 (ie spike) ‘conta’ for contaminants. Let’s replace the first two generic terms by more specific ones (ie ‘Yeast’ and ‘UPS1’) :
## Need to rename $annot[,"SpecType"]
dataPD <- replSpecType(dataPD, replBy=cbind(old=c("mainSpe","species2"), new=c("Yeast","UPS1")))
#> useLi 1 2 3 4 5 6
#> useLi 24 25 36 37 45 70
dataMQ <- replSpecType(dataMQ, replBy=cbind(old=c("mainSpe","species2"), new=c("Yeast","UPS1")))
#> useLi 1 12 13 14 15 16
#> useLi 4 11 21 27 39 47
dataPL <- replSpecType(dataPL, replBy=cbind(old=c("mainSpe","species2"), new=c("Yeast","UPS1")))
#> useLi 1 2 3 4 5 6
#> useLi 20 22 48 53 54 71
## Need to address missing ProteinNames (UPS1) due to missing tags in Fasta
dataPD <- replMissingProtNames(dataPD)
#> -> replreplMissingProtNames : ..trying to replace 1296 'EntryName'
dataMQ <- replMissingProtNames(dataMQ)
#> -> replreplMissingProtNames : ..trying to replace 5 'EntryName'
dataPL <- replMissingProtNames(dataPL) ## extract names of quantified UPS1-proteins
NamesUpsPD <- dataPD$annot[which(dataPD$annot[,"SpecType"]=="UPS1"),"Accession"]
NamesUpsMQ <- dataMQ$annot[which(dataMQ$annot[,"SpecType"]=="UPS1"),"Accession"]
NamesUpsPL <- dataPL$annot[which(dataPL$annot[,"SpecType"]=="UPS1"),"Accession"]tabS <- mergeVectors(PD=table(dataPD$annot[,"SpecType"]), MQ=table(dataMQ$annot[,"SpecType"]), PL=table(dataPL$annot[,"SpecType"]))
tabT <- mergeVectors(PD=table(dataPD$annot[,"Species"]), MQ=table(dataMQ$annot[,"Species"]), PL=table(dataPL$annot[,"Species"]))
tabT[which(is.na(tabT))] <- 0
kable(cbind(tabS[,2:1], tabT), caption="Number of proteins identified, by custom tags, species and software")| Yeast | UPS1 | Gallus gallus | Homo sapiens | Mus musculus | Saccharomyces cerevisiae | Sus scrofa | |
|---|---|---|---|---|---|---|---|
| PD | 1239 | 48 | 1 | 47 | 0 | 1239 | 0 |
| MQ | 1047 | 48 | 0 | 49 | 1 | 1047 | 1 |
| PL | 1137 | 48 | 0 | 48 | 0 | 1137 | 1 |
The initial fasta file also contained the yeast strain number, this has been stripped off when using default parameters.
