inst/doc/MultiDataSet_Extending_Proteome.R

In MultiDataSet: Implementation of MultiDataSet and ResultSet

## ----create_data, eval = TRUE, echo = FALSE-----------------------------------
#setwd("/home/kuragari/Downloads/use case")
#setwd("C:/Users/chernandez/shared/Projects/multidataset/Tutorials")
fdata <- data.frame(
    feature = c("Adiponectin", "CRP", "APO.A1", "APO.B", "APO.E"),
    LoD.T = c(4.0000E+07, 6.8333E+06, 1.6250E+09, 7.5000E+08, 3.0556E+08),
    LoD.B = c(1.6461E+05, 2.8121E+04,	2.0062E+07,	9.2593E+06,	3.7723E+06),
    unit = c("pg/mL", "pg/mL", "pg/mL", "pg/mL", "pg/mL")
)


pc <- textConnection("sample,gender,plate,kit
sp001,male,P01,5
sp002,male,P01,5
sp003,female,P01,5
sp004,male,P01,5
sp005,female,P01,5
sp006,female,P01,5
sp007,male,P01,5
sp009,male,P01,5
sp010,male,P01,5
sp012,male,P01,5
sp013,male,P01,5
sp014,female,P01,5
sp015,male,P01,5
sp016,male,P01,5
sp017,male,P01,5
sp018,male,P01,5
sp019,male,P01,5
sp020,male,P01,5
sp021,female,P01,5
sp022,male,P01,5
sp023,male,P01,5
sp024,female,P01,5
sp025,male,P01,5
sp026,male,P01,5
sp027,female,P01,5
sp028,male,P01,5
sp029,female,P01,5
sp030,female,P01,5
sp031,male,P01,5
sp032,female,P01,5
sp033,male,P01,5
sp034,male,P01,5
sp035,female,P01,5")

pdata <- read.csv(pc)

ac <- textConnection("sample,Adiponectin,CRP,APO.A1,APO.B,APO.E
sp001,15078717.84,256079.1978,166102189.6,81989413.94,12088576.71
sp002,13216433.82,59688.53882,122910159.9,75276925.92,8662032.005
sp003,25639817.37,272556.6707,132103162.9,58156884.1,11317489.32
sp004,28042622.37,981180.998,163470320.6,72249107.96,9534858.626
sp005,37507509.6,1434113.656,186415820,85977072.27,20772027.96
sp006,27201037.83,7140776.893,180018819.1,115016660.7,19020968.34
sp007,19734082.42,849104.2897,95178150.98,74531939.7,12575510.74
sp009,24337465.91,2083163.005,119909512.6,75162215.53,13472623.87
sp010,29346604.23,387120.4506,122159943.6,79429725.86,16315028.32
sp012,23471847.37,313730.4201,162530447.6,61607843.14,6950787.103
sp013,26464829.9,231195.7943,182276383.6,91789932.63,10467398.75
sp014,15188363.38,13795018.08,179266346.9,142227471.5,10822788.47
sp015,28585011.4,48724.14524,182652665.8,98888608.08,7252478.096
sp016,19105671.52,110506.2048,140362398.3,64931137.97,8072973.965
sp017,26570792.93,991854.9577,145432469.7,80533037.64,16914321.56
sp018,23968478.51,234826.4158,127787432.4,65209657.88,7924737.258
sp019,22881995.44,9321008.038,144681259.3,111528556.5,8072973.965
sp020,23118772.5,16891530.05,196202782.9,80823891.6,14734492.79
sp021,22942701.12,139509.27,92932170.07,59302740.96,4131144.701
sp022,25317279.56,103616.9905,115034724.1,75678689.62,7252478.096
sp023,26502576.75,79198.05423,209006876.8,110535951,20643082.4
sp024,25540633.93,334148.4346,182276383.6,84506590.67,14632785.98
sp025,34258505.11,654613.9767,115034724.1,90954853.78,18563191.06
sp026,17832612.2,229986.4403,111660816.1,61993740.5,10964463.6
sp027,20439508.67,534718.8069,157455877.1,56311320.18,7327611.264
sp028,31333461.27,179608.5461,113535112.2,69359486.56,12852515.92
sp029,18419548.11,271333.7529,100420374.9,41698956.49,10110231.2
sp030,15699158.6,1011889.525,147310628.7,54048254.2,8735244.722
sp031,20726453.63,479871.0422,98922372.45,58156884.1,10503016.43
sp032,33833352.39,532158.9756,214657722.7,84800277.57,19866893.4
sp033,14598125.04,179608.5461,84887644.72,86212822.7,13952043.6
sp034,24457846.67,1219055.651,124410700.1,85918155,18235193.56
sp035,11560113.29,134838.8042,92557881.06,34647582.51,6950787.103")

adata <- read.csv(ac)
rm(pc, ac)


write.table(fdata, file="feature_data.tsv", sep="\t", row.names=FALSE, col.names=TRUE, quote=FALSE)
write.table(pdata, file="pheno_data.tsv", sep="\t", row.names=FALSE, col.names=TRUE, quote=FALSE)
write.table(adata, file="assay_data.tsv", sep="\t", row.names=FALSE, col.names=TRUE, quote=FALSE)
rm(fdata, pdata, adata)

## ----Define Proteome, message=FALSE-------------------------------------------
library(Biobase)
setClass (
    Class = "ProteomeSet",
    contains = "eSet", 
    prototype = prototype(new("VersionedBiobase",
                              versions = c(classVersion("eSet"), ProteomeSet = "1.0.0")))
)

## ----Define Proteome Validity, message=FALSE, results="hide"------------------
setValidity("ProteomeSet", function(object) {
    
    ## Check that object has the slots 'prots' and 'raw' in assayData
    msg <- validMsg(NULL, assayDataValidMembers(assayData(object), c("prots", "raw")))
    
    ## Check that objects has the columns 'LoD.T' and 'LoD.B' in featureData
    msg <- validMsg(msg, ifelse(any(!c("LoD.T", "LoD.B") %in% fvarLabels(object)), "fData must contain columns 'LoD.T' and 'LoD.B'", FALSE))
    if (is.null(msg)){
        TRUE
    } else{
        msg
    }
}) 

## ----read_ldset---------------------------------------------------------------
read_ldset <- function(assayFile, phenoFile, featureFile, sep="\t") {
    ## Load the threes files that will be used to create the ProteomeSet
    adata <- read.delim(assayFile, sep=sep, header=TRUE, row.names="sample")
    pdata <- read.delim(phenoFile, sep=sep, header=TRUE, row.names="sample")
    fdata <- read.delim(featureFile, sep=sep, header=TRUE, row.names="feature")
    ## /
    
    ## Check that proteins in assay data are the same in feature data
    if(!identical(colnames(adata), rownames(fdata))) {
        stop("Features names in assay data (columns) are not equal to ",
             "features names in feature data (rownames).")
    }
    ##/
    
    ## Check that feature data include columns LoD.B and LoD.T
    if(sum(c("LoD.T", "LoD.B") %in% colnames(fdata)) != 2) {
        stop("Feature data must contain two columns labeled 'LoD.T' (top ",
             "limit of dectection) and 'LoD.B (bottom limit of dectection)")
    }
    ## /
    
    ## Perform the transformation of the protein level of expression
    low <- fdata[colnames(adata), "LoD.B"]
    up <- fdata[colnames(adata), "LoD.T"]
    names(low) <- names(up) <- rownames(fdata)
    faux <- function(x, low, up) {
        x[x < low] <- as.double(low / 2)
        x[x > up] <- as.double(up * 1.5)
        x
    }
    tadata <- mapply(FUN = faux, x = as.list(adata), low = as.list(low), up = as.list(up))
    dimnames(tadata) <- dimnames(adata)
    ## /
    
    ## Create the ExpressionSet with the two matrices
    prot <- new("ProteomeSet",
                assayData = assayDataNew("environment", prots = t(tadata), raw = t(adata)),
                phenoData = AnnotatedDataFrame(pdata),
                featureData = AnnotatedDataFrame(fdata)
    )
    ## /
    
    ## Check that the new ProteomeSet is valid
    validObject(prot)
    ## /
    
    return(prot)
}

## ----create_generic-----------------------------------------------------------
setGeneric("add_prot", function(object, protSet, warnings = TRUE, ...)
    standardGeneric("add_prot")
)

## ---- message = FALSE---------------------------------------------------------
library(MultiDataSet)
setMethod(
    f = "add_prot",
    signature = c("MultiDataSet", "ProteomeSet"),
    definition = function(object, protSet, warnings = TRUE, ...) {
        ## Add given ProteomeSet as 'proteome'
        object <- MultiDataSet::add_eset(object, protSet, dataset.type = "proteome", GRanges = NA, ...)
        ## /
        return(object)
    }
)

## ----load_ps------------------------------------------------------------------
## Create a ProteomeSet with protein data
ps <- read_ldset(assayFile="assay_data.tsv", 
                 phenoFile="pheno_data.tsv",
                 featureFile="feature_data.tsv"
)
ps

## ----add_prot_md_1, warning=FALSE---------------------------------------------
md <- createMultiDataSet()
md <- add_prot(md, ps)

## ----show_md_1----------------------------------------------------------------
names(md)

## ----show_md_2----------------------------------------------------------------
md