Protein Identification and Initial Quantification
MaxQuant
MaxQuant is free software provided by the Max-Planck-Insutute, see Tyanova et al 2016. Typically MaxQuant exports by default quantitation data on level of consensus-proteins as a folder called txt with a file called proteinGroups.txt . So in a standard case one needs only to provide the path to this file.
path1 <- system.file("extdata", package="wrProteo")
fiNaMQ <- "proteinGroups.txt.gz"
specPrefMQ <- list(conta="CON_|LYSC_CHICK", mainSpecies="OS=Saccharomyces cerevisiae", spike=UPS1ac)
dataMQ <- readMaxQuantFile(path1, file=fiNaMQ, specPref=specPrefMQ, refLi="mainSpe")
#> -> readMaxQuantFile : Note: Found 11 proteins marked as 'REV_' (reverse peptide identification) - Removing
#> -> readMaxQuantFile : Note: 8 proteins with unknown species
#> by species : conta: 1 mainSpe: 1047 species2: 48
#> -> readMaxQuantFile : Found 1 species name(s) appearing inside other ones, assume as truncated (eg Saccharomyces cerevis)
#> -> readMaxQuantFile : Found 75 composite accession-numbers, truncating (eg P00761;CON__P00761)
#> -> readMaxQuantFile : Use column 'Accession' (has fewest, ie 0 duplicated entries) as rownames
#> -> readMaxQuantFile : normalize using subset of 1047
The data were imported and median-normalized, the protein annotation was parsed to automatically extract IDs, protein-names and species information.
In addition we need to correct the accesion-names for attached ‘ups’-tags (originating from the fasta-file) :
## 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 :115Confirming the presence of UPS1 proteins by MaxQuant: In sum 48 UPS1 proteins were found, 0 are missing.
ProteomeDiscoverer
ProteomeDiscoverer is commercial software from ThermoFisher (www.thermofisher.com). In this case, the identification was performed using the XCalibur module of ProteomeDiscoverer. In ProteomeDiscoverer quantitation data on level of consensus-proteins should be exported to tabulated text files, which can be treated by this function. The resultant data were export as ‘Proteins’ in tablulated format (the option R-headers was checked, however data can also be read when this option was not chosen).
path1 <- system.file("extdata", package="wrProteo")
fiNaPd <- "pxd001819_PD24_Proteins.txt.gz"
## Note: data exported from ProteomeDiscoverer do not have proper column-names, providing names here
UPSconc <- c(50,125,250,500,2500,5000,12500,25000,50000)
sampNa <- paste(rep(UPSconc, each=3),"amol_",rep(1:3,length(UPSconc)),sep="")
specPrefPD <- list(conta="Bos tauris|Gallus", mainSpecies="OS=Saccharomyces cerevisiae", spike=UPS1ac)
dataPD <- readPDExport(file=fiNaPd, path=path1, sampleNames=sampNa, refLi="mainSpe", specPref=specPrefPD)
#> -> readPDExport : Trouble ahead, expecting tabulated text file (this file might not be right format) !!
#> -> readPDExport : setting 'annotCol' to export of 'R-friendly' colnames
#> -> readPDExport : Note: 10 (out of 1297) unrecognized species
#> -> readPDExport : Count by 'specPref' : Gallus gallus: 1 ; Homo sapiens: 47 ; Saccharomyces cerevisiae: 1239 ;
#> -> readPDExport : normalize using subset of 1239
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, much care should be taken on the order when preparing the vector 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: In sum 48 UPS1 proteins were found, 0 are missing.
Proline
Proline is free software provided by the Profi-consortium, see Ramus et al 2016 and Bouyssié et al 2020 (Proline: an efficient and user-friendly software suite for large-scale proteomics. Bioinformatics. 2020, PMID: 32096818, DOI: 10.1093/bioinformatics/btaa118).
Underneith identifications in Proline are 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.
In Proline quantitation data at the level of consensus-proteins can be exported as xlsx or tabulated text files, both formats can be treated by the functions of this package.
## shifted for not printing
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")
#> -> readProlineFile : Found sheets 'Search settings and infos', 'Import and filters', 'Quant config' and 'Protein sets'
#> New names:
#> * `` -> ...2
#> -> readProlineFile : normalize using subset of 1137
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)## 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 :51Confirming the presence of UPS1 proteins by Proline : In sum 48 UPS1 proteins were found, 0 are missing.
Uniform Re-Arranging of Data
In order 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 correctThe from the protein annotation the membership to 3 groups was extracted : yeast (matrix) as “‘main Spe’”, UPS1 (spike) as ‘species2’ and other contaminants (‘conta’). The first two terms will be replace by more specific ones (‘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 addres missing ProteinNames (UPS1) due to missing tags in Fasta
dataPD <- replNAProtNames(dataPD)
#> ..replacing 1296 entry-names
dataMQ <- replNAProtNames(dataMQ)
#> ..replacing 1104 entry-names
dataPL <- replNAProtNames(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"]))
#> yy
knitr::kable(tabS, caption="Number of proteins identified, by custom tags and software")| UPS1 | Yeast | conta | |
|---|---|---|---|
| PD | 48 | 1239 | 1 |
| MQ | 48 | 1047 | 1 |
| PL | 48 | 1137 | 1 |
tabT <- mergeVectors(PD=table(dataPD$annot[,"Species"]), MQ=table(dataMQ$annot[,"Species"]), PL=table(dataPL$annot[,"Species"]))
#> yy
knitr::kable(tabT, caption="Number of proteins identified, by species and software")| Gallus gallus | Homo sapiens | Mus musculus | Saccharomyces cerevisiae | Saccharomyces cerevisiae (strain ATCC 204508 / S288c) | Sus scrofa | |
|---|---|---|---|---|---|---|
| PD | 1 | 47 | 0 | 1239 | 0 | 0 |
| MQ | 1 | 49 | 1 | 1047 | 0 | 1 |
| PL | 0 | 48 | 0 | 0 | 1137 | 1 |
Based on PCA one can see that the comparison with concentrations >= 250 aMol may actually be better to detect differences, as also confirmed by ROC part later.
On this PCA one can see a rather distinct group of low concentration-pairs (mostly containg a 50-250), low/med conc pairs (containing 2500, ev 5000) & rest (med+high to any)