-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathcode_for_SE.R
More file actions
169 lines (115 loc) · 5.29 KB
/
Copy pathcode_for_SE.R
File metadata and controls
169 lines (115 loc) · 5.29 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
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
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
# FILE: lipidomics_analysis_pipeline_psoriasis.R
# DATE: 11 October 2022
# AUTHOR: Nisha Stephan
# PURPOSE: reading metabolon lipidomics data and processing it in autonomics and/or maplet
# This is setup under github project https://github.com/wcmq-abmc/Psoriasis
# cleanup
rm(list=ls())
# libraries
library(tidyverse)
library(autonomics)
library(maplet)
library(magrittr)
library(grid)
library(gridExtra)
library(data.table)
library(rstatix)
library(ggpubr)
library(ggplot2)
library(ggdendro)
library(dplyr)
library(cowplot)
library(plotly)
#library(xlsx)
source("mt_plots_volcano_log_pvalue.R")
source("mt_plots_pls.R")
source("mt_plots_boxplot.R")
source("mt_reporting_html_resize.R")
source("mt_plots_volcano_log_pvalue.R")
source("mt_plots_box_scatter_test.R")
old <- theme_set(theme_bw())
theme_set(theme_bw())
# define basename for output
basename = "lipidomics_analysis_psoriasis"
subgroupvar = 'GROUP_ID'
# open a log file
logfile = paste(basename, ".log", sep = "")
options(warn=-1); sink(); options(warn=0)
sink(file = logfile, split = TRUE)
cat("Start logging to", logfile, "\n")
print(Sys.time())
######################################################################################################
# # load the data from excel file "WCQA-01-22MLCLP PLASMA CLP DATA TABLES.xlsx"
######################################################################################################
file_lipidomics_data <- "WCQA-01-22MLCLP-PLASMA-CLP-DATA-TABLES.xlsx"
pheno_data <- "pheno_data.xlsx"
# To get the checksum
# library(openssl)
# as.character(md5(file(file_lipidomics_data, open = "rb")))
D_CLP <-
# validate checksum
mt_load_checksum(file=file_lipidomics_data, checksum = "83c4beae42b164f609faf6218485c7ec") %>%
# load data
mt_load_metabolon_lipidomics(file=file_lipidomics_data, sheet_list = c("Lipid Class Concentrations", "Species Concentrations","Fatty Acid Concentrations")) %>%
# # log assay dimensions and number of columns for both metabolite and clincial annotations
mt_reporting_data() %>%
# start timing
mt_reporting_tic() %>%
{.}
# the make.names function adds X to the rownames
rowData(D_CLP)$CHEMICAL_NAME = make.names(rowData(D_CLP)$name)
rownames(D_CLP) = rowData(D_CLP)$CHEMICAL_NAME
# Requirements for autonomics
rowData(D_CLP)$name = rownames(D_CLP)
rowData(D_CLP)$BIOCHEMICAL = rownames(D_CLP)
rowData(D_CLP)$feature_id = rownames(D_CLP)
colData(D_CLP)$sample_id = colData(D_CLP)$CLIENT_SAMPLE_ID
colnames(D_CLP) = colData(D_CLP)$CLIENT_SAMPLE_ID
# Read additional phenotypes
pheno = readxl::read_excel(path = pheno_data, sheet = "Sheet1", col_names = TRUE,
na = c("", "NA"))
pheno$Age = pheno$`Age (in years)`
pheno$`Age (in years)` = NULL
pheno$CLIENT_SAMPLE_ID = pheno$`Sample Code/Client sample ID (NP11-XXX)`
pheno$`Sample Code/Client sample ID (NP11-XXX)` = NULL
pheno$Ethnicity = pheno$`Ethnicity (Arab, Asian, Caucasian)`
pheno$`Ethnicity (Arab, Asian, Caucasian)`= NULL
pheno$Gender = pheno$"Gender (male, female)"
pheno$"Gender (male, female)" = NULL
pheno$BMI = pheno$"BMI (kg*m-2)"
pheno$"BMI (kg*m-2)" = NULL
pheno_add = pheno[,c("CLIENT_SAMPLE_ID","Gender", "BMI", "Age","Ethnicity", "HbA1C")]
data_join <- merge(as.data.frame(colData(D_CLP)), pheno_add,by.x = "CLIENT_SAMPLE_ID", by.y = "CLIENT_SAMPLE_ID" ,all = TRUE)
data_join_select <- data_join %>% select(CLIENT_SAMPLE_ID, Gender, GENDER, BMI.x,BMI.y,Age,Ethnicity,RACE_ETHNICITY,HbA1C) %>%
mutate(gen = ifelse(Gender == GENDER, TRUE, FALSE)) %>%
mutate(bm = ifelse(BMI.x == BMI.y, TRUE, FALSE)) %>%
mutate(eth = ifelse(Ethnicity == RACE_ETHNICITY, TRUE, FALSE))
data_join_select %>% filter(gen== FALSE| bm == FALSE | eth == FALSE |(is.na(Gender) & !is.na(GENDER)))
#colData(D_CLP)$AGE = NULL
colData(D_CLP)$AGE = data_join$Age
colData(D_CLP)$HbA1C = as.numeric(data_join$HbA1C)
#colData(D_CLP)$GENDER = NULL
colData(D_CLP)$GENDER = data_join$Gender
colData(D_CLP)$RACE_ETHNICITY = NULL
colData(D_CLP)$ETHNICITY = data_join$Ethnicity
#colData(D_CLP)$BMI = NULL
colData(D_CLP)$BMI = data_join$BMI.y
colData(D_CLP)$GROUP_ID = colData(D_CLP)$Group
#############################################################################
# PART 2 - DATA CLEANING ------------plot and decide on sample and feature missingness % to filter later
##############################################################################
source("mt_plots_volcano_log_pvalue.R")
D <- D_CLP
D %>% mt_reporting_heading(heading = "Data Clean-up", lvl = 1) %>%
# filter samples
mt_modify_filter_samples(filter = !is.na(GROUP_ID)) %>%
# modify variable to factor
mt_anno_apply(anno_type = "samples", col_name = "GROUP_ID", fun = as.factor) %>%
mt_reporting_data()
# PART 3.1 - PREPROCESSING: ANALYSIS MISSING VALUES ----------------------------------------------------
colData(D)$Diabetic_status <- case_when(
(colData(D)$GROUP_ID == "Healthy" | colData(D)$GROUP_ID == "Psoriasis_Naïve" |colData(D)$GROUP_ID == "Psoriasis_FollowUp") ~ "Non-diabetic",
(colData(D)$GROUP_ID == "Psoriasis_Prediab_Naïve" | colData(D)$GROUP_ID == "Psoriasis_Prediab_FollowUp") ~ "Prediab",
TRUE ~ "diabetic"
)
save(D, file="SE.Rdata")