{"cells":[{"metadata":{},"cell_type":"markdown","source":"This algorithm uses the sklearn IncrementalPCA implementation in combination with pytorch data loader to study which principal components contain confounding batch effects. For a notebook on PCA see https://www.kaggle.com/mariakesa/pcaforvisualizingbatcheffects IncrementalPCA is more memory efficient than PCA and it could be used inside a training loop, though transforming the data should be done in a loop separate from training the model as we do here as I noticed that there were artefacts due to the incremental training of the model (e.g. the first batches looked different when transformed within the training loop). \n\nThis analysis can be performed on both the train and test set because PCA is unsupervised. Using this approach we can discard the PCA compoennts that contain confounding effects, reconstruct the images and train a conv net on the reconstructed image that does not contain the confounding information (coming up in a subsequent notebook). \n\n*Hey! I'm interested in joining a team. I'm currently working as a programmer at the Howard Hughes Medical Institute at Janelia, processing large-scale neural calcium imaging recordings https://www.janelia.org/lab/pachitariu-lab\nI have a master's degree in Applied Mathematics and I've worked as a data scientist in a team focused on big data. I've also done an internship in applying convolutional neural networks on satellite images. My github is here https://github.com/mariakesa If you would be interested in taking me into your team shoot me an email at maria.kesa@gmail.com"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nfrom sklearn.decomposition import IncrementalPCA\n\nfrom PIL import Image\n\nimport torch\nimport torch.nn as nn\nimport torch.utils.data as D\nimport torch.nn.functional as F\n\nimport torchvision\nfrom torchvision import transforms as T\n\nfrom tqdm import tqdm\n\nimport matplotlib.colors as colors\nfrom matplotlib import cm\nimport matplotlib.gridspec as gridspec\nimport gc\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_path='../input'\ncsv_path='../'","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"#ImagesDS was taken from https://www.kaggle.com/leighplt/densenet121-pytorch\nclass ImagesDS(D.Dataset):\n    def __init__(self, csv_file, img_dir, mode='train', site=1, channels=[1,2,3,4,5,6]):\n        \n        df = pd.read_csv(csv_file)\n        self.records = df.to_records(index=False)\n        self.channels = channels\n        self.site = site\n        self.mode = mode\n        self.img_dir = img_dir\n        self.len = df.shape[0]\n        \n    @staticmethod\n    def _load_img_as_tensor(file_name):\n        with Image.open(file_name) as img:\n            return T.ToTensor()(img)\n\n    def _get_img_path(self, index, channel):\n        experiment, well, plate = self.records[index].experiment, self.records[index].well, self.records[index].plate\n        return '/'.join([self.img_dir,self.mode,experiment,f'Plate{plate}',f'{well}_s{self.site}_w{channel}.png'])\n        \n    def __getitem__(self, index):\n        paths = [self._get_img_path(index, ch) for ch in self.channels]\n        img = torch.cat([self._load_img_as_tensor(img_path) for img_path in paths])\n        \n        if self.mode == 'train':\n            return img, self.records[index].sirna\n        else:\n            return img, self.records[index].id_code\n\n    def __len__(self):\n        \"\"\"\n        Total number of samples in the dataset\n        \"\"\"\n        return self.len\n    \nclass ProcessWithIncrementalPCA():\n    def __init__(self, data_path,csv_path):\n        self.data_path=data_path\n        self.csv_path=csv_path\n        self.batch_names={'HEPG2':['HEPG2-01','HEPG2-02','HEPG2-03','HEPG2-04','HEPG2-05','HEPG2-06','HEPG2-07'],\n                          'HUVEC':['HUVEC-01','HUVEC-02','HUVEC-03','HUVEC-04','HUVEC-05','HUVEC-06','HUVEC-07',\n                                  'HUVEC-08','HUVEC-09','HUVEC-10','HUVEC-11','HUVEC-12','HUVEC-13','HUVEC-14',\n                                  'HUVEC-15','HUVEC-16'],\n                          'RPE':['RPE-01','RPE-02','RPE-03','RPE-04','RPE-05','RPE-06','RPE-07'],\n                          'U2OS':['U2OS-01','U2OS-02','U2OS-03']                            \n        }\n\n    def create_csv_for_pca(self,nr_of_samples,cell_type):\n        df=pd.read_csv(self.data_path+'/train.csv')\n        if cell_type=='HEPG2':\n            batches=self.batch_names['HEPG2']\n        if cell_type=='HUVEC':\n            batches=self.batch_names['HUVEC']\n        if cell_type=='RPE':\n            batches=self.batch_names['RPE']\n        if cell_type=='U2OS':\n            batches=self.batch_names['U2OS']\n        df_lst=[]\n        for batch in batches:\n            ind=df['experiment']==batch\n            sub=np.array(df[ind].index)\n            sub_=list(df[ind].index)\n            nr_of_samples_in_batch=len(sub_)\n            generate_random_numbers=np.random.randint(0,nr_of_samples_in_batch,nr_of_samples)\n            sub=sub[generate_random_numbers]\n            df_lst=df_lst+list(sub)\n        df_=df.loc[df_lst,:]\n        df_.to_csv(self.csv_path+'train_pca.csv')\n\n    \n    def create_loader(self,nr_of_samples):\n        dataset = ImagesDS(self.csv_path+'train_pca.csv', self.data_path)\n        loader = D.DataLoader(dataset, batch_size=nr_of_samples, shuffle=False, num_workers=1)\n        return loader\n        \n    def incremental_PCA(self,cell_type,nr_of_samples,nr_components):\n        self.create_csv_for_pca(nr_of_samples,cell_type)\n        loader=self.create_loader(nr_of_samples)\n        pca_res=np.zeros((1,nr_components))\n        ipca = IncrementalPCA(n_components=nr_components)\n        for x, y in tqdm(loader):\n            #Flatten into 1D vector of features for PCA\n            x=x.flatten().view(nr_of_samples,1572864).numpy()\n            ipca.partial_fit(x)\n        loader=self.create_loader(nr_of_samples)\n        for x, y in tqdm(loader):\n            x=x.flatten().view(nr_of_samples,1572864).numpy()\n            tr=ipca.transform(x)\n            pca_res=np.vstack((pca_res,tr))\n        pca_array=np.array(pca_res)\n        components=ipca.components_\n        return pca_array[1:,:],components\n    \n    def plot_PCs(self,pca_array,nr_components):\n        data_to_plot=[]\n        for i in range(0,nr_components):\n            fig = plt.figure(1, figsize=(10, 5))\n            for j in range(0,int(pca_array.shape[0]/nr_components)):\n                # Create an axes instance\n                ax = fig.add_subplot(111)\n                # Create the boxplot\n                data_to_plot.append(pca_array[j*100:(j+1)*100,i])\n            bp = ax.boxplot(data_to_plot)\n            plt.title('PC'+str(i))\n            data_to_plot=[]\n            plt.show()\n            \n    def plot_remove_dimensions(self,pca_array,components,remove_dims):\n        for dim in remove_dims:\n            components=np.delete(components,dim,axis=0)\n            pca_array=np.delete(pca_array,dim,axis=1)\n        reconstruction=components.T@pca_array.T\n        reconstruction=reconstruction.reshape((6,512,512,pca_array.shape[0]))\n        gs = gridspec.GridSpec(6, 5)\n        fig = plt.figure(figsize=(24,20))\n        for sample in range(0,5):\n            for channel_dim in range(0,6):\n                ax=fig.add_subplot(gs[channel_dim, sample])\n                ax.imshow(reconstruction[channel_dim,:,:,sample].reshape(512,512))\n                ax.set_yticks([])\n                ax.set_xticks([])\n        plt.show()\n                \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"proc=ProcessWithIncrementalPCA(data_path,csv_path)\npca_array,components=proc.incremental_PCA('HEPG2',100,100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Plot the distribution of PC's according to batch. \nproc.plot_PCs(pca_array,100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Plot the full reconstructions\nproc.plot_remove_dimensions(pca_array,components,[])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Plot the reconstruction with the first PC removed\nproc.plot_remove_dimensions(pca_array,components,[0])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The reconstructed images look the same, so removing the first component from the reconstruction should not hinder a convnet, but this removes the largest batch effects that are contained in the first dimension. "}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}