Francisella tularensis dataset (design with pseudoreplication)
The bacterium Francisella tularensis is a facultative intracellular parasite of macrophages and must import host cell arginine in order to be able to multiply inside the host cell. Ramond et al. (2015) investigated the effect of gene deletion of a newly identified arginine transporter called ArgP. Their data is deposited in the PRIDE repository at http://www.ebi.ac.uk/pride/archive/projects/PXD001584. We used the authors’ supplied peptides.txt file which is the result of preprocessing with MaxQuant version 1.4.1.2. As the authors used now outdated GI numbers, we modified this file to incorporate NCBI accessions and protein names (September 25, 2016).
1. Importing data
MSqRob data input can come from any kind of source, as long as feature-level or peptide-level data is available. I will here focus on a MaxQuant example, but in fact, output from any search engine can be used.
1.1. Importing MaxQuant data
First, specify the location of the peptides.txt MaxQuant output file. By default, MaxQuant creates a “combined” folder in the same folder as the raw files are (or if the raw files are in different directories, in the directory of the first raw file). The peptides.txt file can then be found in “path_to_raw_files/combined/txt/peptides.txt”. Just provide the path to where the file is located on your computer. The system.file function is used here only to access the example tab-delimited peptides.txt file that is delivered with our package.
file_peptides_txt <- system.file("extdata/Francisella", "peptides.txt", package = "MSqRobData")
file_peptides_txt## [1] "/Library/Frameworks/R.framework/Versions/3.5/Resources/library/MSqRobData/extdata/Francisella/peptides.txt"Import the data using the import2MSnSet function. The pattern argument should contain a text string that is present in all column names that contain peptide intensities but not in the other column names. This defaults to "Intensity\ " for MaxQuant data. The default setting remove_pattern=TRUE then removes the pattern from these column names so that the names equal the names that were entered in MaxQuant’s “Experiment” column when conducting the search. In this particular Francisella experiment, our column names start with a number (either “1” or “3”), which is not a valid column name in R. Therefore, an “X” is prepended to these names.
peptidesFranc <- import2MSnSet(file_peptides_txt, filetype = "MaxQuant", remove_pattern = TRUE)import2MSnSet creates an MSnSet object with 3 slots:
exprs: contains a matrix of peptide intensities of which the rows correspond to the peptides present in the dataset and the columns correspond to the different mass spec runs.fData: contains a data frame containing all additional information (such as peptide sequences, protein names, scores, number of missed cleavages, etc.) present in the input file.pData: contains a data frame with 0 columns and as many rows as there are mass spec runs. This slot will be used to store the experimental annotation.
If you want to assess whether importing happened in a correct way, you can always check the exprs, fData and pData slots.
The exprs slot should contain all the peptide intensities grouped per column.
head(exprs(peptidesFranc))## X1WT_20_2h_n3_1 X1WT_20_2h_n3_2 X1WT_20_2h_n3_3 X1WT_20_2h_n4_1
## 1 0 0 9290400 92996000
## 2 0 20679000 0 0
## 3 28853000 28091000 30436000 19476000
## 4 339830000 420390000 393930000 235060000
## 5 240280000 0 241890000 222250000
## 6 0 31069000 31448000 0
## X1WT_20_2h_n4_2 X1WT_20_2h_n4_3 X1WT_20_2h_n5_1 X1WT_20_2h_n5_2
## 1 0 98059000 0 80803000
## 2 0 0 0 17358000
## 3 22050000 20687000 36588000 28654000
## 4 381090000 360270000 255710000 311370000
## 5 189930000 0 160440000 147230000
## 6 27721000 28657000 0 0
## X1WT_20_2h_n5_3 X3D8_20_2h_n3_1 X3D8_20_2h_n3_2 X3D8_20_2h_n3_3
## 1 0 0 95433000 0
## 2 13841000 12278000 10712000 0
## 3 7009100 14174000 8548000 7859200
## 4 271440000 225140000 289880000 282500000
## 5 143340000 135910000 185740000 168620000
## 6 0 0 0 0
## X3D8_20_2h_n4_1 X3D8_20_2h_n4_2 X3D8_20_2h_n4_3 X3D8_20_2h_n5_1
## 1 0 0 0 0
## 2 9500700 11111000 8723500 8125200
## 3 14871000 11690000 4470500 0
## 4 209490000 223610000 219700000 169120000
## 5 160240000 158000000 0 92542000
## 6 0 18253000 0 0
## X3D8_20_2h_n5_2 X3D8_20_2h_n5_3
## 1 0 0
## 2 0 0
## 3 4048100 3422300
## 4 178150000 170400000
## 5 102350000 100090000
## 6 0 0The fData slot contains all other columns.
head(fData(peptidesFranc))## Sequence Amino.acid.before First.amino.acid
## 1 AAAEELDTR K A
## 2 AAAGFVITASHNK R A
## 3 AAANEYELALAYSIEEVAPDLHK K A
## 4 AAANNPQLEAFK K A
## 5 AAASAGLVDEK K A
## 6 AAASLDLYSYPK K A
## Second.amino.acid Second.last.amino.acid Last.amino.acid
## 1 A T R
## 2 A N K
## 3 A H K
## 4 A F K
## 5 A E K
## 6 A P K
## Amino.acid.after A.Count R.Count N.Count D.Count C.Count Q.Count E.Count
## 1 K 3 1 0 1 0 0 2
## 2 F 4 0 1 0 0 0 0
## 3 Y 6 0 1 1 0 0 4
## 4 K 4 0 2 0 0 1 1
## 5 A 4 0 0 1 0 0 1
## 6 V 3 0 0 1 0 0 0
## G.Count H.Count I.Count L.Count K.Count M.Count F.Count P.Count S.Count
## 1 0 0 0 1 0 0 0 0 0
## 2 1 1 1 0 1 0 1 0 1
## 3 0 1 1 3 1 0 0 1 1
## 4 0 0 0 1 1 0 1 1 0
## 5 1 0 0 1 1 0 0 0 1
## 6 0 0 0 2 1 0 0 1 2
## T.Count W.Count Y.Count V.Count U.Count Length Missed.cleavages
## 1 1 0 0 0 0 9 0
## 2 1 0 0 1 0 13 0
## 3 0 0 2 1 0 23 0
## 4 0 0 0 0 0 12 0
## 5 0 0 0 1 0 11 0
## 6 0 0 2 0 0 12 0
## Mass Proteins GI.number Leading.razor.protein
## 1 974.4669 WP_003038655 gi|118497196 gi|118497196
## 2 1285.6779 WP_003041237 gi|118498194 gi|118498194
## 3 2516.2435 WP_003038915 gi|118497331 gi|118497331
## 4 1272.6463 WP_003038264 gi|118496879 gi|118496879
## 5 1030.5295 WP_003035781 gi|118497152 gi|118497152
## 6 1297.6554 WP_003039212 gi|118497492 gi|118497492
## Protein.names
## 1 hypothetical protein [Francisella tularensis]
## 2 phosphoglucosamine mutase [Francisella tularensis]
## 3 glycine--tRNA ligase subunit beta [Francisella tularensis]
## 4 molecular chaperone HtpG [Francisella tularensis]
## 5 delta-aminolevulinic acid dehydratase [Francisella tularensis]
## 6 phosphorylase [Francisella tularensis]
## Unique..Groups. Unique..Proteins. Charges PEP Score
## 1 yes yes 1,2 1.2531e-04 59.35
## 2 yes yes 2 5.3020e-05 45.825
## 3 yes yes 3 1.0517e-45 84.327
## 4 yes yes 2 8.7549e-24 108.56
## 5 yes yes 2 9.5062e-05 67.385
## 6 yes yes 2 1.2563e-02 31.675
## Experiment.1WT_20_2h_n3_1 Experiment.1WT_20_2h_n3_2
## 1 NA NA
## 2 NA 1
## 3 1 1
## 4 1 1
## 5 1 NA
## 6 NA 1
## Experiment.1WT_20_2h_n3_3 Experiment.1WT_20_2h_n4_1
## 1 1 1
## 2 NA NA
## 3 2 1
## 4 1 1
## 5 1 1
## 6 1 NA
## Experiment.1WT_20_2h_n4_2 Experiment.1WT_20_2h_n4_3
## 1 NA 1
## 2 NA NA
## 3 1 1
## 4 1 1
## 5 1 NA
## 6 1 1
## Experiment.1WT_20_2h_n5_1 Experiment.1WT_20_2h_n5_2
## 1 NA 1
## 2 NA 1
## 3 2 2
## 4 1 1
## 5 1 1
## 6 NA NA
## Experiment.1WT_20_2h_n5_3 Experiment.3D8_20_2h_n3_1
## 1 NA NA
## 2 1 1
## 3 1 1
## 4 1 1
## 5 1 1
## 6 NA NA
## Experiment.3D8_20_2h_n3_2 Experiment.3D8_20_2h_n3_3
## 1 1 NA
## 2 1 NA
## 3 2 1
## 4 1 1
## 5 1 1
## 6 NA NA
## Experiment.3D8_20_2h_n4_1 Experiment.3D8_20_2h_n4_2
## 1 NA NA
## 2 1 1
## 3 1 1
## 4 1 1
## 5 1 1
## 6 NA 1
## Experiment.3D8_20_2h_n4_3 Experiment.3D8_20_2h_n5_1
## 1 NA NA
## 2 1 1
## 3 1 NA
## 4 1 1
## 5 NA 1
## 6 NA NA
## Experiment.3D8_20_2h_n5_2 Experiment.3D8_20_2h_n5_3 Intensity Reverse
## 1 NA NA 1.3295e+09
## 2 NA NA 1.6719e+08
## 3 1 1 5.5675e+08
## 4 1 1 1.4830e+10
## 5 1 1 4.8209e+09
## 6 NA NA 1.3715e+08
## Contaminant id Protein.group.IDs Mod..peptide.IDs
## 1 0 419 0
## 2 1 1150 1
## 3 2 513 2
## 4 3 199 3
## 5 4 388 4
## 6 5 624 5
## Evidence.IDs
## 1 0;1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18
## 2 19;20;21;22;23;24;25;26;27;28;29;30;31;32;33
## 3 34;35;36;37;38;39;40;41;42;43;44;45;46;47;48;49;50;51;52;53;54;55;56;57;58;59;60;61;62;63;64;65;66;67;68;69;70;71;72;73;74;75;76
## 4 77;78;79;80;81;82;83;84;85;86;87;88;89;90;91;92;93;94;95;96;97;98;99;100;101;102;103;104;105;106;107;108;109;110;111;112;113;114;115;116;117;118;119;120;121;122;123;124
## 5 125;126;127;128;129;130;131;132;133;134;135;136;137;138;139;140;141;142;143;144;145;146;147;148;149;150;151;152;153;154;155
## 6 156;157;158;159;160
## MS.MS.IDs
## 1 0;1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18
## 2 19;20;21;22;23;24;25;26;27;28;29;30;31;32;33
## 3 34;35;36;37;38;39;40;41;42;43;44;45;46;47;48;49;50;51;52;53;54;55;56;57;58;59;60;61;62;63;64;65;66;67;68;69;70;71;72;73;74;75;76;77;78;79
## 4 80;81;82;83;84;85;86;87;88;89;90;91;92;93;94;95;96;97;98;99;100;101;102;103;104;105;106;107;108;109;110;111;112;113;114;115;116;117;118;119;120;121;122;123;124;125;126;127;128;129;130;131;132;133;134;135;136;137;138;139;140;141;142;143;144;145;146;147;148;149;150;151;152;153;154;155;156;157;158;159;160;161;162;163;164;165;166;167;168;169;170;171;172
## 5 173;174;175;176;177;178;179;180;181;182;183;184;185;186;187;188;189;190;191;192;193;194;195;196;197;198;199;200;201;202;203
## 6 204;205;206;207;208
## Best.MS.MS Oxidation..M..site.IDs
## 1 15
## 2 26
## 3 39
## 4 159
## 5 191
## 6 205The pData slot will contain the experimental annotation. Right now, it has 18 rows (equal to the number of mass spec runs) and 0 columns (as the experimental annotation is not yet loaded).
pData(peptidesFranc)## data frame with 0 columns and 18 rows1.2. Importing other data types
The import2MSnSet has built-in support for moFF (Argentini et al., 2016), mzTab (Griss et al., 2014) and Progenesis (Nonlinear Dynamics, Newcastle upon Tyne, U.K.) output. You should then simply set filetype to "moFF", "mzTab" or "Progenesis" respectively.
If you are using yet another data type, you can use MSqRob’s read2MSnSet function to import your data. Hereby, you have to specify (1) your file, (2) a regular expression pattern that indicates the columns with the peptide intensities or simply the indices of these columns and (3) the separator (e.g. "\t", ",", …).
2. Preprocessing
2.1. Create an experimental annotation data frame.
Now, we are going to create a data frame that contains the experimental annotation. This data frame should contain all explanatory variables that groups different runs together. This is important in order to be able to include factors on which you will do inference (such as treatment effects or effects over time, or, in our case: the effect of the genetic knock-out) as well as other factors that might have an influence (such as biological repeats, technical repeats and potential confounders such as batch effects).
Note that the “experimental annotation” can also be given as an Excel file or a tab delimited file! This file should have the column names as a header (i.e. the first row in the file) and the same structure as the exp_annotation data frame we create below. If you have your experimental annotation in a file, just provide the filepath to this annotation file in the preprocess_MaxQuant function.
Also note that exactly one column in the experimental annotation should contain the names of the mass spec runs. These names should be equal to the names you entered in the “Experiment” column in MaxQuant (i.e. the column names of the exprs slot of the peptidesFranc MSnSet object). In our case, we will call this column “run” (see below). By default, MSqRob guesses this column by taking the only column of which its elements correspond to the column names of the exprs slot. If there is no such a column, or there are multiple such columns, MSqRob throws an error.
In this example, we create the experimental annotation data frame ourselves.
Extract the names of the mass spec runs. The mass spec run names are equal to the column names of the exprs slot of the peptidesFranc MSnSet object. Other columns in the experimental annotation data frame will typically be derived from this column.
runs <- sampleNames(peptidesFranc)
runs## [1] "X1WT_20_2h_n3_1" "X1WT_20_2h_n3_2" "X1WT_20_2h_n3_3"
## [4] "X1WT_20_2h_n4_1" "X1WT_20_2h_n4_2" "X1WT_20_2h_n4_3"
## [7] "X1WT_20_2h_n5_1" "X1WT_20_2h_n5_2" "X1WT_20_2h_n5_3"
## [10] "X3D8_20_2h_n3_1" "X3D8_20_2h_n3_2" "X3D8_20_2h_n3_3"
## [13] "X3D8_20_2h_n4_1" "X3D8_20_2h_n4_2" "X3D8_20_2h_n4_3"
## [16] "X3D8_20_2h_n5_1" "X3D8_20_2h_n5_2" "X3D8_20_2h_n5_3"colnames(exprs(peptidesFranc))## [1] "X1WT_20_2h_n3_1" "X1WT_20_2h_n3_2" "X1WT_20_2h_n3_3"
## [4] "X1WT_20_2h_n4_1" "X1WT_20_2h_n4_2" "X1WT_20_2h_n4_3"
## [7] "X1WT_20_2h_n5_1" "X1WT_20_2h_n5_2" "X1WT_20_2h_n5_3"
## [10] "X3D8_20_2h_n3_1" "X3D8_20_2h_n3_2" "X3D8_20_2h_n3_3"
## [13] "X3D8_20_2h_n4_1" "X3D8_20_2h_n4_2" "X3D8_20_2h_n4_3"
## [16] "X3D8_20_2h_n5_1" "X3D8_20_2h_n5_2" "X3D8_20_2h_n5_3"Note again that the sample names have a prepended “X” because a syntactically valid name in R cannot start with a number.
Add a factor for the genetic knock-out, which we will call “genotype” (i.e. wild type (WT) or knock-out (KO)).
#Create the "genotype" vector
genotype <- factor(c(rep("WT",9),rep("KO",9)), levels=c("WT","KO"))
genotype## [1] WT WT WT WT WT WT WT WT WT KO KO KO KO KO KO KO KO KO
## Levels: WT KOAdd a factor for each of the 6 biological replicates (i.e. 3 biological repeats for each genotype).
#Create a factor for the biological replicates
n <- nchar(runs)
biorep <- as.factor(paste0("b_",rep(1:6,each=3)))
biorep## [1] b_1 b_1 b_1 b_2 b_2 b_2 b_3 b_3 b_3 b_4 b_4 b_4 b_5 b_5 b_5 b_6 b_6
## [18] b_6
## Levels: b_1 b_2 b_3 b_4 b_5 b_6#Alternative: extract the biological repeat based on the mass spec run names.
# biorep <- as.factor(paste0("b_",factor(as.numeric(factor(paste0(substr(runs, 1,1),substr(runs, n-2,n-2))))+c(rep(0,9),rep(3,9)))))Each of the 18 technical repeats are here represented by the 18 mass spec runs.
Create an experimental annotation data frame.
exp_annotation <- data.frame(run=runs, genotype=genotype, biorep=biorep)
exp_annotation## run genotype biorep
## 1 X1WT_20_2h_n3_1 WT b_1
## 2 X1WT_20_2h_n3_2 WT b_1
## 3 X1WT_20_2h_n3_3 WT b_1
## 4 X1WT_20_2h_n4_1 WT b_2
## 5 X1WT_20_2h_n4_2 WT b_2
## 6 X1WT_20_2h_n4_3 WT b_2
## 7 X1WT_20_2h_n5_1 WT b_3
## 8 X1WT_20_2h_n5_2 WT b_3
## 9 X1WT_20_2h_n5_3 WT b_3
## 10 X3D8_20_2h_n3_1 KO b_4
## 11 X3D8_20_2h_n3_2 KO b_4
## 12 X3D8_20_2h_n3_3 KO b_4
## 13 X3D8_20_2h_n4_1 KO b_5
## 14 X3D8_20_2h_n4_2 KO b_5
## 15 X3D8_20_2h_n4_3 KO b_5
## 16 X3D8_20_2h_n5_1 KO b_6
## 17 X3D8_20_2h_n5_2 KO b_6
## 18 X3D8_20_2h_n5_3 KO b_6Note that you can also import the experimental annotation from a tab-delimited or Excel file. As an example, we show you here how you can import data from an experimental annotation file delivered with our package. First, provide the path to where the data is saved.
file_exp_annotation <- system.file("extdata/Francisella", "label-free_Francisella_annotation.xlsx", package = "MSqRobData")
file_exp_annotation## [1] "/Library/Frameworks/R.framework/Versions/3.5/Resources/library/MSqRobData/extdata/Francisella/label-free_Francisella_annotation.xlsx"Then, import the data using the read.xlsx function from the openxlsx package and convert the columns of the data frame to factors.
#Import Excel file
exp_annotation_xlsx <- openxlsx::read.xlsx(file_exp_annotation)
#Convert columns of data frame to factors
exp_annotation_xlsx <- as.data.frame(unclass(exp_annotation_xlsx))
exp_annotation_xlsx## run genotype biorep
## 1 1WT_20_2h_n3_1 WT b_1
## 2 1WT_20_2h_n3_2 WT b_1
## 3 1WT_20_2h_n3_3 WT b_1
## 4 1WT_20_2h_n4_1 WT b_2
## 5 1WT_20_2h_n4_2 WT b_2
## 6 1WT_20_2h_n4_3 WT b_2
## 7 1WT_20_2h_n5_1 WT b_3
## 8 1WT_20_2h_n5_2 WT b_3
## 9 1WT_20_2h_n5_3 WT b_3
## 10 3D8_20_2h_n3_1 KO b_4
## 11 3D8_20_2h_n3_2 KO b_4
## 12 3D8_20_2h_n3_3 KO b_4
## 13 3D8_20_2h_n4_1 KO b_5
## 14 3D8_20_2h_n4_2 KO b_5
## 15 3D8_20_2h_n4_3 KO b_5
## 16 3D8_20_2h_n5_1 KO b_6
## 17 3D8_20_2h_n5_2 KO b_6
## 18 3D8_20_2h_n5_3 KO b_6Note that here, there is no “X” prepended to the run names. However, names will be changed automatically to syntactically valid names by MSqRob during preprocessing.
2.2. Preprocess the data using the preprocess_MaxQuant function.
The identification of peptides that carry one or more chemical modifications is often not as reliable as the identification of unmodified peptides. Therefore, a common step during the preprocessing of proteomics data consists of removing the proteins of which all identified peptides have one or more modification sites. If we want to remove these proteins that are only identified by modified peptides from the dataset, we also need the proteinGroups.txt file. This file can be found in the same location as the peptides.txt file (“path_to_raw_files/combined/txt/proteinGroups.txt”). Here, we make use of the proteinGroups.txt file that is delivered with the package.
Note: if you don’t want to remove proteins that are only identified by modified peptides, set the remove_only_site argument of the preprocess_MaxQuant function to FALSE and leave the file_proteinGroups argument at its default value NULL.
Define the path to the proteinGroups.txt file.
file_proteinGroups <- system.file("extdata/Francisella", "proteinGroups.txt", package = "MSqRobData")
file_proteinGroups## [1] "/Library/Frameworks/R.framework/Versions/3.5/Resources/library/MSqRobData/extdata/Francisella/proteinGroups.txt"Add to the useful_properties argument all column names of the fData slot of the peptidesFranc MSnSet object that you would like to keep (such as gene names, protein names, ontologies, etc.). Other columns will be removed for improved efficiency and memory usage. Here, we would like to keep the GI number and the protein names. These 2 columns are present in our peptides.txt file, but note that you can use basic R manipulations to add other columns to the fData slot if you want to. We certainly want to keep the “Sequence” slot, as peptide-specific effects must be included in our model.
peptidesFranc2 <- preprocess_MaxQuant(peptidesFranc, accession="Proteins", exp_annotation=exp_annotation, logtransform=TRUE, base=2, normalisation="quantiles", smallestUniqueGroups=TRUE, useful_properties=c("GI number","Protein.names","Sequence"), filter=c("Contaminant","Reverse"), remove_only_site=TRUE, file_proteinGroups=file_proteinGroups, filter_symbol="+", minIdentified=2)preprocess_MaxQuant internally executes the following preprocessing steps in this order:
Peptide intensities are (log2-)transformed (
logtransform=TRUE, base=2).As a normalization approach, we choose to quantile normalize the peptide intensities (
normalisation="quantiles"). Other options are"sum","max","center.mean","center.median","quantiles.robust","vsn"or"none"(seeMSnbase::normalizefor more information).Handling overlapping protein groups (
smallestUniqueGroups=TRUE): in our approach a peptide can map to multiple accessions (accession="Proteins"), as long as there is none of these proteins present in a smaller subgroup. Peptides that map to protein groups for wich also subgroups are present in the dataset are removed from the dataset. (Proteins in a protein group are separated by a;).Remove contaminants and reverse sequences (
filter=c("Contaminant","Reverse"), filter_symbol="+"): each row with the filter symbol (“+” in our case) in the columns “Contaminant” and “Reverse” will be removed from the data. Note that the “Contaminant” column in older MaxQuant peptides.txt files can also be called “Potential.contaminant”.Remove all proteins that are only identified by peptides carrying a modification site (
remove_only_site=TRUE,file_proteinGroups=file_proteinGroups).Remove all superfluous columns in the
fDataslot, except the “Proteins”, “GI.number”, “Protein.names” and “Sequence” columns (useful_properties=c("GI.number","Protein.names","Sequence")). Note that the “Proteins” column does not have to be provided: theaccessioncolumn is always retained.Remove all peptides that are only identified once in the dataset (
minIdentified=2).Add experimental annotation (
exp_annotation=exp_annotation).
We can now make density plots to inspect the effect of the preprocessing on the distribution of our log2-transformed peptide intensities.
3. Make diagnostic plots
Here we take a closer look at how you can make the diagnostic plots from the Shiny app. We will plot both the density plots and the MDS plot.
#Make sure the order of the rows in your exp_annotation is the same as the columns of the exprs slot in your peptides object.
all(exp_annotation[,"run"]==colnames(exprs(peptidesFranc)))## [1] TRUE#Color by biological repeat
colors <- as.numeric(exp_annotation[,"biorep"])
### Density plot before preprocessing ###
densAll <- apply(log2(exprs(peptidesFranc)),2,density,na.rm=TRUE)
#Automatically determine the plot window size
ymax=max(vapply(densAll,function(d) max(d$y),1))
rangematrix <- vapply(densAll,function(d) range(d$x, na.rm=TRUE), c(1,1)) #no longer
xlim=range(rangematrix,na.rm=TRUE)
ylim=c(0,ymax)
plot(densAll[[1]],col=colors[1],xlim=xlim,ylim=ylim, las=1, frame.plot=FALSE, main="Density before preprocessing")
for (i in 2:length(densAll)) lines(densAll[[i]],col=colors[i])
######
### Density plot after preprocessing ###
densAll <- apply(log2(exprs(peptidesFranc2)),2,density,na.rm=TRUE)
#Automatically determine the plot window size
ymax=max(vapply(densAll,function(d) max(d$y),1))
rangematrix <- vapply(densAll,function(d) range(d$x, na.rm=TRUE), c(1,1)) #no longer
xlim=range(rangematrix,na.rm=TRUE)
ylim=c(0,ymax)
plot(densAll[[1]],col=colors[1],xlim=xlim,ylim=ylim, las=1, frame.plot=FALSE, main="Density after preprocessing")
for (i in 2:length(densAll)) lines(densAll[[i]],col=colors[i])
######We can also construct an MDS plot. This plot will cluster similar runs together.
plotMDS(exprs(peptidesFranc2), col=colors, las=1)
4. Convert the MSnSet object peptidesFranc2 into a protdata object proteinsFranc.
The accession argument denotes by which column the individual peptide observations should be grouped. The annotations argument indicates which columns in peptidesFranc2 are annotations linked to the accessions (in our case: the GI number and the protein names). These annotations will be added to a separate annotation slot.
proteinsFranc <- MSnSet2protdata(peptidesFranc2, accession="Proteins", annotations = c("GI.number", "Protein.names"))Advanced: add information on the design after creating your protdata object
Imagine you forgot a factor in your experimental annotation. Should you redo the whole preprocessing part? No, often, you can still add this factor thanks to the addVarFromVar function. We will now show how we can add an extra (redundant) factor called techrep.
Create a factor for each of the 18 technical replicates (each sample was analyzed in triplicate). This factor is completely redundant, as the “run” column in the experimental annotation already has 18 different levels and is thus already equal to the factor for the technical replicates.
techrep <- as.factor(paste0("t_",1:18))
names(techrep) <- levels(exp_annotation$run)
techrep## X1WT_20_2h_n3_1 X1WT_20_2h_n3_2 X1WT_20_2h_n3_3 X1WT_20_2h_n4_1
## t_1 t_2 t_3 t_4
## X1WT_20_2h_n4_2 X1WT_20_2h_n4_3 X1WT_20_2h_n5_1 X1WT_20_2h_n5_2
## t_5 t_6 t_7 t_8
## X1WT_20_2h_n5_3 X3D8_20_2h_n3_1 X3D8_20_2h_n3_2 X3D8_20_2h_n3_3
## t_9 t_10 t_11 t_12
## X3D8_20_2h_n4_1 X3D8_20_2h_n4_2 X3D8_20_2h_n4_3 X3D8_20_2h_n5_1
## t_13 t_14 t_15 t_16
## X3D8_20_2h_n5_2 X3D8_20_2h_n5_3
## t_17 t_18
## 18 Levels: t_1 t_10 t_11 t_12 t_13 t_14 t_15 t_16 t_17 t_18 t_2 ... t_9proteinsFranc2 <- addVarFromVar(proteinsFranc, basecol="run", name="techrep", vector=techrep)
proteinsFranc2## @accession
##
## [1] "WP_003038655"
##
## @annotation
##
## GI.number
## "gi|118497196"
## Protein.names
## "hypothetical protein [Francisella tularensis]"
##
## @data
##
## quant_value Sequence run genotype biorep techrep
## 1.1 25.11978 AVVVFFDYGCGK X1WT_20_2h_n3_1 WT b_1 t_1
## 1.2 25.34586 DLLSDPATPTVGPQDANK X1WT_20_2h_n3_1 WT b_1 t_1
## 1.3 24.11491 ETNGELTVQDVDNVVK X1WT_20_2h_n3_1 WT b_1 t_1
## 1.4 26.50674 IIAIPDIVK X1WT_20_2h_n3_1 WT b_1 t_1
## 1.5 24.69046 QTNVITDNYQQAAK X1WT_20_2h_n3_1 WT b_1 t_1
## 2.1 24.75196 AVVVFFDYGCGK X1WT_20_2h_n3_2 WT b_1 t_2
##
## 96 more rows...
##
## @accession
##
## [1] "WP_003041237"
##
## @annotation
##
## GI.number
## "gi|118498194"
## Protein.names
## "phosphoglucosamine mutase [Francisella tularensis]"
##
## @data
##
## quant_value Sequence run genotype biorep techrep
## 1.1 25.81732 ADLGISLDGDADR X1WT_20_2h_n3_1 WT b_1 t_1
## 1.2 27.36716 FVIVGQDTR X1WT_20_2h_n3_1 WT b_1 t_1
## 1.3 24.91571 GEVANSTITAEFTQK X1WT_20_2h_n3_1 WT b_1 t_1
## 1.4 27.38539 ILANAIDEYIESIHSR X1WT_20_2h_n3_1 WT b_1 t_1
## 1.6 27.93398 SLATNEAEYLVEK X1WT_20_2h_n3_1 WT b_1 t_1
## 1.7 26.59163 VASDVNDVEK X1WT_20_2h_n3_1 WT b_1 t_1
##
## 187 more rows...
##
## @accession
##
## [1] "WP_003038915"
##
## @annotation
##
## GI.number
## "gi|118497331"
## Protein.names
## "glycine--tRNA ligase subunit beta [Francisella tularensis]"
##
## @data
##
## quant_value Sequence run genotype biorep
## 1 24.37301 AAANEYELALAYSIEEVAPDLHK X1WT_20_2h_n3_1 WT b_1
## 1.3 22.60222 EEGIAVDIFEAINNTNYDSIK X1WT_20_2h_n3_1 WT b_1
## 1.4 25.34002 EGEPTQVGLGFAK X1WT_20_2h_n3_1 WT b_1
## 1.5 25.13163 EYSYALELLTCLDK X1WT_20_2h_n3_1 WT b_1
## 1.6 24.09977 FHHPQAIVIDNISDYIK X1WT_20_2h_n3_1 WT b_1
## 1.8 25.75712 IIDITLESYK X1WT_20_2h_n3_1 WT b_1
## techrep
## 1 t_1
## 1.3 t_1
## 1.4 t_1
## 1.5 t_1
## 1.6 t_1
## 1.8 t_1
##
## 215 more rows...
##
## @accession
##
## [1] "WP_003038264"
##
## @annotation
##
## GI.number
## "gi|118496879"
## Protein.names
## "molecular chaperone HtpG [Francisella tularensis]"
##
## @data
##
## quant_value Sequence run genotype biorep techrep
## 1 27.92076 AAANNPQLEAFK X1WT_20_2h_n3_1 WT b_1 t_1
## 1.2 26.48741 EGQDTIYYITSDSYK X1WT_20_2h_n3_1 WT b_1 t_1
## 1.4 27.33299 EHLDLLEYHVLK X1WT_20_2h_n3_1 WT b_1 t_1
## 1.5 26.37656 EILQHNK X1WT_20_2h_n3_1 WT b_1 t_1
## 1.6 28.11488 ELLPPYLR X1WT_20_2h_n3_1 WT b_1 t_1
## 1.7 28.77850 ELVSNSSDAIEK X1WT_20_2h_n3_1 WT b_1 t_1
##
## 338 more rows...
##
## @accession
##
## [1] "WP_003035781"
##
## @annotation
##
## GI.number
## "gi|118497152"
## Protein.names
## "delta-aminolevulinic acid dehydratase [Francisella tularensis]"
##
## @data
##
## quant_value Sequence run genotype biorep techrep
## 1 27.47327 AAASAGLVDEK X1WT_20_2h_n3_1 WT b_1 t_1
## 1.1 24.05872 VILTYTALDVANWLK X1WT_20_2h_n3_1 WT b_1 t_1
## 2.1 23.88890 VILTYTALDVANWLK X1WT_20_2h_n3_2 WT b_1 t_2
## 2.2 25.37749 YASAYYGPFR X1WT_20_2h_n3_2 WT b_1 t_2
## 3 27.47269 AAASAGLVDEK X1WT_20_2h_n3_3 WT b_1 t_3
## 3.1 24.93174 VILTYTALDVANWLK X1WT_20_2h_n3_3 WT b_1 t_3
##
## 41 more rows...
##
## @accession
##
## [1] "WP_003039212"
##
## @annotation
##
## GI.number
## "gi|118497492"
## Protein.names
## "phosphorylase [Francisella tularensis]"
##
## @data
##
## quant_value Sequence run genotype biorep
## 1.2 22.96510 FPAGLNDVEFVAEHIFQHSK X1WT_20_2h_n3_1 WT b_1
## 1.3 25.10257 LQDLLTFIGK X1WT_20_2h_n3_1 WT b_1
## 1.6 24.55598 VVFSVDYR X1WT_20_2h_n3_1 WT b_1
## 2 24.53369 AAASLDLYSYPK X1WT_20_2h_n3_2 WT b_1
## 2.1 24.68029 EVLFLNR X1WT_20_2h_n3_2 WT b_1
## 2.2 23.16038 FPAGLNDVEFVAEHIFQHSK X1WT_20_2h_n3_2 WT b_1
## techrep
## 1.2 t_1
## 1.3 t_1
## 1.6 t_1
## 2 t_2
## 2.1 t_2
## 2.2 t_2
##
## 56 more rows...
##
## 1182 more elements...5. Inspect the data for a specific protein via a detailed protein plot
Here, we illustrate how you can reproduce the detail plots from the Shiny app. We will take a closer look at the protein “WP_003039212”.
protein_name <- "WP_003039212"
data <- getData(proteinsFranc[protein_name])
#We plot the peptide sequences as different point symbols
points <- as.numeric(droplevels(as.factor(data$Sequence)))
#We don't want to use repeated plot symbols or invisible symbols
pch_vals <- c(0:25,32:127)
#Repeat pch_vals if there would be more than 122 unique levels
pch_vals <- rep_len(pch_vals, length(unique(points)))
pch_vals <- pch_vals[points]
#We color by mass spec run
colors <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(8,"Spectral"))(length(unique(droplevels(as.factor(data$run)))))
colors <- colors[as.numeric(droplevels(as.factor(data$run)))]
#We make a different boxplot for each genotype
boxplot(data$quant_value~as.factor(data$genotype), outline=FALSE, ylim=c(min(data$quant_value)-0.2,max(data$quant_value)+0.2), ylab="preprocessed peptide intensity", xlab="", main=paste0("Detail plot for protein ", protein_name), las=2, frame.plot=FALSE, frame=FALSE, col="grey", pars=list(boxcol="white"))
#Add the data as separate points
points(jitter((as.numeric(data$genotype)), factor=2),data$quant_value, col=colors, pch=pch_vals)
6. Fit the model
Fit the robust ridge models for each protein. The fit.model function will fit a model to each protein in the protdata object proteinsFranc. The default setting reg="ridge" indicates that ridge regression should be used. The default setting weights="Huber" indicates that an M estimation algorithm should be used that assigns Huber weights to each peptide observation based on their residuals.
Caution: fitting the models takes some time (less than 15 minutes for the Francisella tularensis dataset on our system).
system.time(protLMFranc <- fit.model(proteinsFranc, response="quant_value", fixed=c("genotype"), random=c("biorep","run","Sequence"), add.intercept=TRUE, weights="Huber"))## user system elapsed
## 384.990 1.433 396.099The system.time function was only added to output the duration of the model fitting. If you are not interested in this, you can of course omit it.
5.1. Inspect the model
protLMFranc is a protLM object.
class(protLMFranc)## [1] "protLM"
## attr(,"package")
## [1] "MSqRob"Printing out a protLM object returns an overview of the first 6 accessions with their annotations and their corresponding models.
protLMFranc## @accession
##
## [1] "WP_003038655"
##
## @annotation
##
## GI.number
## "gi|118497196"
## Protein.names
## "hypothetical protein [Francisella tularensis]"
##
## @model
##
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 66.4301
## Random effects:
## Groups Name Std.Dev.
## run (Intercept) 0.0000
## Sequence (Intercept) 0.6682
## biorep (Intercept) 0.2247
## ridgeGroup.1 (Intercept) 1.0983
## Residual 0.2524
## Number of obs: 102, groups:
## run, 18; Sequence, 7; biorep, 6; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -260
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
##
## @accession
##
## [1] "WP_003041237"
##
## @annotation
##
## GI.number
## "gi|118498194"
## Protein.names
## "phosphoglucosamine mutase [Francisella tularensis]"
##
## @model
##
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 136.4625
## Random effects:
## Groups Name Std.Dev.
## run (Intercept) 0.00000
## Sequence (Intercept) 1.28720
## biorep (Intercept) 0.09322
## ridgeGroup.1 (Intercept) 0.00000
## Residual 0.27463
## Number of obs: 193, groups:
## run, 18; Sequence, 12; biorep, 6; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -370.9
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
##
## @accession
##
## [1] "WP_003038915"
##
## @annotation
##
## GI.number
## "gi|118497331"
## Protein.names
## "glycine--tRNA ligase subunit beta [Francisella tularensis]"
##
## @model
##
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 340.3599
## Random effects:
## Groups Name Std.Dev.
## Sequence (Intercept) 1.07373
## run (Intercept) 0.00000
## biorep (Intercept) 0.05502
## ridgeGroup.1 (Intercept) 2.11989
## Residual 0.40125
## Number of obs: 221, groups:
## Sequence, 18; run, 18; biorep, 6; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -365.9
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
##
## @accession
##
## [1] "WP_003038264"
##
## @annotation
##
## GI.number
## "gi|118496879"
## Protein.names
## "molecular chaperone HtpG [Francisella tularensis]"
##
## @model
##
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 378.4025
## Random effects:
## Groups Name Std.Dev.
## Sequence (Intercept) 1.02426
## run (Intercept) 0.09019
## biorep (Intercept) 0.11985
## ridgeGroup.1 (Intercept) 0.38438
## Residual 0.31793
## Number of obs: 344, groups:
## Sequence, 23; run, 18; biorep, 6; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -489.8
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
##
## @accession
##
## [1] "WP_003035781"
##
## @annotation
##
## GI.number
## "gi|118497152"
## Protein.names
## "delta-aminolevulinic acid dehydratase [Francisella tularensis]"
##
## @model
##
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 57.1314
## Random effects:
## Groups Name Std.Dev.
## run (Intercept) 0.0000
## biorep (Intercept) 0.2059
## Sequence (Intercept) 1.9048
## ridgeGroup.1 (Intercept) 0.5694
## Residual 0.3522
## Number of obs: 47, groups:
## run, 18; biorep, 6; Sequence, 3; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -173.3
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
##
## @accession
##
## [1] "WP_003039212"
##
## @annotation
##
## GI.number
## "gi|118497492"
## Protein.names
## "phosphorylase [Francisella tularensis]"
##
## @model
##
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 9.2867
## Random effects:
## Groups Name Std.Dev.
## run (Intercept) 0.0000
## Sequence (Intercept) 0.6635
## biorep (Intercept) 0.0000
## ridgeGroup.1 (Intercept) 0.3159
## Residual 0.2026
## Number of obs: 62, groups:
## run, 18; Sequence, 7; biorep, 6; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -192.6
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
##
## 1182 more elements...The numbers in the Std.Dev. column give the estimated standard deviations of the random effects. Effects starting with ridgeGroup indicate shrunken fixed effects (we make use of the random effects framework of the lm4 package in order to estimate Ridge penalties). However, as the fixed effects in the ridgeGroups are orthogonalized using the Gramm-Schmidt procedure, these numbers are not very informative.
A protLM object has three slots:
- an
accessionslot: a vector containing all the accessions - a
modelslot: a list of all fitted models (one for each accession) - an
annotationslot: data frame of which each row (one for each accession) corresponds to the annotations of a different accession.
These slots can be accessed using the getAccessions, getModels and getAnnotations functions respectively. Here we give the accessions, models and annotations for the first 5 models.
getAccessions(protLMFranc)[1:5]## [1] WP_003038655 WP_003041237 WP_003038915 WP_003038264 WP_003035781
## 1188 Levels: WP_003013731 WP_003013860 WP_003013909 ... WP_011733734getModels(protLMFranc)[1:5]## [[1]]
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 66.4301
## Random effects:
## Groups Name Std.Dev.
## run (Intercept) 0.0000
## Sequence (Intercept) 0.6682
## biorep (Intercept) 0.2247
## ridgeGroup.1 (Intercept) 1.0983
## Residual 0.2524
## Number of obs: 102, groups:
## run, 18; Sequence, 7; biorep, 6; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -260
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
##
## [[2]]
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 136.4625
## Random effects:
## Groups Name Std.Dev.
## run (Intercept) 0.00000
## Sequence (Intercept) 1.28720
## biorep (Intercept) 0.09322
## ridgeGroup.1 (Intercept) 0.00000
## Residual 0.27463
## Number of obs: 193, groups:
## run, 18; Sequence, 12; biorep, 6; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -370.9
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
##
## [[3]]
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 340.3599
## Random effects:
## Groups Name Std.Dev.
## Sequence (Intercept) 1.07373
## run (Intercept) 0.00000
## biorep (Intercept) 0.05502
## ridgeGroup.1 (Intercept) 2.11989
## Residual 0.40125
## Number of obs: 221, groups:
## Sequence, 18; run, 18; biorep, 6; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -365.9
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
##
## [[4]]
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 378.4025
## Random effects:
## Groups Name Std.Dev.
## Sequence (Intercept) 1.02426
## run (Intercept) 0.09019
## biorep (Intercept) 0.11985
## ridgeGroup.1 (Intercept) 0.38438
## Residual 0.31793
## Number of obs: 344, groups:
## Sequence, 23; run, 18; biorep, 6; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -489.8
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
##
## [[5]]
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 57.1314
## Random effects:
## Groups Name Std.Dev.
## run (Intercept) 0.0000
## biorep (Intercept) 0.2059
## Sequence (Intercept) 1.9048
## ridgeGroup.1 (Intercept) 0.5694
## Residual 0.3522
## Number of obs: 47, groups:
## run, 18; biorep, 6; Sequence, 3; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -173.3
## fit warnings:
## Some predictor variables are on very different scales: consider rescalinggetAnnotations(protLMFranc)[1:5,]## GI.number
## [1,] "gi|118497196"
## [2,] "gi|118498194"
## [3,] "gi|118497331"
## [4,] "gi|118496879"
## [5,] "gi|118497152"
## Protein.names
## [1,] "hypothetical protein [Francisella tularensis]"
## [2,] "phosphoglucosamine mutase [Francisella tularensis]"
## [3,] "glycine--tRNA ligase subunit beta [Francisella tularensis]"
## [4,] "molecular chaperone HtpG [Francisella tularensis]"
## [5,] "delta-aminolevulinic acid dehydratase [Francisella tularensis]"protLM objects can be subsetted using either their numeric index or their corresponding accession.
Retain only the part corresponding to model 24 and model 50.
protLMFranc[c(24,50)]## @accession
##
## [1] "WP_003035135"
##
## @annotation
##
## GI.number
## "gi|118498204"
## Protein.names
## "superoxide dismutase [Francisella tularensis]"
##
## @model
##
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 130.6537
## Random effects:
## Groups Name Std.Dev.
## run (Intercept) 0.00000
## Sequence (Intercept) 1.48360
## biorep (Intercept) 0.05107
## ridgeGroup.1 (Intercept) 1.40062
## Residual 0.29509
## Number of obs: 148, groups:
## run, 18; Sequence, 10; biorep, 6; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -344.8
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
##
## @accession
##
## [1] "WP_003041024"
##
## @annotation
##
## GI.number
## "gi|118496740"
## Protein.names
## "acetate kinase [Francisella tularensis]"
##
## @model
##
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 26.4739
## Random effects:
## Groups Name Std.Dev.
## run (Intercept) 0.0000
## biorep (Intercept) 0.2982
## Sequence (Intercept) 0.3214
## ridgeGroup.1 (Intercept) 0.0000
## Residual 0.2708
## Number of obs: 36, groups:
## run, 16; biorep, 6; Sequence, 4; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -147.4
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingInspect the protLM object corresponding to accession "WP_003033643".
protLMFranc["WP_003033643"]## @accession
##
## [1] "WP_003033643"
##
## @annotation
##
## GI.number
## "gi|118497379"
## Protein.names
## "isochorismatase [Francisella tularensis]"
##
## @model
##
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 5.8471
## Random effects:
## Groups Name Std.Dev.
## run (Intercept) 0.0000
## biorep (Intercept) 0.1829
## Sequence (Intercept) 1.1305
## ridgeGroup.1 (Intercept) 0.3079
## Residual 0.9620
## Number of obs: 3, groups: run, 3; biorep, 2; Sequence, 1; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -36.64
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingSay we want to inspect cell division protein "WP_003040227". We can then select protein "WP_003040227" from the protLM object. We can extract the model via the getModels function. simplify=TRUE is the default setting and indicates that if you select one model, you only want this model and not a list with as only element this model. Note that if you select multiple models, you automatically receive a list of models.
modelWP_003040227 <- getModels(protLMFranc["WP_003040227"], simplify=TRUE)
modelWP_003040227## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 35.3401
## Random effects:
## Groups Name Std.Dev.
## run (Intercept) 0.0000
## biorep (Intercept) 0.5312
## Sequence (Intercept) 0.0000
## ridgeGroup.1 (Intercept) 0.8636
## Residual 0.4931
## Number of obs: 20, groups:
## run, 15; biorep, 5; Sequence, 2; ridgeGroup.1, 1
## Fixed Effects:
## Intercept_MSqRob
## -111.7
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingThe betaBVcovDf function extracts a list from a model containing the following four elements:
- beta: the estimated parameters: both fixed effect sizes and BLUPs (best linear unbiased predictors) for the random effects
- vcov: the variance-covariance matrix
- df: the residual degrees of freedom
- sigma: the residual standard deviation
- df_exp: if the number of experimental units would be specified: a more conservative estimate for the degrees of freedom based on the number of experimental units minus the degrees of freedom lost due to between-treatment effects. Otherwise:
NULL.
betaBVcovDf <- getBetaVcovDf(modelWP_003040227)betaBVcovDf$beta## [,1]
## (Intercept) 2.478692e+01
## runX1WT_20_2h_n4_1 -3.268340e-19
## runX1WT_20_2h_n4_2 6.312569e-19
## runX1WT_20_2h_n4_3 -5.514836e-19
## runX1WT_20_2h_n5_1 -2.653938e-19
## runX1WT_20_2h_n5_2 1.300820e-19
## runX1WT_20_2h_n5_3 3.002867e-20
## runX3D8_20_2h_n3_1 6.312569e-19
## runX3D8_20_2h_n3_2 5.345092e-20
## runX3D8_20_2h_n3_3 -7.547956e-20
## runX3D8_20_2h_n4_1 3.677317e-19
## runX3D8_20_2h_n4_2 1.418812e-19
## runX3D8_20_2h_n4_3 -3.821573e-19
## runX3D8_20_2h_n5_1 -3.498709e-19
## runX3D8_20_2h_n5_2 -1.035900e-19
## runX3D8_20_2h_n5_3 6.912076e-20
## biorepb_2 -2.866945e-01
## biorepb_3 -1.221728e-01
## biorepb_4 7.069615e-01
## biorepb_5 1.479023e-01
## biorepb_6 -4.459965e-01
## SequenceMVFTFYNNLK 6.400946e-19
## SequenceYTYIYQIASFR -6.400946e-19
## genotypeKO 2.573084e-01betaBVcovDf$vcov## (Intercept) runX1WT_20_2h_n4_1 runX1WT_20_2h_n4_2
## (Intercept) 4.644336e-01 -1.133292e-19 -7.590798e-20
## runX1WT_20_2h_n4_1 -1.133292e-19 1.000000e-18 2.081844e-37
## runX1WT_20_2h_n4_2 -7.590798e-20 2.081844e-37 1.000000e-18
## runX1WT_20_2h_n4_3 -1.133292e-19 3.108154e-37 2.081844e-37
## runX1WT_20_2h_n5_1 -1.036390e-19 2.528956e-38 1.693897e-38
## runX1WT_20_2h_n5_2 -1.036390e-19 2.528956e-38 1.693897e-38
## runX1WT_20_2h_n5_3 -1.036390e-19 2.528956e-38 1.693897e-38
## runX3D8_20_2h_n3_1 -1.868609e-20 4.559704e-39 3.054093e-39
## runX3D8_20_2h_n3_2 -4.832685e-20 1.179252e-38 7.898639e-39
## runX3D8_20_2h_n3_3 -4.832685e-20 1.179252e-38 7.898639e-39
## runX3D8_20_2h_n4_1 -2.678135e-20 6.535074e-39 4.377198e-39
## runX3D8_20_2h_n4_2 -5.356270e-20 1.307015e-38 8.754397e-39
## runX3D8_20_2h_n4_3 -5.356270e-20 1.307015e-38 8.754397e-39
## runX3D8_20_2h_n5_1 -4.575672e-20 1.116537e-38 7.478572e-39
## runX3D8_20_2h_n5_2 -4.575672e-20 1.116537e-38 7.478572e-39
## runX3D8_20_2h_n5_3 -4.575672e-20 1.116537e-38 7.478572e-39
## biorepb_2 -3.511045e-01 -1.974862e-19 -1.322764e-19
## biorepb_3 -3.607947e-01 8.803963e-20 5.896902e-20
## biorepb_4 -1.338428e-01 3.265976e-20 2.187553e-20
## biorepb_5 -1.553883e-01 3.791720e-20 2.539697e-20
## biorepb_6 -1.592912e-01 3.886959e-20 2.603489e-20
## SequenceMVFTFYNNLK -5.668853e-19 -1.056869e-37 -7.078916e-38
## SequenceYTYIYQIASFR -4.331147e-19 1.056869e-37 7.078916e-38
## genotypeKO -2.822640e-01 6.887691e-20 4.613381e-20
## runX1WT_20_2h_n4_3 runX1WT_20_2h_n5_1
## (Intercept) -1.133292e-19 -1.036390e-19
## runX1WT_20_2h_n4_1 3.108154e-37 2.528956e-38
## runX1WT_20_2h_n4_2 2.081844e-37 1.693897e-38
## runX1WT_20_2h_n4_3 1.000000e-18 2.528956e-38
## runX1WT_20_2h_n5_1 2.528956e-38 1.000000e-18
## runX1WT_20_2h_n5_2 2.528956e-38 2.820767e-37
## runX1WT_20_2h_n5_3 2.528956e-38 2.820767e-37
## runX3D8_20_2h_n3_1 4.559704e-39 4.169826e-39
## runX3D8_20_2h_n3_2 1.179252e-38 1.078420e-38
## runX3D8_20_2h_n3_3 1.179252e-38 1.078420e-38
## runX3D8_20_2h_n4_1 6.535074e-39 5.976293e-39
## runX3D8_20_2h_n4_2 1.307015e-38 1.195259e-38
## runX3D8_20_2h_n4_3 1.307015e-38 1.195259e-38
## runX3D8_20_2h_n5_1 1.116537e-38 1.021067e-38
## runX3D8_20_2h_n5_2 1.116537e-38 1.021067e-38
## runX3D8_20_2h_n5_3 1.116537e-38 1.021067e-38
## biorepb_2 -1.974862e-19 7.834942e-20
## biorepb_3 8.803963e-20 -1.784378e-19
## biorepb_4 3.265976e-20 2.986719e-20
## biorepb_5 3.791720e-20 3.467509e-20
## biorepb_6 3.886959e-20 3.554605e-20
## SequenceMVFTFYNNLK -1.056869e-37 1.265012e-37
## SequenceYTYIYQIASFR 1.056869e-37 -1.265012e-37
## genotypeKO 6.887691e-20 6.298759e-20
## runX1WT_20_2h_n5_2 runX1WT_20_2h_n5_3
## (Intercept) -1.036390e-19 -1.036390e-19
## runX1WT_20_2h_n4_1 2.528956e-38 2.528956e-38
## runX1WT_20_2h_n4_2 1.693897e-38 1.693897e-38
## runX1WT_20_2h_n4_3 2.528956e-38 2.528956e-38
## runX1WT_20_2h_n5_1 2.820767e-37 2.820767e-37
## runX1WT_20_2h_n5_2 1.000000e-18 2.820767e-37
## runX1WT_20_2h_n5_3 2.820767e-37 1.000000e-18
## runX3D8_20_2h_n3_1 4.169826e-39 4.169826e-39
## runX3D8_20_2h_n3_2 1.078420e-38 1.078420e-38
## runX3D8_20_2h_n3_3 1.078420e-38 1.078420e-38
## runX3D8_20_2h_n4_1 5.976293e-39 5.976293e-39
## runX3D8_20_2h_n4_2 1.195259e-38 1.195259e-38
## runX3D8_20_2h_n4_3 1.195259e-38 1.195259e-38
## runX3D8_20_2h_n5_1 1.021067e-38 1.021067e-38
## runX3D8_20_2h_n5_2 1.021067e-38 1.021067e-38
## runX3D8_20_2h_n5_3 1.021067e-38 1.021067e-38
## biorepb_2 7.834942e-20 7.834942e-20
## biorepb_3 -1.784378e-19 -1.784378e-19
## biorepb_4 2.986719e-20 2.986719e-20
## biorepb_5 3.467509e-20 3.467509e-20
## biorepb_6 3.554605e-20 3.554605e-20
## SequenceMVFTFYNNLK 1.265012e-37 1.265012e-37
## SequenceYTYIYQIASFR -1.265012e-37 -1.265012e-37
## genotypeKO 6.298759e-20 6.298759e-20
## runX3D8_20_2h_n3_1 runX3D8_20_2h_n3_2
## (Intercept) -1.868609e-20 -4.832685e-20
## runX1WT_20_2h_n4_1 4.559704e-39 1.179252e-38
## runX1WT_20_2h_n4_2 3.054093e-39 7.898639e-39
## runX1WT_20_2h_n4_3 4.559704e-39 1.179252e-38
## runX1WT_20_2h_n5_1 4.169826e-39 1.078420e-38
## runX1WT_20_2h_n5_2 4.169826e-39 1.078420e-38
## runX1WT_20_2h_n5_3 4.169826e-39 1.078420e-38
## runX3D8_20_2h_n3_1 1.000000e-18 1.287695e-37
## runX3D8_20_2h_n3_2 1.287695e-37 1.000000e-18
## runX3D8_20_2h_n3_3 1.287695e-37 3.330298e-37
## runX3D8_20_2h_n4_1 5.397077e-39 1.395817e-38
## runX3D8_20_2h_n4_2 1.079415e-38 2.791634e-38
## runX3D8_20_2h_n4_3 1.079415e-38 2.791634e-38
## runX3D8_20_2h_n5_1 9.221064e-39 2.384795e-38
## runX3D8_20_2h_n5_2 9.221064e-39 2.384795e-38
## runX3D8_20_2h_n5_3 9.221064e-39 2.384795e-38
## biorepb_2 1.412639e-20 3.653433e-20
## biorepb_3 1.451627e-20 3.754265e-20
## biorepb_4 -9.205803e-20 -2.380848e-19
## biorepb_5 3.131442e-20 8.098680e-20
## biorepb_6 3.210096e-20 8.302100e-20
## SequenceMVFTFYNNLK -3.713523e-38 -9.604086e-38
## SequenceYTYIYQIASFR 3.713523e-38 9.604086e-38
## genotypeKO -1.802540e-20 -4.661813e-20
## runX3D8_20_2h_n3_3 runX3D8_20_2h_n4_1
## (Intercept) -4.832685e-20 -2.678135e-20
## runX1WT_20_2h_n4_1 1.179252e-38 6.535074e-39
## runX1WT_20_2h_n4_2 7.898639e-39 4.377198e-39
## runX1WT_20_2h_n4_3 1.179252e-38 6.535074e-39
## runX1WT_20_2h_n5_1 1.078420e-38 5.976293e-39
## runX1WT_20_2h_n5_2 1.078420e-38 5.976293e-39
## runX1WT_20_2h_n5_3 1.078420e-38 5.976293e-39
## runX3D8_20_2h_n3_1 1.287695e-37 5.397077e-39
## runX3D8_20_2h_n3_2 3.330298e-37 1.395817e-38
## runX3D8_20_2h_n3_3 1.000000e-18 1.395817e-38
## runX3D8_20_2h_n4_1 1.395817e-38 1.000000e-18
## runX3D8_20_2h_n4_2 2.791634e-38 3.566651e-37
## runX3D8_20_2h_n4_3 2.791634e-38 3.566651e-37
## runX3D8_20_2h_n5_1 2.384795e-38 1.321585e-38
## runX3D8_20_2h_n5_2 2.384795e-38 1.321585e-38
## runX3D8_20_2h_n5_3 2.384795e-38 1.321585e-38
## biorepb_2 3.653433e-20 2.024628e-20
## biorepb_3 3.754265e-20 2.080506e-20
## biorepb_4 -2.380848e-19 3.865760e-20
## biorepb_5 8.098680e-20 -1.257168e-19
## biorepb_6 8.302100e-20 4.600785e-20
## SequenceMVFTFYNNLK -9.604086e-38 -3.944178e-37
## SequenceYTYIYQIASFR 9.604086e-38 3.944178e-37
## genotypeKO -4.661813e-20 -2.583443e-20
## runX3D8_20_2h_n4_2 runX3D8_20_2h_n4_3
## (Intercept) -5.356270e-20 -5.356270e-20
## runX1WT_20_2h_n4_1 1.307015e-38 1.307015e-38
## runX1WT_20_2h_n4_2 8.754397e-39 8.754397e-39
## runX1WT_20_2h_n4_3 1.307015e-38 1.307015e-38
## runX1WT_20_2h_n5_1 1.195259e-38 1.195259e-38
## runX1WT_20_2h_n5_2 1.195259e-38 1.195259e-38
## runX1WT_20_2h_n5_3 1.195259e-38 1.195259e-38
## runX3D8_20_2h_n3_1 1.079415e-38 1.079415e-38
## runX3D8_20_2h_n3_2 2.791634e-38 2.791634e-38
## runX3D8_20_2h_n3_3 2.791634e-38 2.791634e-38
## runX3D8_20_2h_n4_1 3.566651e-37 3.566651e-37
## runX3D8_20_2h_n4_2 1.000000e-18 7.133303e-37
## runX3D8_20_2h_n4_3 7.133303e-37 1.000000e-18
## runX3D8_20_2h_n5_1 2.643169e-38 2.643169e-38
## runX3D8_20_2h_n5_2 2.643169e-38 2.643169e-38
## runX3D8_20_2h_n5_3 2.643169e-38 2.643169e-38
## biorepb_2 4.049255e-20 4.049255e-20
## biorepb_3 4.161011e-20 4.161011e-20
## biorepb_4 7.731521e-20 7.731521e-20
## biorepb_5 -2.514336e-19 -2.514336e-19
## biorepb_6 9.201570e-20 9.201570e-20
## SequenceMVFTFYNNLK 2.111644e-37 2.111644e-37
## SequenceYTYIYQIASFR -2.111644e-37 -2.111644e-37
## genotypeKO -5.166885e-20 -5.166885e-20
## runX3D8_20_2h_n5_1 runX3D8_20_2h_n5_2
## (Intercept) -4.575672e-20 -4.575672e-20
## runX1WT_20_2h_n4_1 1.116537e-38 1.116537e-38
## runX1WT_20_2h_n4_2 7.478572e-39 7.478572e-39
## runX1WT_20_2h_n4_3 1.116537e-38 1.116537e-38
## runX1WT_20_2h_n5_1 1.021067e-38 1.021067e-38
## runX1WT_20_2h_n5_2 1.021067e-38 1.021067e-38
## runX1WT_20_2h_n5_3 1.021067e-38 1.021067e-38
## runX3D8_20_2h_n3_1 9.221064e-39 9.221064e-39
## runX3D8_20_2h_n3_2 2.384795e-38 2.384795e-38
## runX3D8_20_2h_n3_3 2.384795e-38 2.384795e-38
## runX3D8_20_2h_n4_1 1.321585e-38 1.321585e-38
## runX3D8_20_2h_n4_2 2.643169e-38 2.643169e-38
## runX3D8_20_2h_n4_3 2.643169e-38 2.643169e-38
## runX3D8_20_2h_n5_1 1.000000e-18 6.055208e-37
## runX3D8_20_2h_n5_2 6.055208e-37 1.000000e-18
## runX3D8_20_2h_n5_3 6.055208e-37 6.055208e-37
## biorepb_2 3.459136e-20 3.459136e-20
## biorepb_3 3.554605e-20 3.554605e-20
## biorepb_4 6.604765e-20 6.604765e-20
## biorepb_5 7.667975e-20 7.667975e-20
## biorepb_6 -2.128648e-19 -2.128648e-19
## SequenceMVFTFYNNLK 3.465504e-38 3.465504e-38
## SequenceYTYIYQIASFR -3.465504e-38 -3.465504e-38
## genotypeKO -4.413888e-20 -4.413888e-20
## runX3D8_20_2h_n5_3 biorepb_2 biorepb_3
## (Intercept) -4.575672e-20 -3.511045e-01 -3.607947e-01
## runX1WT_20_2h_n4_1 1.116537e-38 -1.974862e-19 8.803963e-20
## runX1WT_20_2h_n4_2 7.478572e-39 -1.322764e-19 5.896902e-20
## runX1WT_20_2h_n4_3 1.116537e-38 -1.974862e-19 8.803963e-20
## runX1WT_20_2h_n5_1 1.021067e-38 7.834942e-20 -1.784378e-19
## runX1WT_20_2h_n5_2 1.021067e-38 7.834942e-20 -1.784378e-19
## runX1WT_20_2h_n5_3 1.021067e-38 7.834942e-20 -1.784378e-19
## runX3D8_20_2h_n3_1 9.221064e-39 1.412639e-20 1.451627e-20
## runX3D8_20_2h_n3_2 2.384795e-38 3.653433e-20 3.754265e-20
## runX3D8_20_2h_n3_3 2.384795e-38 3.653433e-20 3.754265e-20
## runX3D8_20_2h_n4_1 1.321585e-38 2.024628e-20 2.080506e-20
## runX3D8_20_2h_n4_2 2.643169e-38 4.049255e-20 4.161011e-20
## runX3D8_20_2h_n4_3 2.643169e-38 4.049255e-20 4.161011e-20
## runX3D8_20_2h_n5_1 6.055208e-37 3.459136e-20 3.554605e-20
## runX3D8_20_2h_n5_2 6.055208e-37 3.459136e-20 3.554605e-20
## runX3D8_20_2h_n5_3 1.000000e-18 3.459136e-20 3.554605e-20
## biorepb_2 3.459136e-20 5.485906e-01 2.727550e-01
## biorepb_3 3.554605e-20 2.727550e-01 5.392324e-01
## biorepb_4 6.604765e-20 1.011830e-01 1.039756e-01
## biorepb_5 7.667975e-20 1.174711e-01 1.207132e-01
## biorepb_6 -2.128648e-19 1.204216e-01 1.237452e-01
## SequenceMVFTFYNNLK 3.465504e-38 -3.274279e-19 4.403841e-19
## SequenceYTYIYQIASFR -3.465504e-38 3.274279e-19 -4.403841e-19
## genotypeKO -4.413888e-20 2.133871e-01 2.192765e-01
## biorepb_4 biorepb_5 biorepb_6
## (Intercept) -1.338428e-01 -1.553883e-01 -1.592912e-01
## runX1WT_20_2h_n4_1 3.265976e-20 3.791720e-20 3.886959e-20
## runX1WT_20_2h_n4_2 2.187553e-20 2.539697e-20 2.603489e-20
## runX1WT_20_2h_n4_3 3.265976e-20 3.791720e-20 3.886959e-20
## runX1WT_20_2h_n5_1 2.986719e-20 3.467509e-20 3.554605e-20
## runX1WT_20_2h_n5_2 2.986719e-20 3.467509e-20 3.554605e-20
## runX1WT_20_2h_n5_3 2.986719e-20 3.467509e-20 3.554605e-20
## runX3D8_20_2h_n3_1 -9.205803e-20 3.131442e-20 3.210096e-20
## runX3D8_20_2h_n3_2 -2.380848e-19 8.098680e-20 8.302100e-20
## runX3D8_20_2h_n3_3 -2.380848e-19 8.098680e-20 8.302100e-20
## runX3D8_20_2h_n4_1 3.865760e-20 -1.257168e-19 4.600785e-20
## runX3D8_20_2h_n4_2 7.731521e-20 -2.514336e-19 9.201570e-20
## runX3D8_20_2h_n4_3 7.731521e-20 -2.514336e-19 9.201570e-20
## runX3D8_20_2h_n5_1 6.604765e-20 7.667975e-20 -2.128648e-19
## runX3D8_20_2h_n5_2 6.604765e-20 7.667975e-20 -2.128648e-19
## runX3D8_20_2h_n5_3 6.604765e-20 7.667975e-20 -2.128648e-19
## biorepb_2 1.011830e-01 1.174711e-01 1.204216e-01
## biorepb_3 1.039756e-01 1.207132e-01 1.237452e-01
## biorepb_4 5.010379e-01 2.242956e-01 2.299293e-01
## biorepb_5 2.242956e-01 4.309991e-01 2.669425e-01
## biorepb_6 2.299293e-01 2.669425e-01 4.193827e-01
## SequenceMVFTFYNNLK -2.659883e-19 3.238866e-20 1.206434e-19
## SequenceYTYIYQIASFR 2.659883e-19 -3.238866e-20 -1.206434e-19
## genotypeKO -1.291104e-01 -1.498941e-01 -1.536591e-01
## SequenceMVFTFYNNLK SequenceYTYIYQIASFR genotypeKO
## (Intercept) -5.668853e-19 -4.331147e-19 -2.822640e-01
## runX1WT_20_2h_n4_1 -1.056869e-37 1.056869e-37 6.887691e-20
## runX1WT_20_2h_n4_2 -7.078916e-38 7.078916e-38 4.613381e-20
## runX1WT_20_2h_n4_3 -1.056869e-37 1.056869e-37 6.887691e-20
## runX1WT_20_2h_n5_1 1.265012e-37 -1.265012e-37 6.298759e-20
## runX1WT_20_2h_n5_2 1.265012e-37 -1.265012e-37 6.298759e-20
## runX1WT_20_2h_n5_3 1.265012e-37 -1.265012e-37 6.298759e-20
## runX3D8_20_2h_n3_1 -3.713523e-38 3.713523e-38 -1.802540e-20
## runX3D8_20_2h_n3_2 -9.604086e-38 9.604086e-38 -4.661813e-20
## runX3D8_20_2h_n3_3 -9.604086e-38 9.604086e-38 -4.661813e-20
## runX3D8_20_2h_n4_1 -3.944178e-37 3.944178e-37 -2.583443e-20
## runX3D8_20_2h_n4_2 2.111644e-37 -2.111644e-37 -5.166885e-20
## runX3D8_20_2h_n4_3 2.111644e-37 -2.111644e-37 -5.166885e-20
## runX3D8_20_2h_n5_1 3.465504e-38 -3.465504e-38 -4.413888e-20
## runX3D8_20_2h_n5_2 3.465504e-38 -3.465504e-38 -4.413888e-20
## runX3D8_20_2h_n5_3 3.465504e-38 -3.465504e-38 -4.413888e-20
## biorepb_2 -3.274279e-19 3.274279e-19 2.133871e-01
## biorepb_3 4.403841e-19 -4.403841e-19 2.192765e-01
## biorepb_4 -2.659883e-19 2.659883e-19 -1.291104e-01
## biorepb_5 3.238866e-20 -3.238866e-20 -1.498941e-01
## biorepb_6 1.206434e-19 -1.206434e-19 -1.536591e-01
## SequenceMVFTFYNNLK 1.000000e-18 3.142651e-36 -7.108561e-20
## SequenceYTYIYQIASFR 3.142651e-36 1.000000e-18 7.108561e-20
## genotypeKO -7.108561e-20 7.108561e-20 4.579926e-01betaBVcovDf$df## [1] 18.69393betaBVcovDf$df_exp## NULLbetaBVcovDf$sigma## [1] 0.4618902The getBetaVcovDfList function can be applied on (a subset of) the protLM object protLMFranc. This returns a list containing one or more lists similar to the betaBVcovDf list.
betaVcovDfList <- getBetaVcovDfList(protLMFranc[1:2])
#Inspect the structure of the list without printing it out compeletely:
str(betaVcovDfList)## List of 2
## $ WP_003038655:List of 6
## ..$ beta : num [1:33, 1] 2.57e+01 -9.42e-20 1.96e-19 -5.17e-19 -1.58e-19 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:33] "(Intercept)" "runX1WT_20_2h_n3_1" "runX1WT_20_2h_n3_2" "runX1WT_20_2h_n3_3" ...
## .. .. ..$ : NULL
## ..$ vcov : num [1:33, 1:33] 1.23 -7.08e-20 -1.04e-19 -8.46e-20 -9.38e-20 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:33] "(Intercept)" "runX1WT_20_2h_n3_1" "runX1WT_20_2h_n3_2" "runX1WT_20_2h_n3_3" ...
## .. .. ..$ : chr [1:33] "(Intercept)" "runX1WT_20_2h_n3_1" "runX1WT_20_2h_n3_2" "runX1WT_20_2h_n3_3" ...
## ..$ df : num 93.4
## ..$ sigma : num 0.252
## ..$ df_exp : NULL
## ..$ df_pars: Named num [1:33] 1.00 3.19e-18 3.70e-18 3.48e-18 3.31e-18 ...
## .. ..- attr(*, "names")= chr [1:33] "(Intercept)" "runX1WT_20_2h_n3_1" "runX1WT_20_2h_n3_2" "runX1WT_20_2h_n3_3" ...
## $ WP_003041237:List of 6
## ..$ beta : num [1:38, 1] 2.67e+01 -6.09e-19 -2.51e-20 3.31e-19 -1.80e-18 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:38] "(Intercept)" "runX1WT_20_2h_n3_1" "runX1WT_20_2h_n3_2" "runX1WT_20_2h_n3_3" ...
## .. .. ..$ : NULL
## ..$ vcov : num [1:38, 1:38] 1.86 -5.44e-20 -6.04e-20 -5.10e-20 -3.68e-20 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:38] "(Intercept)" "runX1WT_20_2h_n3_1" "runX1WT_20_2h_n3_2" "runX1WT_20_2h_n3_3" ...
## .. .. ..$ : chr [1:38] "(Intercept)" "runX1WT_20_2h_n3_1" "runX1WT_20_2h_n3_2" "runX1WT_20_2h_n3_3" ...
## ..$ df : num 180
## ..$ sigma : num 0.274
## ..$ df_exp : NULL
## ..$ df_pars: Named num [1:38] 1.00 7.14e-18 7.28e-18 6.98e-18 5.14e-18 ...
## .. ..- attr(*, "names")= chr [1:38] "(Intercept)" "runX1WT_20_2h_n3_1" "runX1WT_20_2h_n3_2" "runX1WT_20_2h_n3_3" ...6.2 Save your results
You can save your proteinsFranc and protLMFranc objects to a .RDatas object that can be loaded in R or in the Shiny graphical user interface. If you want your .RDatas object to be loaded into the Shiny App, you should also specify the levelOptions, random and plot2DependentVars objects. These objects are settings that are saved automatically when using the MSqRob Shiny App. In this example, I save the .RDatas object to my desktop.
#Change to eval=TRUE to execute this piece of code.
result_files_new <- list()
result_files_new$proteins <- proteinsFranc
result_files_new$models <- protLMFranc
### This part is only needed if you want to load your .RDatas file in the MSqRob Shiny App ###
#Should contain all contrast options for the fixed effects (see also 5. Statistical inference)
result_files_new$levelOptions <- c("genotypeWT","genotypeKO")
#Should contain all random effects
result_files_new$random <- c("Sequence","biorep","run")
#Should contain a list with all fixed and random effects
result_files_new$plot2DependentVars <- as.list(c("genotype","Sequence","biorep","run"))
#############################################################################################
saves_MSqRob(result_files, file="/Users/lgoeminn/Desktop/project_Francisella.RDatas") #Change file path to your desired save location!If you want to load your saved file, first specify the folder where it is saved.
#Change to eval=TRUE to execute this piece of code.
result_path <- "/Users/lgoeminn/Desktop/project_Francisella.RDatas"To inspect the different objects in the file, you can use the inspect_loads_MSqRob function.
#Change to eval=TRUE to execute this piece of code.
inspect_loads_MSqRob(result_path)To load this file and assign new protdata and protLM object to it, execute the following piece of code.
#Change to eval=TRUE to execute this piece of code.
#To load the "project_Francisella.RDatas" file
result_files <- loads_MSqRob(result_path)
#assign a new protdata object to the saved proteins
proteinsFranc_result <- result_files$proteins
#assign a new protLM object to the saved models
protLMFranc_result <- result_files$models7. Statistical inference
We are interested in which proteins differ significantly in abundance between ArgP mutant and wild type Francisella tularensis. Via the "MSqRob_levels" attribute, you can check which factor levels are present in each model. Pick a model for which there are no missing levels for the factors for which you want to test contrasts and inspect the names. Note that these names are always the combination of the parameter name and the factor level.
attr(modelWP_003040227,"MSqRob_levels")## [1] "(Intercept)" "genotypeWT" "genotypeKO"
## [4] "runX1WT_20_2h_n4_1" "runX1WT_20_2h_n4_2" "runX1WT_20_2h_n4_3"
## [7] "runX1WT_20_2h_n5_1" "runX1WT_20_2h_n5_2" "runX1WT_20_2h_n5_3"
## [10] "runX3D8_20_2h_n3_1" "runX3D8_20_2h_n3_2" "runX3D8_20_2h_n3_3"
## [13] "runX3D8_20_2h_n4_1" "runX3D8_20_2h_n4_2" "runX3D8_20_2h_n4_3"
## [16] "runX3D8_20_2h_n5_1" "runX3D8_20_2h_n5_2" "runX3D8_20_2h_n5_3"
## [19] "biorepb_2" "biorepb_3" "biorepb_4"
## [22] "biorepb_5" "biorepb_6" "SequenceMVFTFYNNLK"
## [25] "SequenceYTYIYQIASFR"Here, we notice that the mutant is encoded as "genotypeKO" and the wild type as "genotypeWT". We construct a contrast matrix L for inferring on the fold change of interest. We are interested in the log2 fold change between mutant and wild type. As all intensities are log2-transformed, we should take the difference between "genotypeKO" and "genotypeWT" (because of the following property of logarithms: log2(a/b)=log2(a)-log2(b)).
L <- makeContrast(contrasts=c("genotypeKO-genotypeWT"),
levels=c("genotypeWT","genotypeKO"))
L## Contrasts
## Levels genotypeKO-genotypeWT
## genotypeWT -1
## genotypeKO 1Imagine we would want to test the average log2 protein intensity between mutant and wild type, then L would be constructed as follows:
L2 <- makeContrast(contrasts=c("genotypeKO/2+genotypeWT/2"),
levels=c("genotypeWT","genotypeKO"))
L2## Contrasts
## Levels genotypeKO/2+genotypeWT/2
## genotypeWT 0.5
## genotypeKO 0.5Important notification: do not use the estimated parameters beta to determine the factor levels! You will notice for example that the level "genotypeWT" is missing in the column names of the betas. This is because "genotypeWT" is the reference class, giving the parameter corresponding to "genotypeKO" the interpretation of the effect of "genotypeKO" relative to "genotypeWT". It is important however to include "genotypeWT" in the contrast matrix because for proteins where factor levels are missing, these parameters can get different interpretations. The makeContrast function takes this into account and will return NA when a contrast cannot be properly estimated.
betaBVcovDf$beta## [,1]
## (Intercept) 2.478692e+01
## runX1WT_20_2h_n4_1 -3.268340e-19
## runX1WT_20_2h_n4_2 6.312569e-19
## runX1WT_20_2h_n4_3 -5.514836e-19
## runX1WT_20_2h_n5_1 -2.653938e-19
## runX1WT_20_2h_n5_2 1.300820e-19
## runX1WT_20_2h_n5_3 3.002867e-20
## runX3D8_20_2h_n3_1 6.312569e-19
## runX3D8_20_2h_n3_2 5.345092e-20
## runX3D8_20_2h_n3_3 -7.547956e-20
## runX3D8_20_2h_n4_1 3.677317e-19
## runX3D8_20_2h_n4_2 1.418812e-19
## runX3D8_20_2h_n4_3 -3.821573e-19
## runX3D8_20_2h_n5_1 -3.498709e-19
## runX3D8_20_2h_n5_2 -1.035900e-19
## runX3D8_20_2h_n5_3 6.912076e-20
## biorepb_2 -2.866945e-01
## biorepb_3 -1.221728e-01
## biorepb_4 7.069615e-01
## biorepb_5 1.479023e-01
## biorepb_6 -4.459965e-01
## SequenceMVFTFYNNLK 6.400946e-19
## SequenceYTYIYQIASFR -6.400946e-19
## genotypeKO 2.573084e-01Build a data frame containing estimate, standard error, degrees of freedom, T-value and p-value for each protein in the protLM object protLMFranc, by using the test.protLMcontrast function.
resultFranc <- test.protLMcontrast(protLMFranc, L)The prot.p.adjust function adds a new column "qval" to the data frame based on an existing "pval" column. The new column contains adjusted p-values using one of the methods in the p.adjust.methods vector.
resultFranc <- prot.p.adjust(resultFranc, method="fdr")The prot.signif function needs a matrix (or a list of matrices) containing a column named “pval” and a column named "qval". The "pval" column will be used to sort the matrix. The function generates a new column "signif" containing 1 for all values of "qval" with a value lower than the specified FDR level (default: 0.05) and 0 otherwise.
resultFranc <- prot.signif(resultFranc, level = 0.05)Inspect the first 6 significant proteins:
head(resultFranc,6)## GI.number
## WP_003040894 gi|118496676
## WP_003026016 gi|118497407
## WP_003023392 gi|118497654
## WP_003040849 gi|118496659
## WP_003020026 gi|118496882
## WP_003038350 gi|118496950
## Protein.names
## WP_003040894 3-isopropylmalate dehydratase large subunit [Francisella tularensis]
## WP_003026016 biotin synthase [Francisella tularensis]
## WP_003023392 tRNA (N6-isopentenyl adenosine(37)-C2)-methylthiotransferase MiaB [Francisella tularensis]
## WP_003040849 hypothetical protein [Francisella tularensis]
## WP_003020026 CTP synthetase [Francisella tularensis]
## WP_003038350 hypothetical protein [Francisella tularensis]
## estimate se df Tval pval
## WP_003040894 -0.9057261 0.05555518 121.14412 -16.30318 2.457599e-32
## WP_003026016 -0.9617749 0.07478739 133.04725 -12.86012 4.225036e-25
## WP_003023392 -1.3965407 0.08530562 44.43458 -16.37103 1.941754e-20
## WP_003040849 -1.1356799 0.09816434 101.51013 -11.56917 3.023500e-20
## WP_003020026 -0.5163597 0.04825840 150.49382 -10.69989 3.184831e-20
## WP_003038350 0.8107740 0.06684619 72.40727 12.12895 4.127913e-19
## qval signif
## WP_003040894 2.831154e-29 TRUE
## WP_003026016 2.433621e-22 TRUE
## WP_003023392 7.337851e-18 TRUE
## WP_003040849 7.337851e-18 TRUE
## WP_003020026 7.337851e-18 TRUE
## WP_003038350 7.925594e-17 TRUEAdvanced: importing as a data frame
It is also possible to import your data as a simple data frame and make use of some custom preprocessing pipeline. Here, we will show the same preprocessing pipeline for the Francisella dataset as outlined above, but note that adding, removing or modifying steps in the preprocessing pipeline is now very simple as the data is in data frame format.
First, we import the Francisella dataset as a data frame via the base R function read.table.
pepFrancDf <- read.table(file = file_peptides_txt, sep = "\t", header = TRUE, quote="", comment.char = "")Of course, other filetypes can also be imported. Say for example your input file is a comma separated file stored in My Documents on a Windows computer, you could use for example:
#Change to eval=TRUE to execute this piece of code.
pepFrancDf <- read.table("C:/Users/Ludger/Documents/peptides.csv", sep = ",", dec=".", header = TRUE, quote="", comment.char = "")The data frame pepFrancD is in wide format.
head(pepFrancDf)## Sequence Amino.acid.before First.amino.acid
## 1 AAAEELDTR K A
## 2 AAAGFVITASHNK R A
## 3 AAANEYELALAYSIEEVAPDLHK K A
## 4 AAANNPQLEAFK K A
## 5 AAASAGLVDEK K A
## 6 AAASLDLYSYPK K A
## Second.amino.acid Second.last.amino.acid Last.amino.acid
## 1 A T R
## 2 A N K
## 3 A H K
## 4 A F K
## 5 A E K
## 6 A P K
## Amino.acid.after A.Count R.Count N.Count D.Count C.Count Q.Count E.Count
## 1 K 3 1 0 1 0 0 2
## 2 F 4 0 1 0 0 0 0
## 3 Y 6 0 1 1 0 0 4
## 4 K 4 0 2 0 0 1 1
## 5 A 4 0 0 1 0 0 1
## 6 V 3 0 0 1 0 0 0
## G.Count H.Count I.Count L.Count K.Count M.Count F.Count P.Count S.Count
## 1 0 0 0 1 0 0 0 0 0
## 2 1 1 1 0 1 0 1 0 1
## 3 0 1 1 3 1 0 0 1 1
## 4 0 0 0 1 1 0 1 1 0
## 5 1 0 0 1 1 0 0 0 1
## 6 0 0 0 2 1 0 0 1 2
## T.Count W.Count Y.Count V.Count U.Count Length Missed.cleavages
## 1 1 0 0 0 0 9 0
## 2 1 0 0 1 0 13 0
## 3 0 0 2 1 0 23 0
## 4 0 0 0 0 0 12 0
## 5 0 0 0 1 0 11 0
## 6 0 0 2 0 0 12 0
## Mass Proteins GI.number Leading.razor.protein
## 1 974.4669 WP_003038655 gi|118497196 gi|118497196
## 2 1285.6779 WP_003041237 gi|118498194 gi|118498194
## 3 2516.2435 WP_003038915 gi|118497331 gi|118497331
## 4 1272.6463 WP_003038264 gi|118496879 gi|118496879
## 5 1030.5295 WP_003035781 gi|118497152 gi|118497152
## 6 1297.6554 WP_003039212 gi|118497492 gi|118497492
## Protein.names
## 1 hypothetical protein [Francisella tularensis]
## 2 phosphoglucosamine mutase [Francisella tularensis]
## 3 glycine--tRNA ligase subunit beta [Francisella tularensis]
## 4 molecular chaperone HtpG [Francisella tularensis]
## 5 delta-aminolevulinic acid dehydratase [Francisella tularensis]
## 6 phosphorylase [Francisella tularensis]
## Unique..Groups. Unique..Proteins. Charges PEP Score
## 1 yes yes 1,2 1.2531e-04 59.35
## 2 yes yes 2 5.3020e-05 45.825
## 3 yes yes 3 1.0517e-45 84.327
## 4 yes yes 2 8.7549e-24 108.56
## 5 yes yes 2 9.5062e-05 67.385
## 6 yes yes 2 1.2563e-02 31.675
## Experiment.1WT_20_2h_n3_1 Experiment.1WT_20_2h_n3_2
## 1 NA NA
## 2 NA 1
## 3 1 1
## 4 1 1
## 5 1 NA
## 6 NA 1
## Experiment.1WT_20_2h_n3_3 Experiment.1WT_20_2h_n4_1
## 1 1 1
## 2 NA NA
## 3 2 1
## 4 1 1
## 5 1 1
## 6 1 NA
## Experiment.1WT_20_2h_n4_2 Experiment.1WT_20_2h_n4_3
## 1 NA 1
## 2 NA NA
## 3 1 1
## 4 1 1
## 5 1 NA
## 6 1 1
## Experiment.1WT_20_2h_n5_1 Experiment.1WT_20_2h_n5_2
## 1 NA 1
## 2 NA 1
## 3 2 2
## 4 1 1
## 5 1 1
## 6 NA NA
## Experiment.1WT_20_2h_n5_3 Experiment.3D8_20_2h_n3_1
## 1 NA NA
## 2 1 1
## 3 1 1
## 4 1 1
## 5 1 1
## 6 NA NA
## Experiment.3D8_20_2h_n3_2 Experiment.3D8_20_2h_n3_3
## 1 1 NA
## 2 1 NA
## 3 2 1
## 4 1 1
## 5 1 1
## 6 NA NA
## Experiment.3D8_20_2h_n4_1 Experiment.3D8_20_2h_n4_2
## 1 NA NA
## 2 1 1
## 3 1 1
## 4 1 1
## 5 1 1
## 6 NA 1
## Experiment.3D8_20_2h_n4_3 Experiment.3D8_20_2h_n5_1
## 1 NA NA
## 2 1 1
## 3 1 NA
## 4 1 1
## 5 NA 1
## 6 NA NA
## Experiment.3D8_20_2h_n5_2 Experiment.3D8_20_2h_n5_3 Intensity
## 1 NA NA 1.3295e+09
## 2 NA NA 1.6719e+08
## 3 1 1 5.5675e+08
## 4 1 1 1.4830e+10
## 5 1 1 4.8209e+09
## 6 NA NA 1.3715e+08
## Intensity.1WT_20_2h_n3_1 Intensity.1WT_20_2h_n3_2
## 1 0 0
## 2 0 20679000
## 3 28853000 28091000
## 4 339830000 420390000
## 5 240280000 0
## 6 0 31069000
## Intensity.1WT_20_2h_n3_3 Intensity.1WT_20_2h_n4_1
## 1 9290400 92996000
## 2 0 0
## 3 30436000 19476000
## 4 393930000 235060000
## 5 241890000 222250000
## 6 31448000 0
## Intensity.1WT_20_2h_n4_2 Intensity.1WT_20_2h_n4_3
## 1 0 98059000
## 2 0 0
## 3 22050000 20687000
## 4 381090000 360270000
## 5 189930000 0
## 6 27721000 28657000
## Intensity.1WT_20_2h_n5_1 Intensity.1WT_20_2h_n5_2
## 1 0 80803000
## 2 0 17358000
## 3 36588000 28654000
## 4 255710000 311370000
## 5 160440000 147230000
## 6 0 0
## Intensity.1WT_20_2h_n5_3 Intensity.3D8_20_2h_n3_1
## 1 0 0
## 2 13841000 12278000
## 3 7009100 14174000
## 4 271440000 225140000
## 5 143340000 135910000
## 6 0 0
## Intensity.3D8_20_2h_n3_2 Intensity.3D8_20_2h_n3_3
## 1 95433000 0
## 2 10712000 0
## 3 8548000 7859200
## 4 289880000 282500000
## 5 185740000 168620000
## 6 0 0
## Intensity.3D8_20_2h_n4_1 Intensity.3D8_20_2h_n4_2
## 1 0 0
## 2 9500700 11111000
## 3 14871000 11690000
## 4 209490000 223610000
## 5 160240000 158000000
## 6 0 18253000
## Intensity.3D8_20_2h_n4_3 Intensity.3D8_20_2h_n5_1
## 1 0 0
## 2 8723500 8125200
## 3 4470500 0
## 4 219700000 169120000
## 5 0 92542000
## 6 0 0
## Intensity.3D8_20_2h_n5_2 Intensity.3D8_20_2h_n5_3 Reverse Contaminant id
## 1 0 0 0
## 2 0 0 1
## 3 4048100 3422300 2
## 4 178150000 170400000 3
## 5 102350000 100090000 4
## 6 0 0 5
## Protein.group.IDs Mod..peptide.IDs
## 1 419 0
## 2 1150 1
## 3 513 2
## 4 199 3
## 5 388 4
## 6 624 5
## Evidence.IDs
## 1 0;1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18
## 2 19;20;21;22;23;24;25;26;27;28;29;30;31;32;33
## 3 34;35;36;37;38;39;40;41;42;43;44;45;46;47;48;49;50;51;52;53;54;55;56;57;58;59;60;61;62;63;64;65;66;67;68;69;70;71;72;73;74;75;76
## 4 77;78;79;80;81;82;83;84;85;86;87;88;89;90;91;92;93;94;95;96;97;98;99;100;101;102;103;104;105;106;107;108;109;110;111;112;113;114;115;116;117;118;119;120;121;122;123;124
## 5 125;126;127;128;129;130;131;132;133;134;135;136;137;138;139;140;141;142;143;144;145;146;147;148;149;150;151;152;153;154;155
## 6 156;157;158;159;160
## MS.MS.IDs
## 1 0;1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18
## 2 19;20;21;22;23;24;25;26;27;28;29;30;31;32;33
## 3 34;35;36;37;38;39;40;41;42;43;44;45;46;47;48;49;50;51;52;53;54;55;56;57;58;59;60;61;62;63;64;65;66;67;68;69;70;71;72;73;74;75;76;77;78;79
## 4 80;81;82;83;84;85;86;87;88;89;90;91;92;93;94;95;96;97;98;99;100;101;102;103;104;105;106;107;108;109;110;111;112;113;114;115;116;117;118;119;120;121;122;123;124;125;126;127;128;129;130;131;132;133;134;135;136;137;138;139;140;141;142;143;144;145;146;147;148;149;150;151;152;153;154;155;156;157;158;159;160;161;162;163;164;165;166;167;168;169;170;171;172
## 5 173;174;175;176;177;178;179;180;181;182;183;184;185;186;187;188;189;190;191;192;193;194;195;196;197;198;199;200;201;202;203
## 6 204;205;206;207;208
## Best.MS.MS Oxidation..M..site.IDs
## 1 15
## 2 26
## 3 39
## 4 159
## 5 191
## 6 205If we want to mimick the MaxQuant preprocessing pipeline, we need to add a factor that indicates whether a protein is only identified by modified peptides. This requires some basic R programming. The information about which proteins are only identified by modified peptides can be found in the proteinGroups.txt file. We again make use of the proteinGroups.txt file that is delivered with the package.
#Import the proteinGroups.txt file delivered with the MSqRobData package
file_proteinGroups <- system.file("extdata/Francisella", "proteinGroups.txt", package = "MSqRobData")
proteinGroups <- read.table(file_proteinGroups, sep="\t", header=TRUE, quote="", comment.char = "")
#Extract the column that indicates which proteins are only identified by a modification site
only_site <- proteinGroups$Only.identified.by.site
#If there are no such proteins (as is the case here), the column will be completely empty and R will import this by default as "NA". However, we want the an empty value instead of "NA" to keep this column consistent with the "Contaminant" and "Reverse" columns.
only_site[is.na(only_site)] <- ""
#Select the protein accessions that are only identified by one or more modified peptides
filter_symbol <- "+"
removed_proteins <- proteinGroups$Protein.IDs[only_site==filter_symbol]
#Create a logical variable "rem" that holds "TRUE" if a row in the pepFrancDf data frame should be removed and FALSE otherwise.
accession <- "Proteins"
rem <- as.character(pepFrancDf[,accession]) %in% as.character(removed_proteins)
#Create a new column "only_site" in "pepFrancDf" that indicates with a "+" which rows should be removed.
pepFrancDf$only_site <- ""
pepFrancDf$only_site[rem] <- "+"Note that in this particular case, the previous chunck of code was superfluous, as there are no proteins in the dataset that are only identified by modified peptides.
We can easily preprocess the data frame using the preprocess_wide function (for data in “long” format, use the preprocess_long function). This function uses our standard preprocessing pipeline on data frames in “wide” format (i.e. the data frame has one observation row per subject with each measurement present as a different variable).
We already made the exp_annotation object before (see “2.1. Create an experimental annotation data frame.”), so we will not repeat this code again here. Importantly, however, contrary to the exprs slot in the peptidesFranc2 object created before, the names of all intensity columns in data frame pepFrancDf start with “Intensity.”.
This is because the default setting remove_pattern=TRUE in the import2MSnSet function removes the prefix "Intensity\ " from the columns containing the intensities as this allows users that are unfamiliar with R to just copy and paste the names they provided in MaxQuant’s “Experiment” column as an input. Here, we just imported the peptides.txt file as it is. We thus create a new experimental annotation data frame called exp_annotation2 with run names starting with “Intensity.” (a space is not syntactically valid in R, so it is automatically converted to a dot).
runs2 <- colnames(pepFrancDf)[grepl("Intensity.",colnames(pepFrancDf))]
runs2## [1] "Intensity.1WT_20_2h_n3_1" "Intensity.1WT_20_2h_n3_2"
## [3] "Intensity.1WT_20_2h_n3_3" "Intensity.1WT_20_2h_n4_1"
## [5] "Intensity.1WT_20_2h_n4_2" "Intensity.1WT_20_2h_n4_3"
## [7] "Intensity.1WT_20_2h_n5_1" "Intensity.1WT_20_2h_n5_2"
## [9] "Intensity.1WT_20_2h_n5_3" "Intensity.3D8_20_2h_n3_1"
## [11] "Intensity.3D8_20_2h_n3_2" "Intensity.3D8_20_2h_n3_3"
## [13] "Intensity.3D8_20_2h_n4_1" "Intensity.3D8_20_2h_n4_2"
## [15] "Intensity.3D8_20_2h_n4_3" "Intensity.3D8_20_2h_n5_1"
## [17] "Intensity.3D8_20_2h_n5_2" "Intensity.3D8_20_2h_n5_3"#Extract the relevant part of the mass spec run names that points to the genotype.
genotype2 <- factor(substr(runs2, 11,13))
genotype2## [1] 1WT 1WT 1WT 1WT 1WT 1WT 1WT 1WT 1WT 3D8 3D8 3D8 3D8 3D8 3D8 3D8 3D8
## [18] 3D8
## Levels: 1WT 3D8n <- nchar(runs2)
biorep2 <- as.factor(paste0("b_",factor(as.numeric(factor(paste0(substr(runs2, 10,10),substr(runs2, n-2,n-2))))+c(rep(0,9),rep(3,9)))))
biorep## [1] b_1 b_1 b_1 b_2 b_2 b_2 b_3 b_3 b_3 b_4 b_4 b_4 b_5 b_5 b_5 b_6 b_6
## [18] b_6
## Levels: b_1 b_2 b_3 b_4 b_5 b_6exp_annotation2 <- data.frame(run=runs2, genotype=genotype2, biorep=biorep2)
exp_annotation2## run genotype biorep
## 1 Intensity.1WT_20_2h_n3_1 1WT b_1
## 2 Intensity.1WT_20_2h_n3_2 1WT b_1
## 3 Intensity.1WT_20_2h_n3_3 1WT b_1
## 4 Intensity.1WT_20_2h_n4_1 1WT b_2
## 5 Intensity.1WT_20_2h_n4_2 1WT b_2
## 6 Intensity.1WT_20_2h_n4_3 1WT b_2
## 7 Intensity.1WT_20_2h_n5_1 1WT b_3
## 8 Intensity.1WT_20_2h_n5_2 1WT b_3
## 9 Intensity.1WT_20_2h_n5_3 1WT b_3
## 10 Intensity.3D8_20_2h_n3_1 3D8 b_4
## 11 Intensity.3D8_20_2h_n3_2 3D8 b_4
## 12 Intensity.3D8_20_2h_n3_3 3D8 b_4
## 13 Intensity.3D8_20_2h_n4_1 3D8 b_5
## 14 Intensity.3D8_20_2h_n4_2 3D8 b_5
## 15 Intensity.3D8_20_2h_n4_3 3D8 b_5
## 16 Intensity.3D8_20_2h_n5_1 3D8 b_6
## 17 Intensity.3D8_20_2h_n5_2 3D8 b_6
## 18 Intensity.3D8_20_2h_n5_3 3D8 b_6The preprocess_wide function allows for the same preprocessing as the preprocess_MaxQuant function on a regular data frame. If you inspect this function, you will notice that the same 8 preprocessing steps are executed. Remember that you can do any preprocessing you like using basic R data frame manipulations. The split argument indicates which character string is used to separate the accession groups. The exp_annotation argument, like before, contains the experimental annotation data frame or a character string indicating the filepath where the experimental annotation is saved as an Excel file or a tab-delimited file. The quant_cols argument contains a vector of names or a vector of numeric indices pointing to the columns containing the quantitative values of interest (mostly intensities or areas under the curve). If quant_cols contains only one character string, each column containing this pattern will be regarded as a column containing the quantitative values of interest. aggr_by indicates the column by which the data should be aggregated. Here, the data is already aggregated to the peptide level (i.e. over different charge states and modification statuses). We need no further aggregation, thus we provide the “Sequence” column (the level to which aggregation has already been done). aggr_function will be superfluous here, but if further aggregation would have been necessary, raw intensities would be summed up (aggr_function="sum"). The other parameters are similar to the preprocess_MaxQuant function.
pepFrancDf2 <- preprocess_wide(pepFrancDf, accession="Proteins", split=";", exp_annotation=exp_annotation2, quant_cols="Intensity.", aggr_by="Sequence", aggr_function="sum", logtransform=TRUE, base=2, normalisation="quantiles", smallestUniqueGroups=TRUE, useful_properties=c("GI.number","Protein.names","Sequence"), filter=c("Contaminant","Reverse","only_site"), filter_symbol="+", minIdentified=2)Note that the experimental annotation is now present as an attribute of the pepFrancDf2 data frame.
attr(pepFrancDf2,"MSqRob_exp_annotation")## run genotype biorep
## 1 Intensity.1WT_20_2h_n3_1 1WT b_1
## 2 Intensity.1WT_20_2h_n3_2 1WT b_1
## 3 Intensity.1WT_20_2h_n3_3 1WT b_1
## 4 Intensity.1WT_20_2h_n4_1 1WT b_2
## 5 Intensity.1WT_20_2h_n4_2 1WT b_2
## 6 Intensity.1WT_20_2h_n4_3 1WT b_2
## 7 Intensity.1WT_20_2h_n5_1 1WT b_3
## 8 Intensity.1WT_20_2h_n5_2 1WT b_3
## 9 Intensity.1WT_20_2h_n5_3 1WT b_3
## 10 Intensity.3D8_20_2h_n3_1 3D8 b_4
## 11 Intensity.3D8_20_2h_n3_2 3D8 b_4
## 12 Intensity.3D8_20_2h_n3_3 3D8 b_4
## 13 Intensity.3D8_20_2h_n4_1 3D8 b_5
## 14 Intensity.3D8_20_2h_n4_2 3D8 b_5
## 15 Intensity.3D8_20_2h_n4_3 3D8 b_5
## 16 Intensity.3D8_20_2h_n5_1 3D8 b_6
## 17 Intensity.3D8_20_2h_n5_2 3D8 b_6
## 18 Intensity.3D8_20_2h_n5_3 3D8 b_6After custom preprocessing, you can transform the data to a protdata object via the df2protdata function.
The acc_col argument in the df2protdata function indicates the column index or column name of the accessions (i.e. the protein (group) identifiers). The quant_cols argument is a vector of indices pointing to all columns containing peptide intensities. In the example below, the column named “Proteins” will be added to the accession slot while all columns with indices in the quant_cols vector will be added as a matrix to the intensities slot. quant_name indicates the name that will be given to the column containing the (preprocessed) intensities and run_name indicates the name that will be given to the column containing the mass spec run names in the new protdata object. Like before, we keep “GI.number” and “Protein.names” as annotations.
#Select columns with intensities.
quant_cols <- which(grepl("Intensity.", colnames(pepFrancDf2)))
quant_cols## [1] 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22#Call df2protdata function
proteinsFrancDf <- df2protdata(pepFrancDf2, acc_col="Proteins", quant_cols=quant_cols, quant_name="quant_value", run_name="run", annotations=c("GI.number","Protein.names"))
#Inspect proteinsFrancDf protdata object.
proteinsFrancDf## @accession
##
## [1] "WP_003038655"
##
## @annotation
##
## GI.number
## "gi|118497196"
## Protein.names
## "hypothetical protein [Francisella tularensis]"
##
## @data
##
## quant_value Sequence genotype biorep run
## 2 25.11978 AVVVFFDYGCGK 1WT b_1 Intensity.1WT_20_2h_n3_1
## 3 25.34586 DLLSDPATPTVGPQDANK 1WT b_1 Intensity.1WT_20_2h_n3_1
## 4 24.11491 ETNGELTVQDVDNVVK 1WT b_1 Intensity.1WT_20_2h_n3_1
## 5 26.50674 IIAIPDIVK 1WT b_1 Intensity.1WT_20_2h_n3_1
## 6 24.69046 QTNVITDNYQQAAK 1WT b_1 Intensity.1WT_20_2h_n3_1
## 9 24.75196 AVVVFFDYGCGK 1WT b_1 Intensity.1WT_20_2h_n3_2
##
## 96 more rows...
##
## @accession
##
## [1] "WP_003041237"
##
## @annotation
##
## GI.number
## "gi|118498194"
## Protein.names
## "phosphoglucosamine mutase [Francisella tularensis]"
##
## @data
##
## quant_value Sequence genotype biorep run
## 2 25.81732 ADLGISLDGDADR 1WT b_1 Intensity.1WT_20_2h_n3_1
## 3 27.36716 FVIVGQDTR 1WT b_1 Intensity.1WT_20_2h_n3_1
## 4 24.91571 GEVANSTITAEFTQK 1WT b_1 Intensity.1WT_20_2h_n3_1
## 5 27.38539 ILANAIDEYIESIHSR 1WT b_1 Intensity.1WT_20_2h_n3_1
## 7 27.93398 SLATNEAEYLVEK 1WT b_1 Intensity.1WT_20_2h_n3_1
## 8 26.59163 VASDVNDVEK 1WT b_1 Intensity.1WT_20_2h_n3_1
##
## 187 more rows...
##
## @accession
##
## [1] "WP_003038915"
##
## @annotation
##
## GI.number
## "gi|118497331"
## Protein.names
## "glycine--tRNA ligase subunit beta [Francisella tularensis]"
##
## @data
##
## quant_value Sequence genotype biorep
## 1 24.37301 AAANEYELALAYSIEEVAPDLHK 1WT b_1
## 4 22.60222 EEGIAVDIFEAINNTNYDSIK 1WT b_1
## 5 25.34002 EGEPTQVGLGFAK 1WT b_1
## 6 25.13163 EYSYALELLTCLDK 1WT b_1
## 7 24.09977 FHHPQAIVIDNISDYIK 1WT b_1
## 9 25.75712 IIDITLESYK 1WT b_1
## run
## 1 Intensity.1WT_20_2h_n3_1
## 4 Intensity.1WT_20_2h_n3_1
## 5 Intensity.1WT_20_2h_n3_1
## 6 Intensity.1WT_20_2h_n3_1
## 7 Intensity.1WT_20_2h_n3_1
## 9 Intensity.1WT_20_2h_n3_1
##
## 215 more rows...
##
## @accession
##
## [1] "WP_003038264"
##
## @annotation
##
## GI.number
## "gi|118496879"
## Protein.names
## "molecular chaperone HtpG [Francisella tularensis]"
##
## @data
##
## quant_value Sequence genotype biorep run
## 1 27.92076 AAANNPQLEAFK 1WT b_1 Intensity.1WT_20_2h_n3_1
## 3 26.48741 EGQDTIYYITSDSYK 1WT b_1 Intensity.1WT_20_2h_n3_1
## 5 27.33299 EHLDLLEYHVLK 1WT b_1 Intensity.1WT_20_2h_n3_1
## 6 26.37656 EILQHNK 1WT b_1 Intensity.1WT_20_2h_n3_1
## 7 28.11488 ELLPPYLR 1WT b_1 Intensity.1WT_20_2h_n3_1
## 8 28.77850 ELVSNSSDAIEK 1WT b_1 Intensity.1WT_20_2h_n3_1
##
## 338 more rows...
##
## @accession
##
## [1] "WP_003035781"
##
## @annotation
##
## GI.number
## "gi|118497152"
## Protein.names
## "delta-aminolevulinic acid dehydratase [Francisella tularensis]"
##
## @data
##
## quant_value Sequence genotype biorep run
## 1 27.47327 AAASAGLVDEK 1WT b_1 Intensity.1WT_20_2h_n3_1
## 2 24.05872 VILTYTALDVANWLK 1WT b_1 Intensity.1WT_20_2h_n3_1
## 5 23.88890 VILTYTALDVANWLK 1WT b_1 Intensity.1WT_20_2h_n3_2
## 6 25.37749 YASAYYGPFR 1WT b_1 Intensity.1WT_20_2h_n3_2
## 7 27.47269 AAASAGLVDEK 1WT b_1 Intensity.1WT_20_2h_n3_3
## 8 24.93174 VILTYTALDVANWLK 1WT b_1 Intensity.1WT_20_2h_n3_3
##
## 41 more rows...
##
## @accession
##
## [1] "WP_003039212"
##
## @annotation
##
## GI.number
## "gi|118497492"
## Protein.names
## "phosphorylase [Francisella tularensis]"
##
## @data
##
## quant_value Sequence genotype biorep
## 3 22.96510 FPAGLNDVEFVAEHIFQHSK 1WT b_1
## 4 25.10257 LQDLLTFIGK 1WT b_1
## 7 24.55598 VVFSVDYR 1WT b_1
## 8 24.53369 AAASLDLYSYPK 1WT b_1
## 9 24.68029 EVLFLNR 1WT b_1
## 10 23.16038 FPAGLNDVEFVAEHIFQHSK 1WT b_1
## run
## 3 Intensity.1WT_20_2h_n3_1
## 4 Intensity.1WT_20_2h_n3_1
## 7 Intensity.1WT_20_2h_n3_1
## 8 Intensity.1WT_20_2h_n3_2
## 9 Intensity.1WT_20_2h_n3_2
## 10 Intensity.1WT_20_2h_n3_2
##
## 56 more rows...
##
## 1182 more elements...Note that the proteinsFrancDf object is the same as the proteinsFranc object, only with a different order of columns in their data slot.
