Set01 Analysis

suppressMessages(library(limma))
suppressMessages(library(qvalue))
suppressMessages(library(data.table))
suppressMessages(library(cowplot))
suppressMessages(library(tidyverse))
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

dat <- fread('../data/ALA_Set01_All_Fractions_Peptide_Abundances.csv') 
colnames(dat)

 [1] "V1"                                            "Protein FDR Confidence"                        "Master"                                        "Accession"                                    
 [5] "Description"                                   "Exp. q-value"                                  "Sum PEP Score"                                 "Coverage"                                     
 [9] "# Peptides"                                    "# PSMs"                                        "# Unique Peptides"                             "# Protein Groups"                             
[13] "# AAs"                                         "MW [kDa]"                                      "calc. pI"                                      "Found in Sample: [S1] F1: 114, Control"       
[17] "Found in Sample: [S2] F1: 115, Sample"         "Found in Sample: [S3] F1: 116, Sample"         "Found in Sample: [S4] F1: 117, Sample"         "Modifications"                                
[21] "Abundance Ratio: (F1, 115) / (F1, 114)"        "Abundance Ratio: (F1, 116) / (F1, 114)"        "Abundance Ratio: (F1, 117) / (F1, 114)"        "Abundance Ratio (log2): (F1, 115) / (F1, 114)"
[25] "Abundance Ratio (log2): (F1, 116) / (F1, 114)" "Abundance Ratio (log2): (F1, 117) / (F1, 114)" "Abundances (Scaled): F1: 114, Control"         "Abundances (Scaled): F1: 115, Sample"         
[29] "Abundances (Scaled): F1: 116, Sample"          "Abundances (Scaled): F1: 117, Sample"          "Abundance: F1: 114, Control"                   "Abundance: F1: 115, Sample"                   
[33] "Abundance: F1: 116, Sample"                    "Abundance: F1: 117, Sample"                    "emPAI"                                         "# Razor Peptides"                             
[37] "Score Sequest HT"                              "# Peptides Sequest HT"                        

dat <- subset(dat, V1==TRUE)
setnames(dat, 
         old = c('V1', 
                 'Abundance: F1: 114, Control', 
                 'Abundance: F1: 115, Sample', 
                 'Abundance: F1: 116, Sample',
                 'Abundance: F1: 117, Sample'), 
         new = c('isProtein', 
                 'X114',
                 'X115',
                 'X116',
                 'X117'))
dat <- dat[dat$isProtein == TRUE,]
dat <- dat[, c('Accession', 'Description', 'Exp. q-value', 
                'Coverage', '# Peptides', '# PSMs', '# Unique Peptides', 
                '# Protein Groups', '# AAs', 'X114', 'X115', 'X116', 'X117')]
dat <-  transform(dat, 
                   X114 = as.numeric(X114),
                   X115 = as.numeric(X115),
                   X116 = as.numeric(X116),
                   X117 = as.numeric(X117)
                   )
dat

dat.melt <- melt(dat)

Intensity plots

Already normalized?

ggplot(dat.melt, aes(x=variable, y=value, fill=factor(variable))) + 
  geom_boxplot(position = "dodge") + 
  scale_fill_manual(values = c(cbPalette[2],cbPalette[4],cbPalette[6],cbPalette[7])) + 
  scale_y_log10()+ labs(fill='Channel') + 
  ggtitle('Intensitites') +
  xlab('Channel') + 
  ylab('log2 Intensity')
ggsave('Internsities.pdf')

Saving 7.29 x 4.5 in image

Moderated t-tests

design <- model.matrix(~factor(c(1,2,3,4)))
colnames(design) <- c("Intercept", "C2vsC1", "C3vsC1", "C4vsC1")
datx <- dat[, c('X114', 'X115', 'X116', 'X117')]
fit <- lmFit(datx, design)

Partial NA coefficients for 6 probe(s)

fit

An object of class "MArrayLM"
$coefficients
     Intercept    C2vsC1    C3vsC1    C4vsC1
[1,]  181159.9 -140399.4 -124813.0 -144177.2
[2,]      73.6     -55.8     -60.5     -60.5
[3,]   37028.4    9143.8   26218.4  -15385.9
[4,]   48260.4   -7875.1    3483.1  -25204.9
[5,]   85650.9  -61177.2  -57696.0  -64848.1
3260 more rows ...

$stdev.unscaled
     Intercept   C2vsC1   C3vsC1   C4vsC1
[1,]         1 1.414214 1.414214 1.414214
[2,]         1 1.414214 1.414214 1.414214
[3,]         1 1.414214 1.414214 1.414214
[4,]         1 1.414214 1.414214 1.414214
[5,]         1 1.414214 1.414214 1.414214
3260 more rows ...

$sigma
[1] NA NA NA NA NA
3260 more elements ...

$df.residual
[1] 0 0 0 0 0
3260 more elements ...

$cov.coefficients
          Intercept C2vsC1 C3vsC1 C4vsC1
Intercept         1     -1     -1     -1
C2vsC1           -1      2      1      1
C3vsC1           -1      1      2      1
C4vsC1           -1      1      1      2

$pivot
[1] 1 2 3 4

$rank
[1] 4

$Amean
[1] 78812.50    29.40 42022.47 40861.18 39720.57
3260 more elements ...

$method
[1] "ls"

$design
  Intercept C2vsC1 C3vsC1 C4vsC1
1         1      0      0      0
2         1      1      0      0
3         1      0      1      0
4         1      0      0      1
attr(,"assign")
[1] 0 1 1 1
attr(,"contrasts")
attr(,"contrasts")$`factor(c(1, 2, 3, 4))`
[1] "contr.treatment"