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"
## We need to define the setup of species
specPrefMQ <- list(conta="CON_|LYSC_CHICK", mainSpecies="OS=Saccharomyces cerevisiae", spike=UPS1$ac)
dataMQ <- readMaxQuantFile(path1, file=fiNaMQ, specPref=specPrefMQ, refLi="mainSpe", sdrf="PXD001819", suplAnnotFile=TRUE)
#> -> readMaxQuantFile : Note: Found 11 out of 1115 proteins marked as 'REV_' (reverse peptide identification) - Removing
#> -> readMaxQuantFile : Transform 2845(9.5%) initial '0' values to 'NA'
#> -> readMaxQuantFile : Could not find peptide counts columns (argument 'pepCountCol') matching to 'Unique peptides MS\.MS\.count '
#> -> readMaxQuantFile : Found 1 species name(s) appearing inside other ones, assume as truncated (eg Saccharomyces cerevis)
#> -> readMaxQuantFile : Note: 5 proteins with unknown species
#> data by species : Gallus gallus: 1, Homo sapiens: 49, Mus musculus: 1, Saccharomyces cerevisiae: 1047, Sus scrofa: 1,
#> -> readMaxQuantFile : Found 70 composite accession-numbers (eg P00761;P00761), truncating
#> -> readMaxQuantFile : Use column 'Accession' as identifyer (has fewest, ie 0 duplicated entries) as rownames
#> -> readMaxQuantFile : Normalize using subset of 1047
#> -> readMaxQuantFile -> readSampleMetaData : summaryD exists, but unable to find file-names
The data were imported and median-normalized, the protein annotation was parsed to automatically extract IDs, protein-names and species information. Please note, that in the example above we directly added information about the experimental setup from the sdrf repository.
## The number of lines and colums
dim(dataMQ$quant)
#> [1] 1104 27
## A quick summary of some columns of quantitation data
summary(dataMQ$quant[,1:8]) # the first 8 cols
#> 12500amol_1 12500amol_2 12500amol_3 125amol_1
#> Min. :17.54 Min. :15.66 Min. :14.93 Min. :15.18
#> 1st Qu.:22.51 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.81 3rd Qu.:24.76 3rd Qu.:24.77 3rd Qu.:24.82
#> Max. :30.29 Max. :30.27 Max. :30.32 Max. :30.26
#> NA's :98 NA's :100 NA's :104 NA's :114
#> 125amol_2 125amol_3 25000amol_1 25000amol_2
#> Min. :14.85 Min. :14.93 Min. :15.82 Min. :18.37
#> 1st Qu.:22.36 1st Qu.:22.40 1st Qu.:22.53 1st Qu.:22.53
#> Median :23.42 Median :23.44 Median :23.53 Median :23.53
#> Mean :23.59 Mean :23.62 Mean :23.74 Mean :23.74
#> 3rd Qu.:24.81 3rd Qu.:24.79 3rd Qu.:24.94 3rd Qu.:24.90
#> Max. :30.26 Max. :30.29 Max. :30.27 Max. :30.20
#> 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"
## Next, we define the setup of species
specPrefPD <- list(conta="Bos tauris|Gallus", mainSpecies="OS=Saccharomyces cerevisiae", spike=UPS1$ac)
dataPD <- readProtDiscovFile(file=fiNaPd, path=path1, refLi="mainSpe", specPref=specPrefPD, sdrf="PXD001819")
#> -> 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 -> readSampleMetaData : Sucessfully adjusted order of sdrf to content of pxd001819_PD_InputFiles.txt.gz
The data were imported and median-normalized, the protein annotation was parsed to automatically extract IDs, protein-names and species information.
Please note, that quantitation data exported from ProteomeDiscoverer frequently have very generic column-names (increasing numbers). When calling the import-function they can be replaced by more meaningful names either using the argument sampNa, or from reading the default annotation in the file ‘InputFiles.txt’ or, finally, from the sdrf-annotation. In the example below both the default annotation as file ‘InputFiles.txt’ and sdrf annotation are available and will be integrated to object produced by the import-function.
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 quick summary of some columns of quantitation data
summary(dataPD$quant[,1:8]) # the first 8 cols
#> 50amol_R1 50amol_R2 50amol_R3 125amol_R1
#> Min. :13.62 Min. :11.56 Min. : 9.585 Min. :10.90
#> 1st Qu.:18.14 1st Qu.:18.18 1st Qu.:18.234 1st Qu.:18.22
#> 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.84 3rd Qu.:20.89 3rd Qu.:20.899 3rd Qu.:20.84
#> Max. :26.39 Max. :26.42 Max. :26.435 Max. :26.43
#> NA's :92 NA's :93 NA's :87 NA's :88
#> 125amol_R2 125amol_R3 250amol_R1 250amol_R2
#> Min. :10.48 Min. :11.16 Min. :10.09 Min. :10.23
#> 1st Qu.:18.20 1st Qu.:18.20 1st Qu.:18.16 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.77 3rd Qu.:20.77
#> Max. :26.39 Max. :26.46 Max. :26.46 Max. :26.43
#> 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", sdrf="PXD001819", silent=TRUE)
#> -> readProlineFile -> readSampleMetaData : summaryD exists, but unable to find file-names
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"
dataPL <- cleanListCoNames(dataPL, rem=c("abundance_","Levure2ug+ UPS1"))
## let's check the result
head(colnames(dataPL$raw),8)
#> [1] "_100amol-1" "_100amol-2" "_100amol-3" "_250amol-1" "_250amol-2"
#> [6] "_250amol-3" "_500amol-1" "_500amol-2"## The number of lines and colums
dim(dataPL$quant)
#> [1] 1186 27
## A quick summary of some columns of quantitation data
summary(dataPL$quant[,1:8]) # the first 8 cols
#> _100amol-1 _100amol-2 _100amol-3 _250amol-1
#> Min. :14.29 Min. :14.52 Min. :13.78 Min. :14.43
#> 1st Qu.:19.64 1st Qu.:19.66 1st Qu.:19.70 1st Qu.:19.76
#> 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.62 3rd Qu.:23.62 3rd Qu.:23.70
#> Max. :29.52 Max. :29.42 Max. :29.43 Max. :29.48
#> NA's :44 NA's :40 NA's :40 NA's :46
#> _250amol-2 _250amol-3 _500amol-1 _500amol-2
#> Min. :13.54 Min. :14.77 Min. :14.55 Min. :14.73
#> 1st Qu.:19.77 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.68 3rd Qu.:23.65 3rd Qu.:23.62 3rd Qu.:23.64
#> Max. :29.47 Max. :29.42 Max. :29.43 Max. :29.41
#> 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, however using different software, this is not immediately the case. The basic names of the groups have already been figured out using the sample meta-data, notably the sdfr.
## bring all results (MaxQuant,ProteomeDiscoverer, ...) in same ascending order
## as reference will use the order from ProteomeDiscoverer, it's output is already in a convenient order
sampNa <- colnames(dataPD$quant)
## it is more convenient to re-order columns this way in each project
dataPD <- corColumnOrder(dataPD, sampNames=sampNa) # already in good order
#> -> corColumnOrder : Quant/counting data already in good order ..
dataMQ <- corColumnOrder(dataMQ, sampNames=sampNa, newNames=sub("amol_", "amol_R", colnames(dataMQ$quant)))
dataPL <- corColumnOrder(dataPL, sampNames=sampNa)
#> -> corColumnOrder : Colnames of 'dat$quant' differ from 'sampNames', trying to adjust ..
#> -> corColumnOrder : Failed : Unable to match 27 suggested sampNames ( '50amol_R1', '50amol_R2', '50amol_R3', '125amol_R1', '125amol_R2', '125amol_R3', '250amol_R1', '250amol_R2', '250amol_R3', '500amol_R1', '500amol_R2', '500amol_R3', '2500amol_R1', '2500amol_R2', '2500amol_R3', '5000amol_R1', '5000amol_R2', '5000amol_R3', '12500amol_R1', '12500amol_R2', '12500amol_R3', '25000amol_R1', '25000amol_R2', '25000amol_R3', '50000amol_R1', '50000amol_R2' and '50000amol_R3' )At 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)
## add some grouping/sample annotation to PL
dataPL$sampleSetup$groups <- dataPD$sampleSetup$groups## 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 | 1 | 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.
Let’s extract the information of grouping of samples/replicates
grp9 <- dataMQ$sampleSetup$groups <- dataPD$sampleSetup$groups
head(grp9)
#> 50 50 50 125 125 125
#> 1 1 1 2 2 2
