# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .sos
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: SoS
#     language: sos
#     name: sos
# ---

# + [markdown] kernel="SoS" tags=[]
# ## MWE
# MWE is uploaded in here:https://drive.google.com/drive/u/0/folders/1Ny3fzShONB-kfhvqcvnzL<KEY>J

# + kernel="SoS"
sos run pipeline/per_gene_data_merger.ipynb data_merger \
--genoFile ./mwe.region_plink_files/plink_files_list.txt \
--cwd MWE/rds_per_gene/ \
--region-list MWE.region.list \
--phenoFile MWE.phenotype.list \
--covFile MWE.covar.list &

# + kernel="SoS"
[global]
parameter: genoFile = paths
parameter: phenoFile = path
parameter: covFile = path
parameter: region_list = path
parameter: cwd = path
parameter: name = "demo"
import pandas as pd
region_tbl = pd.read_csv(region_list,sep = "\t")
genoFile = pd.read_csv(genoFile,sep = "\t",names = ["gene_id","path"],header = 0).merge(region_tbl,on = "gene_id").to_dict("records")

# + kernel="SoS"
[data_merger]
input: phenoFile,covFile, for_each = "genoFile"
output: f'{cwd:a}/{name}.{_genoFile["gene_id"]}.rds'
R:  expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
    library("plink2R")
    library("dplyr")
    library("readr")
    library("purrr")
    ## Define function
    read_gene_pheno = function(path){
      arg = paste0(c("tabix -h ",path," ${_genoFile["#chr"]}:${_genoFile["start"]}-${_genoFile["start"]+1}"),collapse = "")
      result = system(arg,intern = T)
      output = read.table(text = result[2], sep = "\t")
      colnames(output) = result[1]%>%stringr::str_split("\t")%>%unlist()
      return(output)
      }
    ## Data Loader
    gene_id = "${_genoFile["gene_id"]}"
    ### Genotype
    geno = read_plink("${path(_genoFile["path"]):n}")
    X = geno$bed
    ### Phenotype
    phenotype_list = read_delim("${_input[0]}","\t")
    covar_list = read_delim("${_input[1]}","\t")
    covar_list = covar_list%>%mutate(covar = map(path, ~read_delim(.x,"\t")%>%select(-`#id`)%>%na.omit%>%t()))
    phenotype_list = inner_join(phenotype_list,covar_list, by = "tissue")
    phenotype_list = phenotype_list%>%mutate(Y = map(path.x, ~read_gene_pheno(.x)%>%select(-c(`#chr`,start,end,gene_id))%>%t%>%as.matrix))%>%mutate(
                                            #### Get residue for each of tissue
                                              Y_resid = map2(Y,covar,~.lm.fit(x = .y, y = .x)$residuals%>%t%>%as_tibble))
    
    Y_resid = phenotype_list%>%select(Y_resid)%>%tidyr::unnest(Y_resid)%>%t%>%as.matrix
    colnames(Y_resid) = phenotype_list$tissue
    list(X,Y_resid)%>%saveRDS("${_output}")

StarcoderdataJupyterScriptsDedupFiltered