{"cells":[{"metadata":{},"cell_type":"markdown","source":"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":{},"cell_type":"markdown","source":"# PCA for visualizing batch effects"},{"metadata":{},"cell_type":"markdown","source":"PCA is a linear dimensionality technique that is useful for visualizing and understanding the structure of high-dimensional data. It relies on the eigendecomposition of the covariance matrix of the data points. Here we plot the two principal components corresponding to the dimensions that explain the most variance of the data for different cell types by batches. We use 50 samples from each batch. The data pipeline developed for working with the data will be extended with more methods and functions in subsequent notebooks."},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np \nimport pandas as pd\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nfrom sklearn.decomposition import PCA\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\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_path='../input'\ncsv_path='../'","execution_count":null,"outputs":[]},{"metadata":{"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 ProcessData():\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_random_loader(self, nr_of_samples):\n        dataset= ImagesDS(data_path+'/train.csv', self.data_path)\n        sampler = SubsetRandomSampler(np.arange(nr_of_samplers, dtype=np.int64))\n        loader = torch.utils.data.DataLoader(dataset, batch_size=1, sampler=sampler, num_workers=2)\n        return loader\n    \n    def create_loader_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(csv_path+'train_pca.csv')\n        dataset = ImagesDS(csv_path+'train_pca.csv', self.data_path)\n        loader = D.DataLoader(dataset, batch_size=1, shuffle=False, num_workers=1)\n        return loader\n        \n    def flatten_data_for_PCA(self,cell_type,nr_of_samples):\n        loader=self.create_loader_for_pca(nr_of_samples,cell_type)\n        arr=torch.zeros((1572864,1))\n        for x, y in tqdm(loader):\n            #Flatten into 1D vector of features for PCA\n            x=x.flatten().view(1572864,1)\n            arr=torch.cat((arr,x),dim=1)\n        return arr[:,1:].numpy()\n    \n    def PCA_for_batches(self,cell_type,n_components,nr_of_samples):\n        '''\n        nr of samples denotes how many samples to take from each batch. \n        '''\n        arr=self.flatten_data_for_PCA(cell_type,nr_of_samples)\n        pca=PCA(n_components=n_components)\n        tr=pca.fit_transform(arr.T)\n        plt.plot(np.cumsum(pca.explained_variance_ratio_))\n        plt.title(cell_type)\n        plt.xlabel('number of components')\n        plt.ylabel('cumulative explained variance')\n        plt.show()\n        labels=[]\n        lab_dict={}\n        colors_=range(0,len(self.batch_names[cell_type]))\n        cmp = cm.get_cmap('viridis',len(self.batch_names[cell_type]))\n        for j in range(0,len(self.batch_names[cell_type])):\n            labels=labels+[j]*nr_of_samples\n            lab_dict[self.batch_names[cell_type][j]]=j\n        labels=np.array(labels)\n        for g in np.unique(self.batch_names[cell_type]):\n            ix=np.where(labels==lab_dict[g])\n            plt.scatter(tr[:,0].flatten()[ix], tr[:,1].flatten()[ix], cmap=cmp, label=g)\n        plt.legend()\n        plt.title(cell_type)\n        plt.xlabel('component 1')\n        plt.ylabel('component 2')\n        plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"proc=ProcessData(data_path,csv_path)\nproc.PCA_for_batches('HEPG2',100,50)\nproc.PCA_for_batches('HUVEC',100,50)\nproc.PCA_for_batches('RPE',100,50)\nproc.PCA_for_batches('U2OS',100,50)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.8"}},"nbformat":4,"nbformat_minor":1}