S04-quant.R

## Quantitative proteomics

##      | Label-free | Labelled    |
## -----+------------+-------------|
## MS1  | XIC        | SILAC, 15N  |
## MS2  | Counting   | TMT (iTRAQ) |

## https://lazear.github.io/sage/
## https://uclouvain-cbio.github.io/sager/

library(tidyverse)
library(msdata)
library(QFeatures)

## Quantitative assay -> SummarizedExperiment
##
## - quantitative data: matrix [features x samples]
## - samples/columns annotation (meta-data)
## - features/rows annotation (meta-data)

## - quantitative data: assay()
## - samples/columns annotation: colData()
## - features/rows annotation: rowData()

## Multiple SummarizedExperiment (PSM, peptides, proteins, ...) are
## combined into a QFeatures object.

data(feat1)
feat1

se <- feat1[[1]]
feat1[["psms"]]

se
assay(se)
rowData(se)
colData(se)

colData(feat1)

colData(feat1)$X <- c("X1", "X2")
feat1$Y <- c("Y1", "Y2")

colData(feat1)

assay(feat1)

feat1
rowData(feat1[[1]])

feat1 <- aggregateFeatures(
    feat1,
    i = "psms",
    fcol = "Sequence",
    name = "peptides",
    fun = colMeans
)

feat1

assay(feat1[[1]])
assay(feat1[[2]])

rowData(feat1[[2]])


feat1 <- aggregateFeatures(
    feat1,
    i = "peptides",
    fcol = "Protein",
    name = "proteins",
    fun = colMeans
)

feat1


feat1["ProtA", , ]

filterFeatures(feat1, ~ pval < 0.05)

filterFeatures(feat1, ~ pval < 0.05, keep = TRUE)

## As an exercise, filter rows of proteins that localise to the
## mitochondrion.

filterFeatures(feat1, ~ location == "Mitochondrion")

## Creating QFeatures objects

## data generated by PD
data(hlpsms)

hl <- readQFeatures(hlpsms,
                    ecol = 1:10,
                    name = "psms")

hl

colData(hl)

colData(hl)$group <- rep(c("CTRL", "COND"),
                         each = 5)
assay(hl[[1]])
rowData(hl[[1]])


se <- readSummarizedExperiment(hlpsms,
                               ecol = 1:10)
se

colData(se)

QFeatures(list(psms = se))

## CPTAC

## tab-delimited file generated by MQ
basename(f <- msdata::quant(full.names = TRUE))


## Read these data in as either a SummarizedExperiment (or a
## QFeatures) object.

read.delim(f) |>
    DT::datatable()


i <- 56:61
i <- grep("Intensity\\.", names(read.delim(f)))

cptac_se <- readSummarizedExperiment(f,
                                     ecol = i,
                                     fname = "Sequence",
                                     sep = "\t")

colData(cptac_se)

colnames(cptac_se) <- sub("I.+\\.", "", colnames(cptac_se))

colnames(cptac_se)

colData(cptac_se)$condition <- rep(c("6A", "6B"), each = 3)
colData(cptac_se)$id <- rep(7:9, 2)

names(rowData(cptac_se))

## --> contaminant database/fasta: CRAP

keep_vars <- c("Sequence", "Proteins", "Leading.razor.protein",
               "PEP", "Score", "Reverse", "Potential.contaminant")

rowData(cptac_se) <- rowData(cptac_se)[, keep_vars]

## Pipeline


assay(cptac_se)

anyNA(assay(cptac_se))

cptac_se <- zeroIsNA(cptac_se)

assay(cptac_se)

nNA(cptac_se)

dim(cptac_se)

barplot(nNA(cptac_se)$nNAcols$nNA)


table(nNA(cptac_se)$nNArows$nNA)

cptac_se <- filterNA(cptac_se, pNA = 4/6)

table(nNA(cptac_se)$nNArows$nNA)

nNA(cptac_se)



