{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# What is batch effect?\n> Batch effect refers to technical variation or non-biological differences between measurements of different groups of samples.\n\n> Although batch effect can be reduced by good experimental design, it is difficult to completely eradicate.\n\n\n> Batch effects can lead to inaccurate conclusions when their causes are correlated with one or more outcomes of interest in an experiment.\n\n- https://en.wikipedia.org/wiki/Batch_effect/\n- https://www.nature.com/articles/s41598-017-11110-6\n\n#### They are most commonly discussed in the context of genomics and high-throughput sequencing research, but they exist in other fields of science as well.\n","metadata":{}},{"cell_type":"markdown","source":"# So... can we identify and remove Batch Effects?\n\n- There's an excellent visual explanation from 10x Genomics.\n\n### Identify and correct Batch Effects. From 10x Genomics\n","metadata":{}},{"cell_type":"markdown","source":"<img src=\"https://cdn.10xgenomics.com/image/upload/v1629241143/analysis-guides/BatchCorrection-Intro.png\"  width=\"600\" height=\"300\">\n\n- As you can see from the image, there were clearly signs of batch effects from different experiment, and the right image shows the batch effect is diminished through processing. ","metadata":{}},{"cell_type":"markdown","source":"# How can one identify batch effect?\n- For the simplest, dimension reduction methods like PCA will do. \n- Let's figure out first 10 batches have batch effects or not!","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport os\nimport pyarrow.parquet as pq\nfrom tqdm import tqdm\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom sklearn.decomposition import PCA\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-02-10T04:25:33.076750Z","iopub.execute_input":"2023-02-10T04:25:33.077897Z","iopub.status.idle":"2023-02-10T04:25:34.087994Z","shell.execute_reply.started":"2023-02-10T04:25:33.077742Z","shell.execute_reply":"2023-02-10T04:25:34.086658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_batchdata(batch_num,batch_dir):\n    sensor_info_df = pq.read_table(batch_dir+'batch_'+str(batch_num)+'.parquet').to_pandas()\n    return(sensor_info_df)","metadata":{"execution":{"iopub.status.busy":"2023-02-10T04:25:34.090746Z","iopub.execute_input":"2023-02-10T04:25:34.091231Z","iopub.status.idle":"2023-02-10T04:25:34.097950Z","shell.execute_reply.started":"2023-02-10T04:25:34.091184Z","shell.execute_reply":"2023-02-10T04:25:34.096683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Getting first 10 batch data","metadata":{}},{"cell_type":"code","source":"batch_dir = '/kaggle/input/icecube-neutrinos-in-deep-ice/train/'\nbatchdat_lst=[]\nfor i in tqdm(range(10)):\n    i = i+1\n    batch_dat = get_batchdata(i,batch_dir)\n    batch_dat = batch_dat.sample(frac=0.01)\n    batch_dat['Batch'] = 'Batch'+str(i)\n    batchdat_lst.append(batch_dat)\nbatch_df = pd.concat(batchdat_lst)","metadata":{"execution":{"iopub.status.busy":"2023-02-10T04:25:46.058139Z","iopub.execute_input":"2023-02-10T04:25:46.058720Z","iopub.status.idle":"2023-02-10T04:26:21.784020Z","shell.execute_reply.started":"2023-02-10T04:25:46.058638Z","shell.execute_reply":"2023-02-10T04:26:21.779450Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plotting PCA, colored by each batch. (First 10 batches)\n- Two columns were selected : charge, time.","metadata":{}},{"cell_type":"code","source":"X = batch_df[['charge','time']]\npca = PCA(n_components=2)\nprintcipalComponents = pca.fit_transform(X)\npcdf = pd.DataFrame(data=printcipalComponents, columns = ['PC1', 'PC2'])\npcdf['Batch'] = batch_df['Batch'].tolist()\nsns.set(font_scale = 1.5,style = 'ticks')\nplt.figure(figsize = (10,5))\nax = sns.scatterplot(data = pcdf,x = 'PC1',y = 'PC2',hue = 'Batch',linewidth =0,s = 10)\nplt.legend(loc = [1.05,0])\nplt.title('Principal components, first 10 batches',weight = 'bold')","metadata":{"execution":{"iopub.status.busy":"2023-02-10T04:26:21.786782Z","iopub.execute_input":"2023-02-10T04:26:21.787174Z","iopub.status.idle":"2023-02-10T04:26:42.378841Z","shell.execute_reply.started":"2023-02-10T04:26:21.787140Z","shell.execute_reply":"2023-02-10T04:26:42.377814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Overall, total 10 batched seems to have similar characteristics.\n- But it seems that red dots (Batch 4) around PC1 ~35,000 is quite isolated from other batches.\n\n","metadata":{"execution":{"iopub.status.busy":"2023-02-09T06:37:53.006734Z","iopub.execute_input":"2023-02-09T06:37:53.007094Z","iopub.status.idle":"2023-02-09T06:37:53.015392Z","shell.execute_reply.started":"2023-02-09T06:37:53.007064Z","shell.execute_reply":"2023-02-09T06:37:53.014206Z"}}},{"cell_type":"markdown","source":"# How to remove batch effects ?\n- There are several renowned methods for removing batch effects in molecular biology ( Mainly for sequencing )\n- We're going to use pyComBat, one of the ComBat Series!\n- https://www.biorxiv.org/content/10.1101/2020.03.17.995431v2\n- https://github.com/epigenelabs/pyComBat","metadata":{}},{"cell_type":"code","source":"# !pip install combat\nfrom combat.pycombat import pycombat","metadata":{"execution":{"iopub.status.busy":"2023-02-10T04:29:09.592093Z","iopub.execute_input":"2023-02-10T04:29:09.592506Z","iopub.status.idle":"2023-02-10T04:29:09.599232Z","shell.execute_reply.started":"2023-02-10T04:29:09.592469Z","shell.execute_reply":"2023-02-10T04:29:09.597701Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Combat input pre-processing.\nbatch = batch_df['Batch'].tolist()\ncombat_input = batch_df.reset_index(drop = True)\ncombat_input = combat_input.drop(['Batch','sensor_id','auxiliary'],axis=1)\ncombat_input = combat_input.T","metadata":{"execution":{"iopub.status.busy":"2023-02-10T04:26:59.384878Z","iopub.execute_input":"2023-02-10T04:26:59.386947Z","iopub.status.idle":"2023-02-10T04:26:59.599019Z","shell.execute_reply.started":"2023-02-10T04:26:59.386890Z","shell.execute_reply":"2023-02-10T04:26:59.597472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Run pyComBat\n","metadata":{}},{"cell_type":"code","source":"df_corrected = pycombat(combat_input,batch)","metadata":{"execution":{"iopub.status.busy":"2023-02-10T04:26:59.603674Z","iopub.execute_input":"2023-02-10T04:26:59.604295Z","iopub.status.idle":"2023-02-10T04:27:34.976020Z","shell.execute_reply.started":"2023-02-10T04:26:59.604260Z","shell.execute_reply":"2023-02-10T04:27:34.974861Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"combat_corrected = df_corrected.T","metadata":{"execution":{"iopub.status.busy":"2023-02-10T04:27:34.977253Z","iopub.execute_input":"2023-02-10T04:27:34.977567Z","iopub.status.idle":"2023-02-10T04:27:35.354312Z","shell.execute_reply.started":"2023-02-10T04:27:34.977539Z","shell.execute_reply":"2023-02-10T04:27:35.352996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Check if the batch effect is removed?","metadata":{}},{"cell_type":"code","source":"X_corrected = combat_corrected\npca = PCA(n_components=2)\nprintcipalComponents_corrected = pca.fit_transform(X_corrected)\npcdf_corrected = pd.DataFrame(data=printcipalComponents_corrected, columns = ['PC1', 'PC2'])\npcdf_corrected['Batch'] = batch_df['Batch'].tolist()\nsns.set(font_scale = 1.5,style = 'ticks')\nplt.figure(figsize = (10,5))\nsns.scatterplot(data = pcdf_corrected,x = 'PC1',y = 'PC2',hue = 'Batch',linewidth =0,s = 10)\nplt.legend(loc = [1.05,0])\nplt.title('After removing batch effect')","metadata":{"execution":{"iopub.status.busy":"2023-02-10T04:27:35.356002Z","iopub.execute_input":"2023-02-10T04:27:35.356348Z","iopub.status.idle":"2023-02-10T04:27:55.884825Z","shell.execute_reply.started":"2023-02-10T04:27:35.356318Z","shell.execute_reply":"2023-02-10T04:27:55.883402Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Looks like batch effect related to the Red Dot ( Batch 4 ) is somewhat diminished!\n- Let's check the result in pair. ","metadata":{}},{"cell_type":"markdown","source":"## Result of batch effect removal","metadata":{}},{"cell_type":"code","source":"fig,axes = plt.subplots(1,2,figsize = (30,10),sharex = True,sharey = True)\n\nsns.scatterplot(data = pcdf,x = 'PC1',y = 'PC2',hue = 'Batch',linewidth =0,s = 10,ax = axes[0])\naxes[0].set_title('Before removing batch effect',weight = 'bold')\nsns.scatterplot(data = pcdf_corrected,x = 'PC1',y = 'PC2',hue = 'Batch',linewidth =0,s = 10,ax = axes[1])\naxes[1].set_title('After removing batch effect',weight = 'bold')\n","metadata":{"execution":{"iopub.status.busy":"2023-02-10T04:27:55.886597Z","iopub.execute_input":"2023-02-10T04:27:55.887352Z","iopub.status.idle":"2023-02-10T04:29:09.589327Z","shell.execute_reply.started":"2023-02-10T04:27:55.887298Z","shell.execute_reply":"2023-02-10T04:29:09.587736Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# There are over 600 batches in total, how can we compare all of them? ( On-going )\n- Plotting all batches in a single figure would cost us an immense amount of time and resource\n- And even if we plot it, the chances are high that we cannot distinguish other batches by dot size, color... etc by a single plot.\n- So why not try PCA explainable values?\n- Retrieve all pairs of 660 batches, and calculate PCA explained variance ratio and get the pair of batches with lowest ratio","metadata":{}},{"cell_type":"code","source":"batch_dir = '/kaggle/input/icecube-neutrinos-in-deep-ice/train/'\nbatchdat_lst=[]\nfor i in tqdm(range(660)):\n    i = i+1\n    batch_dat = get_batchdata(i,batch_dir)\n    batch_dat = batch_dat.sample(frac=0.0001)\n    batch_dat['Batch'] = 'Batch'+str(i)\n    batchdat_lst.append(batch_dat)\n    batch_df = pd.concat(batchdat_lst)\n\nX = batch_df[['charge','time']].to_numpy()\npca = PCA(n_components=2)\nprintcipalComponents = pca.fit_transform(X)\npca_x_ratio = pca.explained_variance_ratio_\n","metadata":{"execution":{"iopub.status.busy":"2023-02-10T04:25:39.900344Z","iopub.status.idle":"2023-02-10T04:25:39.903400Z","shell.execute_reply.started":"2023-02-10T04:25:39.903010Z","shell.execute_reply":"2023-02-10T04:25:39.903048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}