{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# CiteSeq Target Analysis (Blom function) \n\nI would like to perform an analysis from train_input data with the goal of  important feature selection for each variable (day, cell_type, donor ).\n\nIn this case, we attempted to convert the data to Normal score.","metadata":{}},{"cell_type":"code","source":"library(rhdf5)\nh5ls(\"../input/open-problems-multimodal/train_cite_targets.h5\")\nh5ls(\"../input/open-problems-multimodal/train_cite_inputs.h5\")","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:35:12.421330Z","iopub.execute_input":"2022-11-02T11:35:12.423654Z","iopub.status.idle":"2022-11-02T11:35:13.604731Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"protein <- h5read(\"../input/open-problems-multimodal/train_cite_targets.h5\", '/train_cite_targets/axis0')\nhead(protein)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:35:13.609482Z","iopub.execute_input":"2022-11-02T11:35:13.652005Z","iopub.status.idle":"2022-11-02T11:35:13.695969Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"id <- h5read(\"../input/open-problems-multimodal/train_cite_targets.h5\", '/train_cite_targets/axis1')\nhead(id)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:35:13.699519Z","iopub.execute_input":"2022-11-02T11:35:13.706302Z","iopub.status.idle":"2022-11-02T11:35:13.774588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"values <- h5read(\"../input/open-problems-multimodal/train_cite_targets.h5\", '/train_cite_targets/block0_values')\nhead(values)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:35:13.777700Z","iopub.execute_input":"2022-11-02T11:35:13.779189Z","iopub.status.idle":"2022-11-02T11:35:15.039047Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dat <- as.data.frame(values)\nrownames(dat) <- protein\ncolnames(dat) <- id\n\nhead(t(dat))\ndim(t(dat))","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:35:15.043032Z","iopub.execute_input":"2022-11-02T11:35:15.050489Z","iopub.status.idle":"2022-11-02T11:35:19.812419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"metadata <- read.csv('../input/open-problems-multimodal/metadata.csv',row.names=1)\nhead(metadata)\ndim(metadata)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:35:19.819431Z","iopub.execute_input":"2022-11-02T11:35:19.824857Z","iopub.status.idle":"2022-11-02T11:35:21.466887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df <- merge(t(dat) , metadata, by=0)\nhead(df)\ndim(df)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:35:21.474492Z","iopub.execute_input":"2022-11-02T11:35:21.480375Z","iopub.status.idle":"2022-11-02T11:35:25.023834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Df <- df[,2:141]\nrownames(Df) <- df[,1]\nDf$cell_type <- df$cell_type\n\nhead(Df)\ndim(Df)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:38:12.523982Z","iopub.execute_input":"2022-11-02T11:38:12.525655Z","iopub.status.idle":"2022-11-02T11:38:12.590304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"table(Df$cell_type)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:39:16.074278Z","iopub.execute_input":"2022-11-02T11:39:16.076034Z","iopub.status.idle":"2022-11-02T11:39:16.103076Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list <- seq(1,140,1)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T14:11:36.866730Z","iopub.execute_input":"2022-11-02T14:11:36.869581Z","iopub.status.idle":"2022-11-02T14:11:36.890431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"FUNK <- function(c){\n    p <- ks.test(Df[,c],\"pnorm\")\n    return(p[2])\n}","metadata":{"execution":{"iopub.status.busy":"2022-11-02T14:14:21.491205Z","iopub.execute_input":"2022-11-02T14:14:21.496268Z","iopub.status.idle":"2022-11-02T14:14:21.518268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"P_pre <- sapply(list,FUNK)\nsum(P_pre < 0.05)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T14:14:22.614886Z","iopub.execute_input":"2022-11-02T14:14:22.616509Z","iopub.status.idle":"2022-11-02T14:14:25.410814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**The expression levels of all 140 proteins are non-normally distributed.**","metadata":{}},{"cell_type":"markdown","source":"# Blom function (Normal scores transformation)","metadata":{}},{"cell_type":"markdown","source":"**I want to download \"rcompanion\" with the latest version of R, but...**","metadata":{}},{"cell_type":"code","source":"R.version","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:41:46.255560Z","iopub.execute_input":"2022-11-02T11:41:46.257705Z","iopub.status.idle":"2022-11-02T11:41:46.277412Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# install.packages(\"installr\")","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:44:18.916093Z","iopub.execute_input":"2022-11-02T11:44:18.918629Z","iopub.status.idle":"2022-11-02T11:44:34.174492Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# library(installr)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:44:41.293994Z","iopub.execute_input":"2022-11-02T11:44:41.295647Z","iopub.status.idle":"2022-11-02T11:44:41.382403Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# updateR(TRUE)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:54:43.663665Z","iopub.execute_input":"2022-11-02T11:54:43.666653Z","iopub.status.idle":"2022-11-02T11:54:45.889200Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# install.packages(\"rcompanion\")","metadata":{"execution":{"iopub.status.busy":"2022-11-02T11:44:56.602595Z","iopub.execute_input":"2022-11-02T11:44:56.604169Z","iopub.status.idle":"2022-11-02T11:44:57.483999Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**I would appreciate it if Kaggle Staff could bring it up to date, as new libraries come out one after another.**","metadata":{}},{"cell_type":"code","source":"# devtools::install_github(\"cran/rcompanion\", lib = \"/kaggle/working\")","metadata":{"execution":{"iopub.status.busy":"2022-11-02T12:07:05.399035Z","iopub.execute_input":"2022-11-02T12:07:05.400909Z","iopub.status.idle":"2022-11-02T12:09:54.224722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Installation from Github doesn't work either.**","metadata":{}},{"cell_type":"code","source":"blom = function(x, method=\"general\", alpha=pi/8, \n                complete=FALSE, na.last=\"keep\", na.rm=TRUE,\n                adjustN=TRUE,\n                min=1, max=10, ...){\n  if(complete){x=x[complete.cases(x)]}\n  Ranks = rank(x, na.last=na.last, ...)\n  if(adjustN==FALSE){N = length(x)}\n  if(adjustN==TRUE) {N = sum(complete.cases(x))}\n  if(method==\"blom\")   {Score = qnorm((Ranks-0.375)/(N+0.25))}\n  if(method==\"vdw\")    {Score = qnorm((Ranks)/(N+1))}\n  if(method==\"tukey\")  {Score = qnorm((Ranks-1/3)/(N+1/3))}\n  if(method==\"rankit\") {Score = qnorm((Ranks-1/2)/(N))}\n  if(method==\"elfving\"){Score = qnorm((Ranks-pi/8)/(N-pi/4+1))}\n  if(method==\"general\"){Score = qnorm((Ranks-alpha)/(N-2*alpha+1))}\n  if(method==\"zscore\") {Score = (x-mean(x, na.rm=na.rm))/sd(x, na.rm=na.rm)}\n  if(method==\"scale\")  {\n    Score = (((x - min(x, na.rm=na.rm)) * \n            (max - min)) / \n            (max(x, na.rm=na.rm) - min(x, na.rm=na.rm)) + \n             min)\n  }\n  return(Score)\n}","metadata":{"execution":{"iopub.status.busy":"2022-11-02T12:38:59.203342Z","iopub.execute_input":"2022-11-02T12:38:59.205238Z","iopub.status.idle":"2022-11-02T12:38:59.219911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**↑**　\n\n**Copy and paste the script.\nThis is a primitive method.\nIf I can use the blom function, it's good anyway.**","metadata":{}},{"cell_type":"code","source":"list <- seq(1,140,1)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:15:54.085078Z","iopub.execute_input":"2022-11-02T13:15:54.087315Z","iopub.status.idle":"2022-11-02T13:15:54.102867Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"func <- function(c){\n    res <- blom(Df[,c],method='blom')\n    return(res)\n}","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:15:54.819462Z","iopub.execute_input":"2022-11-02T13:15:54.821120Z","iopub.status.idle":"2022-11-02T13:15:54.834738Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DF <- sapply(list,func)\ndim(DF)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:16:20.256311Z","iopub.execute_input":"2022-11-02T13:16:20.258259Z","iopub.status.idle":"2022-11-02T13:16:23.317729Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DF <- as.data.frame(DF)\nrownames(DF) <- rownames(Df)\ncolnames(DF) <- colnames(Df[1:140])","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:48:29.229909Z","iopub.execute_input":"2022-11-02T13:48:29.234225Z","iopub.status.idle":"2022-11-02T13:48:29.280972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(DF)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:48:41.451291Z","iopub.execute_input":"2022-11-02T13:48:41.453145Z","iopub.status.idle":"2022-11-02T13:48:41.499908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Checking \"blom function (Normal scores transformation)\"**","metadata":{}},{"cell_type":"code","source":"Func <- function(c){\n    p <- ks.test(DF[,c],\"pnorm\")\n    return(p[2])\n}","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:24:46.511912Z","iopub.execute_input":"2022-11-02T13:24:46.515609Z","iopub.status.idle":"2022-11-02T13:24:46.537932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"P <- sapply(list,Func)\nsum(P < 0.05)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:32:34.998747Z","iopub.execute_input":"2022-11-02T13:32:35.000407Z","iopub.status.idle":"2022-11-02T13:32:37.578415Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**All protein expression levels are normally distributed, although a warning message is output in the presence of ties.**","metadata":{}},{"cell_type":"code","source":"hist(Df$CD86,breaks=100,freq=FALSE)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:46:34.372131Z","iopub.execute_input":"2022-11-02T13:46:34.374642Z","iopub.status.idle":"2022-11-02T13:46:34.462007Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**↑**\n\n**Before \"Normal scores transformation\" **","metadata":{}},{"cell_type":"code","source":"hist(DF$CD86,breaks=100,freq=FALSE)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:48:47.152556Z","iopub.execute_input":"2022-11-02T13:48:47.154226Z","iopub.status.idle":"2022-11-02T13:48:47.242175Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**↑**\n\n**After \"Normal scores transformation\" **","metadata":{}},{"cell_type":"code","source":"DF$day <- df$day\nDF$donor <- df$donor\nDF$cell_type <- df$cell_type","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:51:53.361391Z","iopub.execute_input":"2022-11-02T13:51:53.364281Z","iopub.status.idle":"2022-11-02T13:51:53.391234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dim(DF)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:52:00.211462Z","iopub.execute_input":"2022-11-02T13:52:00.214911Z","iopub.status.idle":"2022-11-02T13:52:00.294107Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"write.csv(DF,'CiteSeq_blom.csv',quote=F)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:52:39.497980Z","iopub.execute_input":"2022-11-02T13:52:39.500689Z","iopub.status.idle":"2022-11-02T13:52:53.093582Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**The next step is to proceed with the analysis with this normally distributed data**","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}