se1 <- impute(cptac_se, method = "knn")
assay(se1)

## Imputation with method "MinDet"
se2 <- impute(cptac_se, method = "MinDet")

## Imputation with method "zero"
se3 <- impute(cptac_se, method = "zero")

par(mfrow = c(2, 2))

plot(density(na.omit(log2(assay(cptac_se)))))
plot(density(na.omit(log2(assay(se1)))))
plot(density(na.omit(log2(assay(se2)))), title = "MinDet")
plot(density(na.omit(log2(assay(se3)+1))), title = "zero")

rowData(cptac_se) %>%
    as_tibble() %>%
    ggplot(aes(x = Score, colour = Reverse)) +
    geom_density()


rowData(cptac_se) %>%
    as_tibble() %>%
    ggplot(aes(x = PEP, colour = Reverse)) +
    geom_density()


## connected components -> also with an SE


cptac <- QFeatures(list(peptides = cptac_se))
colData(cptac) <- colData(cptac[[1]])

## Using the filterFeatures() function, filter out the reverse and
## contaminant hits, and also retain those that have a posterior error
## probability smaller than 0.05.


table(rowData(cptac_se)$Reverse)
table(rowData(cptac_se)$Potential.contaminant)

ctpac <- cptac |>
    filterFeatures(~ Reverse != "+") |>
    filterFeatures(~ Potential.contaminant != "+") |>
    filterFeatures(~ PEP < 0.05)

cptac <- cptac |>
    logTransform(i = "peptides",
                 name = "log_peptides") |>
    normalize(i = "log_peptides",
              name = "lognorm_peptides",
              method = "center.median")

## Aggregation: log-normalised peptides into proteins, defined the
## Leading.razor.protein variable.

cptac <-
    aggregateFeatures(cptac,
                      i = "lognorm_peptides",
                      name = "proteins",
                      fcol = "Leading.razor.protein",
                      fun = colMedians,
                      na.rm = TRUE)



table(rowData(cptac[[4]])$.n)

limma::plotDensities(assay(cptac[["peptides"]]))
limma::plotDensities(assay(cptac[["log_peptides"]]))
limma::plotDensities(assay(cptac[["lognorm_peptides"]]))


library(factoextra)

cptac[[3]] |>
    filterNA() |>
    assay() |>
    t() |>
    prcomp(scale = TRUE, center = TRUE) |>
    fviz_pca_ind(habillage = cptac$condition)

cptac[[3]] |>
    filterNA() |>
    assay() |>
    t() |>
    prcomp(scale = TRUE, center = TRUE) |>
    fviz_screeplot()

cptac[[4]] |>
    filterNA() |>
    assay() |>
    t() |>
    prcomp(scale = TRUE, center = TRUE) |>
    fviz_pca_ind(habillage = cptac$condition)

cptac[[4]] |>
    filterNA() |>
    assay() |>
    t() |>
    prcomp(scale = TRUE, center = TRUE) |>
    fviz_pca_ind(habillage = cptac$condition,
                 axes = c(1, 3))


cptac[[4]] |>
    filterNA() |>
    assay() |>
    t() |>
    prcomp(scale = TRUE, center = TRUE) |>
    fviz_pca(habillage = cptac$condition)


cptac["P02787ups|TRFE_HUMAN_UPS", ,
      c("lognorm_peptides", "proteins")] |>
    longFormat() |>
    as_tibble() |>
    ggplot(aes(x = colname,
               y = value,
               colour = rowname)) +
    geom_point() +
    geom_line(aes(group = rowname)) +
    facet_wrap(~ assay)

plot(cptac)


normalize(cptac, "log_peptides",
          name = "logquantiles_peptides",
          method = "quantiles") |>
    aggregateFeatures(
        "logquantiles_peptides",
        name = "proteins2",
        fcol = "Leading.razor.protein",
        fun = colMedians,
        na.rm = TRUE) |>
    plot()


se <- getWithColData(cptac, "proteins")