.Rhistory

)
any(duplicated((ranking_df$syndromic$RankAllAvg))
)
sum(duplicated((ranking_df$syndromic$RankAllAvg)))
head((ranking_df$syndromic$RankAllAvg[order((ranking_df$syndromic$RankAllAvg)], 50)
head(ranking_df$syndromic$RankAllAvg[order((ranking_df$syndromic$RankAllAvg)], 50)
head(ranking_df$syndromic$RankAllAvg[order(ranking_df$syndromic$RankAllAvg)], 50)
rank_select = ranking_df$syndromic %>% group_by(RankAllAvg) %>% arrange()
head(rank_select)
rank_select = ranking_df$syndromic %>% arrange(RankAllAvg)
head(rank_select)
rank_select = ranking_df$syndromic %>% arrange(RankAllAvg) %>% mutate(anyknown = trueset|candidate)
head(rank_select)
rank_vals = rank_select$RankAllAvg
rank_labs = rank_selectc$trueset
rank_labs = rank_select$trueset
rank_select = ranking_df$syndromic %>% arrange(RankAllAvg) %>% mutate(anyknown = trueset|candidate, cumsum = cumsum(trueset))
head(rank_select)
ggplot(rank_select, aes(x = RankAllAvg, y = cumsum))
ggplot(rank_select, aes(x = RankAllAvg, y = cumsum)) + geom_point()
ggplot(rank_select, aes(x = RankAllAvg, y = cumsum)) + geom_point() + scale_x_log10()
ggplot(rank_select, aes(x = RankAllAvg, y = cumsum)) + geom_point() + scale_x_log10() + scale_y_log10()
rank_select = ranking_df$syndromic %>% arrange(RankAllAvg) %>% mutate(anyknown = trueset|candidate, cumsum = cumsum(anyknown))
ggplot(rank_select, aes(x = RankAllAvg, y = cumsum)) + geom_point() + scale_x_log10() + scale_y_log10()
ggplot(rank_select, aes(x = RankAllAvg, y = cumsum)) + geom_point() #+ scale_x_log10() + scale_y_log10()
# TPR rate against diffferent ranking system
refset = "anyknown"
rank_select = ranking_df$syndromic %>% arrange(RankAllAvg) %>% mutate(anyknown = trueset|candidate, cumsum = cumsum(!!refset), TPR = cumsum/RankAllAvg)
head(rank_select)
rank_select = ranking_df$syndromic %>% arrange(RankAllAvg) %>% mutate(anyknown = trueset|candidate, cumsum = cumsum(!!!refset), TPR = cumsum/RankAllAvg)
head(rank_select)
rank_select = ranking_df$syndromic %>% arrange(RankAllAvg) %>% mutate(anyknown = trueset|candidate, cumsum = cumsum(refset), TPR = cumsum/RankAllAvg)
rank_select = ranking_df$syndromic %>% arrange(RankAllAvg) %>% mutate(anyknown = trueset|candidate, cumsum = cumsum(get(refset)), TPR = cumsum/RankAllAvg)
head(rank_select)
ggplot(rank_select, aes(x = RankAllAvg, y = TPR)) + geom_point() #+ scale_x_log10() + scale_y_log10()
ggplot(rank_select, aes(x = RankAllAvg, y = TPR)) + geom_point() + scale_x_log10() #+ scale_y_log10()
rank_criteria = "RankAllAvg"
rank_criteria = "RankAllAvg"
rank_select = ranking_df$syndromic %>% arrange(get(rank_criteria)) %>%
mutate(anyknown = trueset|candidate,  # count both known genes and candidate genes
cumsum = cumsum(get(refset)), # cumulative sum of true ref sets
TPR = cumsum/get(rank_criteria)) #true positive based on ranking system used
head(rank_select)
ggplot(rank_select, aes(x = get(rank_criteria), y = cumsum)) + geom_point() #+ scale_x_log10() + scale_y_log10()
ggplot(rank_select, aes(x = get(input$rank_criteria), y = cumsum)) + geom_point() #+ scale_x_log10() + scale_y_log10()
ggplot(rank_select, aes(x = get(rank_criteria), y = cumsum)) + geom_point() #+ scale_x_log10() + scale_y_log10()
ggplot(rank_select, aes(x = get(rank_criteria), y = TPR)) + geom_point()
ggplot(rank_select, aes(x = get(input$rank_criteria), y = TPR)) + geom_point() + scale_x_log10()
ggplot(rank_select, aes(x = get(irank_criteria), y = TPR)) + geom_point() + scale_x_log10()
ggplot(rank_select, aes(x = get(rank_criteria), y = TPR)) + geom_point() + scale_x_log10()
ggplot(rank_select, aes(x = get(rank_criteria), y = cumsum)) + geom_point() #+ scale_x_log10() + scale_y_log10()
ggplot(rank_select, aes(x = get(rank_criteria), y = cumsum)) + geom_point() + scale_x_log10() + scale_y_log10()
ggplot(rank_select, aes(x = get(rank_criteria), y = cumsum)) + geom_point(aes(fill = get(rank_criteria)==200)) + scale_x_log10() + scale_y_log10()
ggplot(rank_select, aes(x = get(rank_criteria), y = cumsum)) + geom_point(aes(col = get(rank_criteria)==200)) + scale_x_log10() + scale_y_log10()
ggplot(rank_select, aes(x = get(rank_criteria), y = TPR)) + geom_point() + scale_x_log10()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
a = ranking_df$syndromic %>% dplyr::filter(RankAllAvg <= 100)
head(a)
a %>% group_by(trueset, candidate) %>% summarise(paste(GeneName, collapse = "; "))
a %>% group_by(trueset, candidate) %>% summarise(n = length(GeneName), genes = paste(GeneName, collapse = "; "))
a %>% group_by(trueset, candidate) %>% summarise(n = length(GeneName), genes = paste(GeneName, collapse = "; "), allgenes = GeneName)
a %>% group_by(trueset, candidate) %>% summarise(n = length(GeneName), genes = paste(GeneName, collapse = "; "), allgenes = c(GeneName))
a %>% group_by(trueset, candidate) %>% summarise(n = length(GeneName), genes = paste(GeneName, collapse = "; "), allgenes = list(GeneName))
runApp()
runApp()
runApp()
library(reticulate)
a %>% group_by(trueset, candidate) %>% summarise(n = length(GeneName), genes = paste(GeneName, collapse = "; "), allgenes = list(GeneName))
a %>% group_by(trueset, candidate) %>% summarise(n = length(GeneName), genes = paste(GeneName, collapse = "; "), allgenes = list(GeneName)) -> b
b %>% dplyr::filter(!trueset & !candidate) %>% pull(allgenes)
b %>% dplyr::filter(!trueset & !candidate) %>% pull(allgenes) %>% unlist %>% paste(., collapse = "\n")
b %>% dplyr::filter(!trueset & !candidate) %>% pull(allgenes) %>% unlist %>% paste(., collapse = "\n") -> genes_str
genes_str
json = import('json')
requests = import('requests')
requests = import('requests')
requests = import('requests.py')
os = import("os")
os$listdir("")
os$listdir()
os$listdir(".")
py_available()
py_config()
requests = import('requests')
library(reticulate)
requests = import('requests')
dict(list = c(NULL, genes_str),
description = c(NULL, description))
description = 'Tomi top genes'
ENRICHR_URL = 'http://amp.pharm.mssm.edu/Enrichr/addList'
payload = dict(list = c(NULL, genes_str),
description = c(NULL, description))
payload
dict(list = c(None, genes_str),
description = c(None, description))
tuple(NULL, 'a')
payload = dict(list = tuple(NULL, genes_str),
description = tuple(NULL, description))
payload
reticulate::use_python(python = '/Library/Frameworks/Python.framework/Versions/3.7/bin/python3', required = T)
reticulate::use_python(python = '/Users/pbuphamalai/anaconda3/bin/python3', required = T)
reticulate::use_python(python = '/Users/pbuphamalai/anaconda3/bin/python3', required = T)
json = import('json')
library(reticulate)
library(pacman)
library(DT)
library(tidyverse)
library(ROCit)
library(plotly)
library(reticulate)
sys <- import("sys")
sys$version
requests = import('requests')
json = import('json')
requests = import('requests')
ENRICHR_URL = 'http://amp.pharm.mssm.edu/Enrichr/addList'
payload
payload = dict(list = tuple(NULL, genes_str),
description = tuple(NULL, description))
payload
response = requests$post(ENRICHR_URL, files = payload)
response
data = json$load(response$text)
response$text
data = json$loads(response$text)
data
data$shortId
a("Google Homepage", href="https://www.google.com/")
a(h4("Open Link"), target = "_blank", href = paste0("http://www.somesite/", some_link))
uiOutput(a(h4("Open Link"), target = "_blank", href = paste0("http://www.somesite/", some_link)))
library(shiny)
uiOutput(a(h4("Open Link"), target = "_blank", href = paste0("http://www.somesite/", some_link)))
uiOutput(a(h4("Open Link"), target = "_blank", href = paste0("http://www.somesite/", "abc")))
uiOutput(shiny::a(h4("Open Link", class = "btn btn-default action-button" ,
style = "fontweight:600"), target = "_blank",
href = paste0("http://www.somesite/", "abc")))
# use reticulate to retrieve API for enrichR
output$enrichr_url = renderUI(
{
json = import('json')
requests = import('requests')
ENRICHR_URL = 'http://amp.pharm.mssm.edu/Enrichr/addList'
# new genes to be discovered
genes_str = filtered_rank()  %>% dplyr::filter(!trueset & !candidate) %>% pull(allgenes) %>% unlist %>% paste(., collapse = "\n")
description = 'Tomi top genes'
payload = dict(list = tuple(NULL, genes_str),
description = tuple(NULL, description))
response = requests$post(ENRICHR_URL, files = payload)
data = json$loads(response$text)
url = paste0("http://amp.pharm.mssm.edu/Enrichr/enrich?dataset=",data$shortId)
tagList("URL link:", url)
}
)
renderUI(
{
json = import('json')
requests = import('requests')
ENRICHR_URL = 'http://amp.pharm.mssm.edu/Enrichr/addList'
# new genes to be discovered
genes_str = "ABCB1"#filtered_rank()  %>% dplyr::filter(!trueset & !candidate) %>% pull(allgenes) %>% unlist %>% paste(., collapse = "\n")
description = 'Tomi top genes'
payload = dict(list = tuple(NULL, genes_str),
description = tuple(NULL, description))
response = requests$post(ENRICHR_URL, files = payload)
data = json$loads(response$text)
url = paste0("http://amp.pharm.mssm.edu/Enrichr/enrich?dataset=",data$shortId)
tagList("URL link:", url)
}
)
runApp()
runApp()
runApp()
install.packages(shinyjs)
install.packages("shinyjs")
runApp()
runApp()
runApp()
?eventReactive()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
shiny::runApp()
librart
library("reticulate")
py_config()
shiny::runApp()
runApp()
runApp()
Sys.getenv("PATH")
runApp()
runApp()
pacman::p_load("vcfR")
my_vcf = read.vcfR("~/Documents/projects/myDNA/data/0820188482964.filtered.snp.vcf.gz", verbose = FALSE )
my_vcf = read.vcfR("~/Documents/projects/myDNA/data/60820188482964.filtered.snp.vcf.gz", verbose = FALSE )
head(my_vcf)
chrom <- create.chromR(name='Supercontig', vcf=my_vcf)
vcf_df = as_tibble(head(my_vcf))
vcf_df
my_vcf[1]
names(my_vcf)
length(my_vcf)
my_vcf[1:5]
my_vcf[1:5] %>% tibble()
my_vcf[1:5] %>% as.data.frame()
head(my_vcf@fix)
head(my_vcf@gt)
vcf_df = my_vcf@gt[,2]
strsplit("0/1:21,4:0.16:25:8,4:13,0:16:50,0,47:15.513,0.12382,49.774:0,34.77,37.77:9,12,3,1:8,13,3,1", ":")
strsplit("0/1:21,4:0.16:25:8,4:13,0:16:50,0,47:15.513,0.12382,49.774:0,34.77,37.77:9,12,3,1:8,13,3,1", ":")[[1]][4]
as.integer(strsplit("0/1:21,4:0.16:25:8,4:13,0:16:50,0,47:15.513,0.12382,49.774:0,34.77,37.77:9,12,3,1:8,13,3,1", ":")[[1]][4])
depth_extract = function(char){
as.integer(strsplit(char, ":"))[[1]][4]
}
depth_extract("0/1:21,4:0.16:25:8,4:13,0:16:50,0,47:15.513,0.12382,49.774:0,34.77,37.77:9,12,3,1:8,13,3,1" )
depth_extract = function(char){
as.integer(strsplit(char, ":")[[1]][4])
}
depth_extract("0/1:21,4:0.16:25:8,4:13,0:16:50,0,47:15.513,0.12382,49.774:0,34.77,37.77:9,12,3,1:8,13,3,1" )
depth_variants = sapply(vvcf_df, depth_extract)
depth_variants = sapply(vcf_df, depth_extract)
warnings()
length(depth_variants)
hist(depth_variants, breaks = 100)
length_df = as_tibble(length = depth_variants)
length_df = tibble(length = depth_variants)
length_df = length_df %>% count(length)
View(length_df)
ggplot(length_df) + geom_col(aes(x = length, y = n))
ggplot(length_df %>% dplyr::filter(!is.na(length))) + geom_col(aes(x = length, y = n)) + scale_x_log10()
View(length_df)
summary(depth_variants)
View(length_df)
##################
# Download ranking data
#################
ranking_df = list()
ranking_df$all = readRDS("~/Documents/projects/ID/IDome/cache/rank_all_ID_df.RDS")
write_tsv(ranking_df$all,  "~/Documents/projects/ID/IDome/cache/rank_all_ID_df.tsv")
ranking_df$severe = readRDS("~/Documents/projects/ID/IDome/cache/rank_severity_df.RDS")
ranking_df$syndromic = readRDS("~/Documents/projects/ID/IDome/cache/rank_syndromicity_df.RDS")
# ID candidate
sysid_genes = read_csv("~/Documents/projects/ID//data/SysID/sysID_list_20190513.csv")
candidate_genes = sysid_genes %>% dplyr::filter(grepl(pattern =  "candidate", x = `Gene group`)) %>% pull(`Gene symbol`) %>% unique
head(ranking_df$syndromic)
disease_df = read_csv("~/Documents/projects/ID/data/SysID/disease-info.csv")
tf_df = read_csv("~/Documents/projects/ID/data/SysID/transcription-factor.csv")
neuroscreen_df = read_csv("~/Documents/projects/ID/data/SysID/neuronal-screen.csv")
wingscreen_df = read_csv("~/Documents/projects/ID/data/SysID/wing-screen.csv")
sysid_genelist = read_csv("~/Documents/projects/ID/data/SysID/human-gene-info-20190619.csv")
# annotate information from the human protein atlas
housekeeping_df = read_tsv("~/Documents/projects/data/HPA/housekeeping_proteome.tsv")
tissue_rna_df = read_tsv("~/Documents/projects/data/HPA/rna_tissue_consensus.tsv")
brain_rna_df = read_tsv("~/Documents/projects/data/HPA/rna_brain_gtex.tsv")
mouse_allen_rna_df = read_tsv("~/Documents/projects/data/HPA/rna_mouse_brain_allen.tsv")
#  expression in combined HPA, GTEx
rna_expression_normalised = "./"
library(readxl)
PSD_exp = read_xlsx("./data/Brain_expression/barreslab_rnaseq.xlsx")
PSD_exp = read_xlsx("~/Documents/projects/ID/data/Brain_expression/barreslab_rnaseq.xlsx")
##############
# modify the expression data set: taking into account those that we are interested
columnnames = c("name","region","expression","data")
brain_rna_df_mod = brain_rna_df %>% dplyr::select(`Gene name`, `Brain region`, NX) %>% mutate(data = "Human brain")
tissue_rna_df_mod = tissue_rna_df %>% dplyr::select(`Gene name`, Tissue, NX) %>% filter(Tissue %in% c("adrenal gland", "bone marrow")) %>% mutate(data = "GTEx")
colnames(brain_rna_df_mod) = columnnames
colnames(tissue_rna_df_mod) = columnnames
expression_df = rbind(brain_rna_df_mod, tissue_rna_df_mod)
expression_df$ID = expression_df$name %in% ID_all$ID
View(expression_df)
apply(expression_df[,-1], 2, unique)
head(mouse_allen_rna_df)
head(disease_df)
colnames(tf_df)
colnames(neuroscreen_df)
# combine all the information
##########
additional_data_df = full_join(tf_df, neuroscreen_df %>% select("Human gene symbol", "Entrez id", "Any phenotype")
, by = c( "Gene symbol", "Human gene symbol"))
# combine all the information
##########
additional_data_df = full_join(tf_df, neuroscreen_df %>% select("Human gene symbol", "Entrez id", "Any phenotype")
, by = c( "Gene symbol"= "Human gene symbol"))
head(wingscreen_df)
colnames(wingscreen_df)
# combine all the information
##########
additional_data_df = full_join(tf_df, neuroscreen_df %>% select("Human gene symbol", "Entrez id", "Any phenotype")
, by = c( "Gene symbol"= "Human gene symbol")) %>%
rename(WingPhenotypes = `Any phenotype`)
head(additional_data_df)
# combine all the information
##########
additional_data_df = full_join(tf_df, neuroscreen_df %>% select("Human gene symbol",  "Any phenotype")
, by = c( "Gene symbol"= "Human gene symbol")) %>%
rename(WingPhenotypes = `Any phenotype`)
head(additional_data_df)
# combine all the information
##########
additional_data_df = full_join(tf_df, neuroscreen_df %>% select("Human gene symbol",  "Any phenotype")
, by = c( "Gene symbol"= "Human gene symbol")) %>%
rename(Neuroscreen_Phenotypes = `Any phenotype`)
additional_data_df = full_join(additional_data_df, wingscreen_df,
by = c( "Gene symbol"= "Any phenotype")) %>%
rename(WingScreen_Phenotypes = `Any phenotype`)
# combine all the information
##########
additional_data_df = full_join(tf_df, neuroscreen_df %>% select("Human gene symbol",  "Any phenotype")
, by = c( "Gene symbol"= "Human gene symbol")) %>%
rename(Neuroscreen_Phenotypes = `Any phenotype`)
additional_data_df = full_join(additional_data_df, wingscreen_df,
by = c( "Gene symbol"= "Human gene symbol")) %>%
rename(WingScreen_Phenotypes = `Any phenotype`)
head(additional_data_df)
# combine all the information
##########
additional_data_df = full_join(tf_df, neuroscreen_df %>% select("Human gene symbol",  "Any phenotype")
, by = c( "Gene symbol"= "Human gene symbol")) %>%
rename(Neuroscreen_Phenotypes = `Any phenotype`)
additional_data_df = full_join(additional_data_df, wingscreen_df %>% select("Human gene symbol",  "Any phenotype"),
by = c( "Gene symbol"= "Human gene symbol")) %>%
rename(WingScreen_Phenotypes = `Any phenotype`)
head(additional_data_df)
head(housekeeping_df)
colnames(housekeeping_df)
nrow(housekeeping_df)
View(housekeeping_df)
additional_data_df = full_join(additional_data_df, housekeeping_df %>% select(Gene) %>% mutate(HouseKeeping = T),
by = c( "Gene symbol"= "Gene"))
View(additional_data_df)
head(tf_df)
head(wingscreen_df)
View(neuroscreen_df)
View(wingscreen_df)
# modify the wing and neuroscreen: whether percent of evidence found to show any phenotypes
wingscreen_df_combn = wingscreen_df %>% group_by(`Human gene symbol`) %>%
summarise(WingPhenotype = mean(`Any phenotype`))
View(wingscreen_df_combn)
Neurocreen_df_combn = neuroscreen_df %>% group_by(`Human gene symbol`) %>%
summarise(NeuroScreenPhenotype = mean(`Any phenotype`))
View(Neurocreen_df_combn)
# modify the wing and neuroscreen: whether percent of evidence found to show any phenotypes
wingscreen_df_combn = wingscreen_df %>% group_by(`Human gene symbol`) %>%
summarise(WingPhenotype = mean(`Any phenotype`)>0)
Neurocreen_df_combn = neuroscreen_df %>% group_by(`Human gene symbol`) %>%
summarise(NeuroScreenPhenotype = mean(`Any phenotype`)>0)
View(Neurocreen_df_combn)
# combine all the information
##########
additional_data_df = full_join(tf_df, Neurocreen_df_combn
, by = c( "Gene symbol"= "Human gene symbol"))
additional_data_df = full_join(additional_data_df, wingscreen_df_combn,
by = c( "Gene symbol"= "Human gene symbol"))
additional_data_df = full_join(additional_data_df, housekeeping_df %>% select(Gene) %>% mutate(HouseKeeping = T),
by = c( "Gene symbol"= "Gene"))
head(expression_df)
# from long to wide
expression_df_wide = spread(expression_df, key = region, value = expression)
?spread
# from long to wide
expression_df_wide = spread(expression_df, key = name, value = expression)
unique(expression_df$region)
expression_df[c(170333, 208651),]
nrow(brain_rna_df_mod)
expression_df = expression_df %>% group_by(name, region) %>% distinct()
ncol(expression_df)
?distinct
expression_df = rbind(brain_rna_df_mod, tissue_rna_df_mod)
expression_df = distinct(expression_dfname, region)
expression_df = distinct(expression_df, name, region)
nrow(expression_df)
# from long to wide
expression_df_wide = spread(expression_df, key = region, value = expression)
# from long to wide
expression_df_wide = spread(expression_df, key = `region`, value = `expression`)
last_error()
rlang::last_error()
head(expression_df)
expression_df = rbind(brain_rna_df_mod, tissue_rna_df_mod)
expression_df = distinct(expression_df, name, region, .keep_all = T)
head(expression_df)
# from long to wide
expression_df_wide = spread(expression_df, key = `region`, value = `expression`)
View(expression_df_wide)
# from long to wide
expression_df_wide = spread(expression_df %>% select(-data), key = `region`, value = `expression`)
View(expression_df_wide)
additional_data_df = full_join(additional_data_df, expression_df_wide,
by = c( "Gene symbol"= "name"))
View(additional_data_df)
library(mygene)
# map gene name to other IDs
gene_data = read_tsv("~/Documents/projects/data/gene_info_reduced.txt")
View(gene_data)
gene_data = read_tsv("~/Documents/projects/data/geneDatabase.txt")
View(gene_data)
gene_data = read_tsv("~/Documents/projects/data/geneDatabase.txt", col_types = 'cc------dc---')
additional_data_df = left_join(additional_data_df, gene_data, by = c('Gene symbol' = "Approved symbol"))
additional_data_df = left_join(additional_data_df, gene_data, by = c('Gene symbol' = "Approved Symbol"))
sum(is.na(additional_data_df$`HGNC ID`))
head(PSD_exp)
colnames(PSD_exp)
View(PSD_exp)
write_tsv(additional_data_df, "~/Documents/projects/ID/IDome/cache/combined_gene_phenotypes_expression.tsv")
?readRDS
write_tsv(additional_data_df, "~/Documents/projects/ID/IDome/cache/combined_gene_phenotypes_expression.tsv")
saveRDS(additional_data_df,  "~/Documents/projects/ID/IDome/cache/combined_gene_phenotypes_expression.tsv")
write_tsv(additional_data_df, "~/Documents/projects/ID/IDome/cache/combined_gene_phenotypes_expression.tsv")
saveRDS(additional_data_df,  "~/Documents/projects/ID/IDome/cache/combined_gene_phenotypes_expression.RDS")
##################
# add gene info
phenodata = readRDS("../IDome/cache/combined_gene_phenotypes_expression.RDS")
head(ranking_df$all)
head(phenodata)
for(i in names(ranking_df)){
ranking_df[[i]] = left_join(ranking_df[[i]], phenodata, by = c("GeneName" = `Gene symbol`))
}
for(i in names(ranking_df)){
ranking_df[[i]] = left_join(ranking_df[[i]], phenodata, by = c("GeneName" = `Gene symbol`))
}
for(i in names(ranking_df)){
ranking_df[[i]] = left_join(ranking_df[[i]], phenodata, by = c("GeneName" = "Gene symbol"))
}
rm(list = ls())
##################
# Download ranking data
#################
ranking_df = list()
ranking_df$all = readRDS("~/Documents/projects/ID/IDome/cache/rank_all_ID_df.RDS")
write_tsv(ranking_df$all,  "~/Documents/projects/ID/IDome/cache/rank_all_ID_df.tsv")
ranking_df$severe = readRDS("~/Documents/projects/ID/IDome/cache/rank_severity_df.RDS")
ranking_df$syndromic = readRDS("~/Documents/projects/ID/IDome/cache/rank_syndromicity_df.RDS")
# ID candidate
sysid_genes = read_csv("~/Documents/projects/ID//data/SysID/sysID_list_20190513.csv")
candidate_genes = sysid_genes %>% dplyr::filter(grepl(pattern =  "candidate", x = `Gene group`)) %>% pull(`Gene symbol`) %>% unique
##################
# add gene info
phenodata = readRDS("../IDome/cache/combined_gene_phenotypes_expression.RDS")
for(i in names(ranking_df)){
ranking_df[[i]] = left_join(ranking_df[[i]], phenodata, by = c("GeneName" = "Gene symbol"))
}
##################
# Perform ROC computation
#################
roc = list()
roc$candidate = list()
roc$trueset = list()
rank_columns = c("RankAllAvg" , "RankEachAvg" ,"RankArithmP" ,"RankGeomP")
for(i in names(ranking_df)){
ranking_df[[i]]$candidate = ranking_df[[i]]$GeneName %in% candidate_genes
ranking_df[[i]] = ranking_df[[i]][order(ranking_df[[i]]$RankGeomP),]
roc$candidate[[i]] = apply(ranking_df[[i]][,rank_columns], 2, function(x) rocit(score=-x,class=ranking_df[[i]][,"candidate"]))
if(i!="all"){
roc$trueset[[i]] = apply(ranking_df[[i]][,rank_columns], 2, function(x) rocit(score=-x,class=ranking_df[[i]][,"trueset"]))
}
}
rm(list = ls())
load("shinydata.RData")
##################
# add gene info
phenodata = readRDS("../IDome/cache/combined_gene_phenotypes_expression.RDS")
for(i in names(ranking_df)){
ranking_df[[i]] = left_join(ranking_df[[i]], phenodata, by = c("GeneName" = "Gene symbol"))
}
rm(list = ls())
load('shinydata.RData')
##################
# Download ranking data
#################
ranking_df = list()
ranking_df$all = readRDS("~/Documents/projects/ID/IDome/cache/rank_all_ID_df.RDS")
write_tsv(ranking_df$all,  "~/Documents/projects/ID/IDome/cache/rank_all_ID_df.tsv")
ranking_df$severe = readRDS("~/Documents/projects/ID/IDome/cache/rank_severity_df.RDS")
ranking_df$syndromic = readRDS("~/Documents/projects/ID/IDome/cache/rank_syndromicity_df.RDS")
##################
# add gene info
phenodata = readRDS("../IDome/cache/combined_gene_phenotypes_expression.RDS")
save.image(file = "shinydata.RData")
head(ranking_df$all)
colnames(ranking_df$all)
##################
# add gene info
phenodata = readRDS("../IDome/cache/combined_gene_phenotypes_expression.RDS")
for(i in names(ranking_df)){
ranking_df[[i]] = left_join(ranking_df[[i]], phenodata, by = c("GeneName" = "Gene symbol"))
}
colnames(ranking_df$all)
save.image(file = "shinydata.RData")
##################
# Download ranking data
#################
ranking_df = list()
ranking_df$all = readRDS("~/Documents/projects/ID/IDome/cache/rank_all_ID_df.RDS")
#write_tsv(ranking_df$all,  "~/Documents/projects/ID/IDome/cache/rank_all_ID_df.tsv")
ranking_df$severe = readRDS("~/Documents/projects/ID/IDome/cache/rank_severity_df.RDS")
ranking_df$syndromic = readRDS("~/Documents/projects/ID/IDome/cache/rank_syndromicity_df.RDS")
names(ranking_df)
for(i in names(ranking_df)){
ranking_df[[i]] = left_join(ranking_df[[i]], phenodata, by = c("GeneName" = "Gene symbol"))
write_tsv(ranking_df[[i]], paste0("../IDome/cache/prediction_result_annotate_seed=",i,".tsv"))
}