Revision 683ceb5e61575ab4a26fb6b8f9b001f3273cbbee authored by LizBMI on 30 May 2019, 16:26:22 UTC, committed by LizBMI on 30 May 2019, 16:26:22 UTC
2 parent s 6b96508 + ddeac21
Raw File
.Rhistory
result[[1]][,1]
colnames(result[[1]])[1]
######Run the model using R glm function (NOT runRM)##########################################
rownames(clindata_directread) == colnames(metabo_directread_fixed)
forglm <- data.frame(row.names = 1:nrow(clindata_directread))
forglm$g = rev(clindata_directread$g)
forglm$type = rev(as.factor(clindata_directread$type))
forglm$Y = as.numeric(t(metabo_directread_fixed[1,]))
lm(formula = form, data = forglm)
rownames(clindata_directread) == colnames(metabo_directread_fixed)
forglm <- data.frame(row.names = 1:nrow(clindata_directread))
forglm$g = rev(clindata_directread$g)
forglm$type = as.factor(clindata_directread$type)
forglm$Y = as.numeric(t(metabo_directread_fixed[1,]))
lm(formula = form, data = forglm)
result2[[1]]
rownames(clindata_directread) == colnames(metabo_directread_fixed)
forglm <- data.frame(row.names = 1:nrow(clindata_directread))
forglm$g = clindata_directread$g
forglm$type = as.factor(clindata_directread$type)
forglm$Y = as.numeric(t(metabo_directread_fixed[1,]))
lm(formula = form, data = forglm)
document()
install()
build()
getstatsOneLM <- function(form, clindata, arraydata) {
#array data is metabolites
#clindata is genes
call=match.call()
print(head(clindata))
print(head(arraydata))
YY <- t(arraydata)                      # the data matrix
#mean of metabolites accross all samples
EY <- apply(YY, 2, mean)                # its mean vector
#sum of squares after centering
SYY <- apply(YY, 2, function(y) {sum(y^2)}) - nrow(YY)*EY^2     # sum of squares after centering
clindata <- data.frame(y=YY[,1], clindata)
dimnames(clindata)[[2]][1] <- 'Y'
print(form)
X <- stats::model.matrix(form, clindata)       # contrasts matrix
N = dim(X)[1]
p <- dim(X)[2]
XtX <- t(X) %*% X
ixtx <- solve(XtX)
bhat <- ixtx %*% t(X) %*% YY            # Use the pseudo-inverse to estimate the parameters
yhat <- X %*% bhat                      # Figure out what is predicted by the model
# Now we partition the sum-of-square errors
rdf <- ncol(X)-1                        # number of parameters in the model
edf <- nrow(YY)-rdf-1                   # additional degrees of freedom
errors <- YY - yhat                     # difference between observed and model predictions
sse <- apply(errors^2, 2, sum)  # sum of squared errors over the samples
mse <- sse/edf                  # mean squared error
ssr <- SYY - sse                        # regression error
msr <- ssr/rdf                  # mean regression error
fval <- msr/mse                 # f-test for the overall regression
pfval <- 1-stats::pf(fval, rdf, edf)           # f-test p-values
stderror.coeff <- sapply(mse,function(x){sqrt(diag(ixtx)*x)})
t.coeff <- bhat/stderror.coeff
p.val.coeff <- 2*stats::pt(-abs(t.coeff),df = (N-p))
#methods::new('IntLimModel', call=call, model=form,
list(# call=call, model=form,
coefficients=bhat,
# predictions=yhat,
#df=c(rdf, edf),
#sse=sse,
#ssr=ssr,
#F.statistics=fval,
#F.p.values=pfval
#std.error.coeff = stderror.coeff,
#t.value.coeff = t.coeff,
p.value.coeff = p.val.coeff # interaction p-value
)
}
result2  = getstatsOneLM(form, clindata = clindata_directread, arraydata = metabo_directread_fixed)
result2
result2[[1]]
result2[[2]]
result2[[1]]
colnames(result2[[1]])[1]
result2[[1]][,1]
lm(formula = form, data = forglm)
clindata_directread
metabolite[1,]
gene = reformatCSVReadinGene(gene_directread, phenotype_directread, gene_of_interest)
metabolite = reformatCSVReadinMetabolite(metabo_directread)
document()
build()
install()
getwd()
result_getStatsOne = getCoeffUsingGetStatsOneLM(gene_directread,
metabo_directread,
phenotype_directread,
gene_of_interest = "117_at",
metabolite_of_interest = "Indoxyl" )
install.packages(devtools)
install_github("mathelab/IntLIM")
library("devtools")
install.packages(devtools)
install_github("mathelab/IntLIM")
library("IntLIM")
document()
bui
build()
install()
gene_file = "Files/Microarray/COPD136Microarray.csv"
gene_data = read.csv(gene_file,row.names=1, stringsAsFactors = FALSE)
colnames(gene_data) = unlist(lapply(colnames(gene_data), function(x) substr(x, 2, 9)))
package_path = getwd()
data_folder_path = paste0(package_path,"/Files/LCMSMetabolomicsUntargeted/")
#Column1: HILIC
HILIC_filename = "COPD131HILICPosCompoundAbundance.xlsx"
HILICFilePath = paste0(data_folder_path, HILIC_filename)
HILIC_asthma_data_tibble <- read_excel(HILICFilePath, col_names=as.character(c(1:392)))
HILIC_values_only = convertToDataFrame(HILIC_asthma_data_tibble)
#fixing errors in naming
HILICsamples = colnames(HILIC_values_only)
HILICsamples[HILICsamples=="101020"]<-"10102O"
HILICsamples[HILICsamples=="103280"]<-"10328O"
colnames(HILIC_values_only) = HILICsamples
#check errors are fixed
which(colnames(HILIC_values_only)=="10166O")
which(colnames(HILIC_values_only)=="10102O")
which(colnames(HILIC_values_only)=="10328O")
which(colnames(gene_data)=="10166O")
which(colnames(gene_data)=="10102O")
which(colnames(gene_data)=="10328O")
HILIC_incommon_with_gene = HILIC_values_only[,colnames(HILIC_values_only) %in% colnames(gene_data)]
gene_incommon_with_HILIC = gene_data[,colnames(gene_data) %in% colnames(HILIC_values_only)]
#samples only in metabolite
unique(colnames(HILIC_values_only)[which(!colnames(HILIC_values_only) %in% colnames(gene_data))])
#samples only in genes
colnames(gene_data)[which(!colnames(gene_data) %in% colnames(HILIC_values_only))]
colnames(HILIC_incommon_with_gene) = colnames(HILIC_values_only)[colnames(HILIC_values_only) %in% colnames(gene_data)]
HILIC_sample_hash = createHashofSampleToRepColumns(HILIC_incommon_with_gene)
getMeans<-function(HILIC_incommon_with_gene, HILIC_sample_hash){
#Get means of duplicates and triplicates.
HILIC_means = data.frame(row.names=1:662)
samples = c()
for(reps in keys(HILIC_sample_hash)){
for_sample_mean =  data.frame(row.names=1:662)
samples = c(samples, reps)
for (each in HILIC_sample_hash[[reps]]){
for_sample_mean = cbind(for_sample_mean, data.frame(HILIC_incommon_with_gene[,each]))
}
for_sample_mean <- data.frame(sapply(for_sample_mean, function(x) as.numeric(as.character(x))))
HILIC_means = cbind(HILIC_means, data.frame(rowMeans(for_sample_mean)))
}
colnames(HILIC_means) = samples
rownames(HILIC_means) = rownames(HILIC_incommon_with_gene)
return(HILIC_means)
}
HILIC_means = getMeans(HILIC_incommon_with_gene, HILIC_sample_hash)
HILIClog2=log2(pseudoValue(HILIC_means))
metabolite_names = data.frame(rownames(HILIClog2),stringsAsFactors = FALSE)
rownames(metabolite_names) = rownames(HILIClog2)
colnames(metabolite_names) = "id"
genes_names = as.data.frame(rownames(gene_incommon_with_HILIC))
rownames(genes_names) = rownames(gene_incommon_with_HILIC)
colnames(genes_names) = "id"
#Path to all files
package_path = getwd()
#Path to all files
metabo_data_folder_path = paste0(package_path,"/Files/LCMSMetabolomicsUntargeted/")
#Path to all files
metadata_file = paste0(metabo_data_folder_path,"COPD131-ClinicalVariables.csv")
meta_data = read.csv(metadata_file,row.names=1, stringsAsFactors = FALSE)
colnames(meta_data) = unlist(lapply(colnames(meta_data), function(x) substr(x, 2, 9)))
#check errors are fixed
which(colnames(meta_data)=="10166O")
which(colnames(meta_data)=="10102O")
which(colnames(meta_data)=="10328O")
FEV1FVC=t(meta_data[6,])
rownames(FEV1FVC) <- NULL
FEV1FVCpheno = cbind(colnames(meta_data), colnames(meta_data), FEV1FVC)
colnames(FEV1FVCpheno) = c("id", "sample", "FEV1FEC")
FEV1FVCpheno_metabolitesamples = FEV1FVCpheno[FEV1FVCpheno[,1] %in% colnames(HILIClog2),]
inputData = data.frame(c("metabolites.csv","genes.csv","metaboliteNames.csv","geneNames.csv","FEV1FVC.csv"))
rownames(inputData) = c("metabData","geneData","metabMetaData","geneMetaData","sampleMetaData")
colnames(inputData) = c("filenames")
rownames(metabolite_names)=metabolite_names[,1]
rownames(HILIClog2) == metabolite_names
getwd()
document()
library("devtooks")
library("devtools")
document()
build()
document()
build()
install()
csvfile <- "/Users/liz/COPDOmics/filesForIntLim/inputData.csv"
inputData <- IntLIM::ReadData(inputFile = csvfile,metabid='id',geneid='id')
source("R/internalfunctions.R")
csvfile <- "/Users/liz/COPDOmics/filesForIntLim/inputData.csv"
inputData <- IntLIM::ReadData(inputFile = csvfile,metabid='id',geneid='id')
myres1 <- RunIntLim(inputData = inputData,stype="FEV1FEC", continuous = TRUE)
mpca = computePCA(metabolite)
library("COPDOmics")
mpca = computePCA(metabolite)
metabo_directread = read.csv("filesForIntLim/metabolites.csv", header=FALSE, colClasses = "character")
metabolite = reformatCSVReadinMetabolite(metabo_directread)
mpca = computePCA(metabolite)
metabo_directread = read.csv("/Users/liz/COPDOmics/filesForIntLim/metabolites.csv", header=FALSE, colClasses = "character")
metabolite_directread
metabo_directread
metabolite = reformatCSVReadinMetabolite(metabo_directread)
mpca = computePCA(metabolite)
mpca
metabo_directread = read.csv("/Users/liz/COPDOmics/filesForIntLim/metabolites.csv", header=FALSE, colClasses = "character")
metabolite = reformatCSVReadinMetabolite(metabo_directread)
mpca = computePCA(metabolite)
myres1 <- RunIntLim(inputData = inputData,stype="FEV1FEC", continuous = TRUE, covar=c(mpca$x[,1],mpca$x[,2]), class.covar="numeric")
createAsthmaHILICCsvForProcessing(subset = TRUE)
source("R/internalfunctions.R")
csvfile <- "/Users/liz/COPDOmics/filesForIntLim/inputData.csv"
inputData <- IntLIM::ReadData(inputFile = csvfile,metabid='id',geneid='id')
IntLIM::ShowStats(IntLimObject = inputData)
incommon <- getCommon(inputData,stype="FEV1FEC")
metabo_directread = read.csv("/Users/liz/COPDOmics/filesForIntLim/metabolites.csv", header=FALSE, colClasses = "character")
metabolite = reformatCSVReadinMetabolite(metabo_directread)
mpca = computePCA(metabolite)
myres1 <- RunIntLim(inputData = inputData,stype="FEV1FEC", continuous = TRUE, covar=c(mpca$x[,1],mpca$x[,2]), class.covar="numeric")
myres1 <- RunIntLim(inputData = inputData,stype="FEV1FEC", continuous = TRUE, covar=c(mpca$x[,1],mpca$x[,2]), class.covar=c("numeric","numeric"))
mpca$x[,1]
mpca
myres1 <- RunIntLim(inputData = inputData,stype="FEV1FEC", continuous = TRUE)
incommon <- getCommon(inputData,stype="FEV1FEC")
incommon
stype="FEV1FEC"
covar=c(mpca$x[,1],mpca$x[,2])
class.covar=c("numeric","numeric")
incommon<-MultiDataSet::commonSamples(inputData)
mp <- Biobase::pData(incommon[["metabolite"]])
gp <- Biobase::pData(incommon[["expression"]])
if(all.equal(mp[,stype],gp[,stype])[1] != TRUE) {
stop(paste("The column", stype,"for the samples in common between the metabolite and gene datasets are not equal.  Please check your input."))
}
gene <- Biobase::assayDataElement(inputData[["expression"]], 'exprs')
metab <- Biobase::assayDataElement(inputData[["metabolite"]], 'metabData')
# Force the order to be the same, in case it isn't
p <- mp[rownames(gp),]
p.0 <- p
metab <- metab[,colnames(gene)]
if(!is.null(stype)) {
p <- p[,stype]
uniqp <- unique(p)
uniqtypes <- unique(p)
# Deal with missing values or ""
if(length(which(p==""))>0) {
new.p <- p[which(p!="")]
metab <- metab[,which(p!="")]
gene <- gene[,which(p!="")]
p <- new.p
}
if(length(which(is.na(p)))>0) {
new.p <- p[which(!is.na(p))]
metab <- metab[,which(!is.na(p))]
gene <- gene[,which(!is.na(p))]
p <- new.p
}
}
if(!is.null(covar)){
if (length(covar %in% colnames(p.0)) != sum(covar %in% colnames(p.0))){
stop("Additional variable names not in pData")
}
covar_matrix <- p.0[colnames(gene),covar, drop = FALSE]
na.covar <- which(is.na(covar_matrix) | covar_matrix == '',arr.ind = TRUE)
na.covar.list <- unique(rownames(na.covar))
new.overall.list <- setdiff(colnames(gene), na.covar.list)
covar_matrix <- covar_matrix[new.overall.list,,drop = FALSE]
class.var <- apply(covar_matrix,2,class)
gene <- gene[,new.overall.list]
metab <- metab[,new.overall.list]
p <- p.0[new.overall.list,stype]
if(!(is.null(class.covar))){
if(length(class.covar) != length(covar)){
stop("lengths of covar and class.covar not the same")
}
len.covar <- length(covar)
for(i in 1:len.covar){
if(class.covar[i] == 'numeric'){
covar_matrix[,i] <- as.numeric(covar_matrix[,i])
}else{
covar_matrix[,i] <- as.factor(as.character(covar_matrix[,i]))
}
}
}
}
covar
colnames(p.0)
covar
covar
colnames(p.0)
covar
colnames(p.0)
p <- mp[rownames(gp),]
p
p.0
rownames(p.0)
covar %in% colnames(p.0)
covar %in% rownames(p.0)
covar
names(covar)
names(covar) %in% rownames()
names(covar) %in% rownames(p.0)
sum(covar) %in% rownames(p.0)
sum(covar)
gp
incommon <- getCommon(inputData,stype="FEV1FEC")
incommon
myres1 <- RunIntLim(inputData = inputData,stype="FEV1FEC", continuous = TRUE)
covar
colnames(p.0)
p.0
colnames(p.0)
covar
(length(names(covar) %in% rownames(p.0)) != sum(names(covar) %in% rownames(p.0))
)
sum(names(covar) %in% rownames(p.0))
length(names(covar) %in% rownames(p.0))
getCommon <- function(inputData,stype=NULL, covar = NULL, class.covar = NULL) {
incommon<-MultiDataSet::commonSamples(inputData)
mp <- Biobase::pData(incommon[["metabolite"]])
gp <- Biobase::pData(incommon[["expression"]])
if(all.equal(mp[,stype],gp[,stype])[1] != TRUE) {
stop(paste("The column", stype,"for the samples in common between the metabolite and gene datasets are not equal.  Please check your input."))
}
gene <- Biobase::assayDataElement(inputData[["expression"]], 'exprs')
metab <- Biobase::assayDataElement(inputData[["metabolite"]], 'metabData')
# Force the order to be the same, in case it isn't
p <- mp[rownames(gp),]
p.0 <- p
metab <- metab[,colnames(gene)]
if(!is.null(stype)) {
p <- p[,stype]
uniqp <- unique(p)
uniqtypes <- unique(p)
# Deal with missing values or ""
if(length(which(p==""))>0) {
new.p <- p[which(p!="")]
metab <- metab[,which(p!="")]
gene <- gene[,which(p!="")]
p <- new.p
}
if(length(which(is.na(p)))>0) {
new.p <- p[which(!is.na(p))]
metab <- metab[,which(!is.na(p))]
gene <- gene[,which(!is.na(p))]
p <- new.p
}
}
if(!is.null(covar)){
if (length(names(covar) %in% rownames(p.0)) != sum(names(covar) %in% rownames(p.0))){
stop("Additional variable names not in pData")
}
covar_matrix <- p.0[colnames(gene),covar, drop = FALSE]
na.covar <- which(is.na(covar_matrix) | covar_matrix == '',arr.ind = TRUE)
na.covar.list <- unique(rownames(na.covar))
new.overall.list <- setdiff(colnames(gene), na.covar.list)
covar_matrix <- covar_matrix[new.overall.list,,drop = FALSE]
class.var <- apply(covar_matrix,2,class)
gene <- gene[,new.overall.list]
metab <- metab[,new.overall.list]
p <- p.0[new.overall.list,stype]
if(!(is.null(class.covar))){
if(length(class.covar) != length(covar)){
stop("lengths of covar and class.covar not the same")
}
len.covar <- length(covar)
for(i in 1:len.covar){
if(class.covar[i] == 'numeric'){
covar_matrix[,i] <- as.numeric(covar_matrix[,i])
}else{
covar_matrix[,i] <- as.factor(as.character(covar_matrix[,i]))
}
}
}
}else{
covar_matrix <- NULL
}
# Check that everything is in right order
if(!all.equal(rownames(mp),rownames(gp)) || !all.equal(colnames(metab),colnames(gene))){
stop("Something went wrong with the merging!  Sample names of input files may not match.")
} else {
out <- list(p=as.factor(as.character(p)),gene=gene,metab=metab, covar_matrix=covar_matrix)
}
return(out)
}
incommon <- getCommon(inputData,stype="FEV1FEC", covar=c(mpca$x[,1],mpca$x[,2]), class.covar=c("numeric","numeric"))
incommon <- getCommon(inputData,stype="FEV1FEC")
if (length(names(covar) %in% rownames(p.0)) != sum(names(covar) %in% rownames(p.0))){
stop("Additional variable names not in pData")
}
covar_matrix <- p.0[colnames(gene),covar, drop = FALSE]
covar_matrix <- p.0[colnames(gene),covar, drop = FALSE]
colnames(gene)
covar
document()
build()
install()
length(covar %in% colnames(p.0)) != sum(covar %in% colnames(p.0))
colnames(p.0)
rownames(p.0)
covar
dir <- system.file("extdata", package="IntLIM", mustWork=TRUE)
dir
covar
colnames(p.0)
covar
covar = c("FEV1FEC")
class.covar = c("FEV1FEC")
class.covar = c("numeric")
incommon<-MultiDataSet::commonSamples(inputData)
mp <- Biobase::pData(incommon[["metabolite"]])
gp <- Biobase::pData(incommon[["expression"]])
if(all.equal(mp[,stype],gp[,stype])[1] != TRUE) {
stop(paste("The column", stype,"for the samples in common between the metabolite and gene datasets are not equal.  Please check your input."))
}
gene <- Biobase::assayDataElement(inputData[["expression"]], 'exprs')
metab <- Biobase::assayDataElement(inputData[["metabolite"]], 'metabData')
# Force the order to be the same, in case it isn't
p <- mp[rownames(gp),]
p.0 <- p
metab <- metab[,colnames(gene)]
if(!is.null(stype)) {
p <- p[,stype]
uniqp <- unique(p)
uniqtypes <- unique(p)
# Deal with missing values or ""
if(length(which(p==""))>0) {
new.p <- p[which(p!="")]
metab <- metab[,which(p!="")]
gene <- gene[,which(p!="")]
p <- new.p
}
if(length(which(is.na(p)))>0) {
new.p <- p[which(!is.na(p))]
metab <- metab[,which(!is.na(p))]
gene <- gene[,which(!is.na(p))]
p <- new.p
}
}
if (length(covar %in% colnames(p.0)) != sum(covar %in% colnames(p.0))){
stop("Additional variable names not in pData")
}
documnent()
library("devtools")
document()
build()
install()
myres2 <- RunIntLim(inputData = inputData,stype="FEV1FEC", continuous = TRUE, covar=c("PC1", "PC2"), class.covar=c("numeric","numeric"))
myres1 <- RunIntLim(inputData = inputData,stype="FEV1FEC", continuous = TRUE)
myres2 <- RunIntLim(inputData = inputData,stype="FEV1FEC", continuous = TRUE, covar=c("PC1", "PC2"), class.covar=c("numeric","numeric"))
setwd("/Users/liz/COPDOmics")
samples = createCSVForIntLIMGenesMetabolites(subset = TRUE, metabolite_file_name="COPD131HILICPosCompoundAbundance.xlsx", output_file_prefix = "HILIC")
gene_directread = read.csv("filesForIntLim/HILICgenes.csv", header=FALSE)
phenotype_directread = read.csv("filesForIntLim/HILICFEV1FEC.csv")
metabo_directread = read.csv("filesForIntLim/HILICmetabolites.csv", header=FALSE, colClasses = "character")
metabolite = reformatCSVReadinMetabolite(metabo_directread)
package_path = "/Users/liz/COPDOmics"
batch_path = paste0(package_path,"/Files/LCMSMetabolomicsUntargeted/COPD_MS_Batch_Info.xlsx")
batch_tibble <- read_excel(batch_path, col_names=as.character(c(1:2)))
batch_dataframe = convertToDataFrame(batch_tibble)
metabolite_label_vector = createLabelVector(batch_dataframe, metabolite)
metabolite_color_vector = createColorVector(metabolite_label_vector)
mpca = computePCA(metabolite)
graphPCA(mpca, "PCA by batch", metabolite_color_vector, metabolite_label_vector)
createCSVForIntLIMPhenotype(output_file_prefix = "HILIC", pca = mpca, samples = samples)
source("R/internalfunctions.R")
csvfile <- "/Users/liz/COPDOmics/filesForIntLim/HILICinputData.csv"
inputData <- IntLIM::ReadData(inputFile = csvfile,metabid='id',geneid='id')
myres1 <- RunIntLim(inputData = inputData,stype="FEV1FEC", continuous = TRUE)
myres2 <- RunIntLim(inputData = inputData,stype="FEV1FEC", continuous = TRUE, covar=c("PC1", "PC2"), class.covar=c("numeric","numeric"))
setwd("/Users/liz/COPDOmics")
samples = createCSVForIntLIMGenesMetabolites(subset = TRUE, metabolite_file_name="COPD131HILICPosCompoundAbundance.xlsx", output_file_prefix = "HILIC")
gene_directread = read.csv("filesForIntLim/HILICgenes.csv", header=FALSE)
phenotype_directread = read.csv("filesForIntLim/HILICFEV1FEC.csv")
metabo_directread = read.csv("filesForIntLim/HILICmetabolites.csv", header=FALSE, colClasses = "character")
metabolite = reformatCSVReadinMetabolite(metabo_directread)
package_path = "/Users/liz/COPDOmics"
batch_path = paste0(package_path,"/Files/LCMSMetabolomicsUntargeted/COPD_MS_Batch_Info.xlsx")
batch_tibble <- read_excel(batch_path, col_names=as.character(c(1:2)))
batch_dataframe = convertToDataFrame(batch_tibble)
library("readxl")
library("VennDiagram")
library("highcharter")
library("dplyr")
library("reshape2")
library("hash")
library("RColorBrewer")
library("IntLIM")
library("COPDOmics")
setwd("/Users/liz/COPDOmics")
samples = createCSVForIntLIMGenesMetabolites(subset = TRUE, metabolite_file_name="COPD131HILICPosCompoundAbundance.xlsx", output_file_prefix = "HILIC")
gene_directread = read.csv("filesForIntLim/HILICgenes.csv", header=FALSE)
phenotype_directread = read.csv("filesForIntLim/HILICFEV1FEC.csv")
metabo_directread = read.csv("filesForIntLim/HILICmetabolites.csv", header=FALSE, colClasses = "character")
metabolite = reformatCSVReadinMetabolite(metabo_directread)
package_path = "/Users/liz/COPDOmics"
batch_path = paste0(package_path,"/Files/LCMSMetabolomicsUntargeted/COPD_MS_Batch_Info.xlsx")
batch_tibble <- read_excel(batch_path, col_names=as.character(c(1:2)))
batch_dataframe = convertToDataFrame(batch_tibble)
metabolite_label_vector = createLabelVector(batch_dataframe, metabolite)
metabolite_color_vector = createColorVector(metabolite_label_vector)
mpca = computePCA(metabolite)
graphPCA(mpca, "PCA by batch", metabolite_color_vector, metabolite_label_vector)
createCSVForIntLIMPhenotype(output_file_prefix = "HILIC", pca = mpca, samples = samples)
source("R/internalfunctions.R")
csvfile <- "/Users/liz/COPDOmics/filesForIntLim/HILICinputData.csv"
inputData <- IntLIM::ReadData(inputFile = csvfile,metabid='id',geneid='id')
myres1 <- RunIntLim(inputData = inputData,stype="FEV1FEC", continuous = TRUE)
myres2 <- RunIntLim(inputData = inputData,stype="FEV1FEC", continuous = TRUE, covar=c("PC1", "PC2"), class.covar=c("numeric","numeric"))
back to top