load R packages
library(dplyr)
library(tibble)
library(data.table)
library(impute)
load phenotype
phen <- read.csv("../../Study/Phenotype/clinical__information.csv")
prof <- fread("../../Study/RawData/samples_expression.txt")
curation
choose index of the phenotype for survival analysis
phen.cln <- phen %>%
mutate(stage_event_tnm_categories = paste0(pathologic_T, pathologic_N, pathologic_M)) %>%
dplyr::select(bcr_patient_barcode, gender, vital_status, days_to_death,
days_to_last_followup, race, person_neoplasm_cancer_status,
age_at_initial_pathologic_diagnosis, laterality,
neoplasm_histologic_grade, pathologic_stage, stage_event_tnm_categories) %>%
distinct(bcr_patient_barcode, .keep_all = TRUE)
# dead_patient
dead_patient <- phen.cln %>%
dplyr::filter(vital_status == "Dead") %>%
dplyr::select(-days_to_last_followup) %>%
reshape::rename(c(bcr_patient_barcode = "Barcode",
gender = "Gender",
vital_status = "OS",
days_to_death = "OS.Time",
race = "Race",
person_neoplasm_cancer_status = "cancer_status",
age_at_initial_pathologic_diagnosis = "Age",
neoplasm_histologic_grade = "Grade",
pathologic_stage = "Stage",
stage_event_tnm_categories = "TNM")) %>%
mutate(OS = ifelse(OS == "Dead", 1, 0)) %>%
mutate(OS.Time = as.numeric(as.character(OS.Time))/365)
# alive_patient
alive_patient <- phen.cln %>%
dplyr::filter(vital_status == "Alive") %>%
dplyr::select(-days_to_death) %>%
reshape::rename(c(bcr_patient_barcode = "Barcode",
gender = "Gender",
vital_status = "OS",
days_to_last_followup = "OS.Time",
race = "Race",
person_neoplasm_cancer_status = "cancer_status",
age_at_initial_pathologic_diagnosis = "Age",
neoplasm_histologic_grade = "Grade",
pathologic_stage = "Stage",
stage_event_tnm_categories = "TNM")) %>%
mutate(OS = ifelse(OS == "Dead", 1, 0)) %>%
mutate(OS.Time = as.numeric(as.character(OS.Time))/365)
# combine data
survival_data <- rbind(dead_patient, alive_patient) %>%
mutate(Barcode = gsub("-", "_", Barcode))
filter expression matrix
-
First, the biological features were removed if having zero value in more than 20% of patients.
-
The samples were removed if missing across more than 20% features.
-
Then we used impute function from R impute package, to fill out the missing values.
-
Lastly, we removed input features with zero values across all samples.
get_profile <- function(profile = prof,
metadata = survival_data,
occurrence = 0.2,
present = 0.2){
# profile = prof
# metadata = survival_data
# occurrence = 0.2
# present = 0.2
colnames(profile) <- c(colnames(profile)[1:3], gsub("-", "_", colnames(profile)[4:ncol(profile)]))
filter_fill_fun <- function(tag){
# tag = data_type[1]
dat <- profile %>% filter(Platform == tag) %>%
dplyr::select(-c("Platform", "Description"))
# feature occurrence
feature_occ <- apply(dat[, -1], 1, function(x){sum(x[!is.na(x)] != 0)/length(x)}) %>%
data.frame() %>%
setNames("feature_occurrence") %>%
mutate(GeneSymbol = dat$GeneSymbol)
remain_feature <- feature_occ %>% filter(feature_occurrence > occurrence)
# sample occurrence
sample_occ <- apply(dat[, -1], 2, function(x){sum(x[!is.na(x)] != 0)/length(x)}) %>%
data.frame() %>%
setNames("sample_occurrence") %>%
rownames_to_column("sample")
remain_sample <- sample_occ %>% filter(sample_occurrence > present)
# remove duplicate features
dat.cln <- dat %>% filter(GeneSymbol%in%as.character(remain_feature$GeneSymbol)) %>%
dplyr::select(c("GeneSymbol", as.character(remain_sample$sample)))
dat.cln$median <- apply(dat.cln[, -1], 1, median)
dat.cln <- with(dat.cln, dat.cln[order(GeneSymbol, median, decreasing = T), ])
dat.symbol <- dat.cln[!duplicated(dat.cln$GeneSymbol), ] %>%
dplyr::select(-median) %>%
column_to_rownames("GeneSymbol")
# impute.knn
set.seed(12345)
saved.state <- .Random.seed
dat.imputed <- impute.knn(as.matrix(dat.symbol))
return(dat.imputed)
}
data_type <- unique(profile$Platform)
#print(data_type)
copyNumber <- filter_fill_fun(data_type[1])
geneExp <- filter_fill_fun(data_type[2])
methylation <- filter_fill_fun(data_type[3])
protein_RPPA <- filter_fill_fun(data_type[4])
copyNumber_prf <- copyNumber$data
geneExp_prf <- geneExp$data
methylation_prf <- methylation$data
protein_RPPA_prf <- protein_RPPA$data
common_sampleid <- intersect(metadata$Barcode,
intersect(colnames(copyNumber_prf),
intersect(colnames(geneExp_prf),
intersect(colnames(methylation_prf), colnames(protein_RPPA_prf))
)
)
)
metadata.cln <- metadata %>% filter(Barcode %in% common_sampleid)
copyNumber_prf.cln <- copyNumber_prf[, colnames(copyNumber_prf)%in%common_sampleid]
geneExp_prf.cln <- geneExp_prf[, colnames(geneExp_prf)%in%common_sampleid]
methylation_prf.cln <- methylation_prf[, colnames(methylation_prf)%in%common_sampleid]
protein_RPPA_prf.cln <- protein_RPPA_prf[, colnames(protein_RPPA_prf)%in%common_sampleid]
res <- list(meta = metadata.cln,
CN = copyNumber_prf.cln,
GE = geneExp_prf.cln,
ML = methylation_prf.cln,
PR = protein_RPPA_prf.cln)
return(res)
}
data_trim_res <- get_profile(profile = prof,
metadata = survival_data)
output
dir_prf <- "../../Result/profile"
if(!dir.exists(dir_prf)){
dir.create(dir_prf)
}
# copyNumber
copyNumber.path <- paste(dir_prf, "copyNumber_filter.tsv", sep = "/")
write.table(data_trim_res$CN, file = copyNumber.path, quote = F, sep = "\t")
# geneExp
geneExp.path <- paste(dir_prf, "geneExp_filter.tsv", sep = "/")
write.table(data_trim_res$GE, file = geneExp.path, quote = F, sep = "\t")
# methylation
methylation.path <- paste(dir_prf, "methylation_filter.tsv", sep = "/")
write.table(data_trim_res$ML, file = methylation.path, quote = F, sep = "\t")
# protein_RPPA
protein_RPPA.path <- paste(dir_prf, "protein_RPPA_filter.tsv", sep = "/")
write.table(data_trim_res$PR, file = protein_RPPA.path, quote = F, sep = "\t")
dir_phe <- "../../Result/phenotype"
if(!dir.exists(dir_phe)){
dir.create(dir_phe)
}
# total
all_survival.path <- paste(dir_phe, "all_survival_data.tsv", sep = "/")
write.table(survival_data, file = all_survival.path, quote = F, sep = "\t", row.names = F)
# total
common_survival.path <- paste(dir_phe, "common_survival_data.tsv", sep = "/")
write.table(data_trim_res$meta, file = common_survival.path, quote = F, sep = "\t", row.names = F)
version
sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS/LAPACK: /data/share/anaconda3/lib/libopenblasp-r0.3.12.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] impute_1.60.0 ggpubr_0.4.0 ggplot2_3.3.3 data.table_1.13.6 tibble_3.0.4 dplyr_1.0.2
loaded via a namespace (and not attached):
[1] Rcpp_1.0.5 lattice_0.20-41 tidyr_1.1.2 rprojroot_2.0.2 assertthat_0.2.1 digest_0.6.27 R6_2.5.0
[8] cellranger_1.1.0 plyr_1.8.6 backports_1.2.1 evaluate_0.14 pillar_1.4.7 rlang_0.4.10 curl_4.3
[15] readxl_1.3.1 rstudioapi_0.13 car_3.0-10 vegan_2.5-7 Matrix_1.3-0 rmarkdown_2.6 desc_1.2.0
[22] splines_3.6.3 stringr_1.4.0 foreign_0.8-76 munsell_0.5.0 broom_0.7.3 compiler_3.6.3 xfun_0.19
[29] pkgconfig_2.0.3 mgcv_1.8-33 htmltools_0.5.0 tidyselect_1.1.0 rio_0.5.16 reshape_0.8.8 fansi_0.4.1
[36] permute_0.9-5 crayon_1.3.4 withr_2.3.0 MASS_7.3-53 grid_3.6.3 nlme_3.1-150 gtable_0.3.0
[43] lifecycle_0.2.0 magrittr_2.0.1 scales_1.1.1 zip_2.1.1 cli_2.2.0 stringi_1.5.3 carData_3.0-4
[50] ggsignif_0.6.0 testthat_3.0.1 ellipsis_0.3.1 generics_0.1.0 vctrs_0.3.6 cowplot_1.1.1 openxlsx_4.2.3
[57] tools_3.6.3 forcats_0.5.0 glue_1.4.2 purrr_0.3.4 hms_0.5.3 abind_1.4-5 pkgload_1.1.0
[64] parallel_3.6.3 yaml_2.2.1 colorspace_2.0-0 cluster_2.1.0 rstatix_0.6.0 knitr_1.30 haven_2.3.1
网友评论