# %% [code] {"_execution_state":"idle","execution":{"iopub.status.busy":"2023-11-25T11:37:58.352872Z","iopub.execute_input":"2023-11-25T11:37:58.355209Z","iopub.status.idle":"2023-11-25T11:37:58.376760Z"}}
# This R environment comes with many helpful analytics packages installed
# It is defined by the kaggle/rstats Docker image: https://github.com/kaggle/docker-rstats
# For example, here's a helpful package to load

library(tidyverse) # metapackage of all tidyverse packages

# Input data files are available in the read-only "../input/" directory
# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory

list.files(path = "../input")

# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using "Save & Run All" 
# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session

# %% [markdown]
# ### Source
#  - https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4402510/

# %% [markdown]
# ### Overview
#  - Goal of this competition is to predict how small molecules change gene expression in different cell types.
#  - Develop methods to predict how cells respond to small molecule drug perturbations, which could have important applications in drug discovery and basic biology.
#  
# ### Submission File
#  - For each id in the evaluation set, you should predict a value for each of the 18,211 genes named in the remaining columns. Each id corresponds to a cell_type / sm_name pair, which you may identify from the id_map.csv file.
#  - submission should contain a header and have the following format:
#  ![image.png](attachment:14a55b0a-9837-43a2-8b8e-d7b5f049dd3e.png)
#  
#  
# ### Dataset Description
#   - Single-cell perturbational dataset in human peripheral blood mononuclear cells (PBMCs). Selected 144 compounds from the Library of Integrated Network-Based Cellular Signatures (LINCS) Connectivity Map dataset (PMID: 29195078) and measured single-cell gene expression profiles after 24 hours of treatment. The experiment was repeated in three healthy human donors, and the compounds were selected based on diverse transcriptional signatures observed in CD34+ hematopoietic stem cells (data not released). 
#   
# ### Technical details about the experiment
# - PBMCs from donors were thawed and plated on 96-well plates. Two columns of the plates were dedicated to positive controls (dabrfenib and belinostat) and one column was dedicated to a negative control (DMSO). The positive controls were selected because they tend to have a large impact on transcription, and the negative control is used as a solvent for the compounds used in this study. The remaining wells on the plate are allocated to each of 72 compounds. The full dataset comprises 2 different compound plates per donor for a total of 6 plates.![image.png](attachment:7bdab59e-7642-4314-b24d-e4b0935414b5.png)
# 
# ### calculate differential expression
# - Prticipants are tasked with modelling differential expression (DE), which enables us to estimate the impact of an experimental perturbation on the expression level of every gene in the transcription (18211 genes in this dataset). 
# - We estimate the impact of each compound by first averaging the raw gene expression counts in each cell of a specific type in each sample, which is called pseudobulking in the single-cell literature. We then fit a linear model to the pseudobulked counts data using Limma and include the library (row), plate, and donor as technical covariates and compound as the experimental covariate. Here, pseudobulked means we summed the raw counts for all cells of a given type for each well in the experiment.Diagram of this process is shown below:![image.png](attachment:6b12e723-0668-4b3a-8a24-f2be48f9301b.png)
# 
# ### Limma
# - limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. 
# 
# 
# ### Differential expression
#  - Differential gene expression, commonly abbreviated as DG or DGE analysis refers to the analysis and interpretation of differences in abundance of gene transcripts within a transcriptome . Lists of genes that differ between 2 sample sets are often provided by RNA-seq data analysis tools, or can be generated manually by statistical testing of data sets. Due to the large number of genes to be tested, (e.g., >20,000 in the human genome), multiple testing correction such as Bonferroni correction is usually applied.

# %% [code] {"execution":{"iopub.status.busy":"2023-11-25T11:38:54.673052Z","iopub.execute_input":"2023-11-25T11:38:54.675008Z","iopub.status.idle":"2023-11-25T11:45:02.902407Z"}}
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")

# BiocManager::install("SingleCellExperiment")


# %% [markdown]
# ### Example below
# - exprMatrix - Each column represents a sample, and each row represents a gene.
# - group variable represents the experimental groups (e.g., Control and Treatment). Modify this according to your experimental design.
# - design matrix is created to model the experimental design. Adjust matrix accordingly.
# - Linear model is fitted using the lmFit function.
# - Empirical Bayes moderated t-statistics are estimated using the eBayes function.
# - topTable function is then used to identify differentially expressed genes. The coef argument specifies the contrast of interest. Adjust this based on your experimental design.

# %% [code] {"execution":{"iopub.status.busy":"2023-11-25T11:46:23.864113Z","iopub.execute_input":"2023-11-25T11:46:23.870298Z","iopub.status.idle":"2023-11-25T11:46:23.973072Z"}}
library(limma)

# Generate example data
set.seed(123)
nGenes <- 1000
nSamples <- 10
exprMatrix <- matrix(rnorm(nGenes * nSamples), nrow = nGenes)
colnames(exprMatrix) <- paste("Sample", 1:nSamples, sep = "")
group <- factor(rep(c("Control", "Treatment"), each = nSamples/2))

# Create a design matrix
design <- model.matrix(~0 + group)

# Create a linear model
fit <- lmFit(exprMatrix, design)

# Estimate the empirical Bayes moderated t-statistics
fit <- eBayes(fit)

# Identify differentially expressed genes
topTable(fit, coef = 1)  # Change the 'coef' argument 

# %% [code] {"execution":{"iopub.status.busy":"2023-11-25T11:57:27.348705Z","iopub.execute_input":"2023-11-25T11:57:27.350607Z","iopub.status.idle":"2023-11-25T11:57:27.663700Z"}}
# library(limma)
# library(SingleCellExperiment)
# counts <- matrix(rpois(100, lambda = 10), ncol=10, nrow=10)
# sce <- SingleCellExperiment(counts)
# sce


# %% [code]