Tutorial/Tutorial_Intra.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "f964c2ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "##Import the required packages\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "import os\n",
    "from sklearn.metrics import accuracy_score,f1_score,recall_score,precision_score\n",
    "import math\n",
    "import pandas as pd\n",
    "import scanpy as sc\n",
    "from tqdm.auto import tqdm\n",
    "from torch.utils.data import (DataLoader,Dataset)\n",
    "torch.set_default_tensor_type(torch.DoubleTensor)\n",
    "import numpy as np\n",
    "import random\n",
    "from sklearn.model_selection import KFold, StratifiedKFold\n",
    "from sklearn import preprocessing\n",
    "##Set random seeds\n",
    "def same_seeds(seed):\n",
    "    random.seed(seed)\n",
    "    # Numpy\n",
    "    np.random.seed(seed)\n",
    "    # Torch\n",
    "    torch.manual_seed(seed)\n",
    "    if torch.cuda.is_available():\n",
    "        torch.cuda.manual_seed(seed)\n",
    "        torch.cuda.manual_seed_all(seed)\n",
    "    torch.backends.cudnn.benchmark = False\n",
    "    torch.backends.cudnn.deterministic = True\n",
    "\n",
    "same_seeds(2021)\n",
    "\n",
    "##function\n",
    "    #getting the input of CIForm\n",
    "##parameters\n",
    "    #gap              :the length of a sub-vector \n",
    "    #data_path        :the paths of training datasets\n",
    "    #topgenes         :the number of highly variable genes(HVGs)\n",
    "    #label_path       :the path of the cell type annotation file(.csv) of training datasets\n",
    "    #T                 \n",
    "    #                  T==True：The rows represent the cells and the columns represent the genes\n",
    "    #                  T=False: The rows represent the genes and the columns represent the cells\n",
    "def getData(gap, data_path,topgenes,label_path=None,T=False):\n",
    "    \n",
    "    #getting the scRNA-seq datasets\n",
    "    adata = sc.read_csv(data_path, delimiter=\",\")\n",
    "    if(T):\n",
    "        adata = adata.T\n",
    "    #Filtering the genes which not express in all cells\n",
    "    sc.pp.filter_genes(adata, min_cells=1)\n",
    "    #Preprocessing the scRNA-seq data\n",
    "    sc.pp.normalize_total(adata, target_sum=1e4)\n",
    "    sc.pp.log1p(adata)\n",
    "    #Obtaining the HVGs\n",
    "    sc.pp.highly_variable_genes(adata, n_top_genes=topgenes)\n",
    "    adata.raw = adata\n",
    "    adata = adata[:, adata.var.highly_variable]\n",
    "\n",
    "    X = adata.X  #getting the gene expression matrix\n",
    "    # Converting the gene expression matrix into sub-vectors\n",
    "    #(n_cells,n_genes) -> (n_cells,gap_num,gap)  gap_num = int(gene_num / gap) + 1\n",
    "    single_cell_list = []\n",
    "    for single_cell in X:\n",
    "        feature = []\n",
    "        length = len(single_cell)\n",
    "        #spliting the gene expression vector into some sub-vectors whose length is gap\n",
    "        for k in range(0, length, gap):\n",
    "            if (k + gap > length):\n",
    "                a = single_cell[length - gap:length]\n",
    "            else:\n",
    "                a = single_cell[k:k + gap]\n",
    "\n",
    "            #scaling each sub-vectors \n",
    "            a = preprocessing.scale(a,axis=0, with_mean=True, with_std=True, copy=True)\n",
    "            feature.append(a)\n",
    "\n",
    "        feature = np.asarray(feature)\n",
    "        single_cell_list.append(feature)\n",
    "\n",
    "    single_cell_list = np.asarray(single_cell_list) #(n_cells,gap_num,gap)\n",
    "    \n",
    "    ##Obtaining cell types and all cell type labels of scRNA-seq data\n",
    "    if(label_path != None):\n",
    "        y_train = pd.read_csv(label_path)\n",
    "        y_train = y_train.T\n",
    "        y_train = y_train.values[0]\n",
    "\n",
    "        cell_types = []\n",
    "        labelss = []\n",
    "        for i in y_train:\n",
    "            i = str(i).upper()\n",
    "            if (not cell_types.__contains__(i)):\n",
    "                cell_types.append(i)\n",
    "            labelss.append(cell_types.index(i))\n",
    "\n",
    "        labelss = np.asarray(labelss)\n",
    "\n",
    "        return single_cell_list, labelss, cell_types\n",
    "    else:\n",
    "        return single_cell_list\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "6e50b1dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "class myDataSet(Dataset):\n",
    "    def __init__(self, data, label):\n",
    "        self.data = data\n",
    "        self.label = label\n",
    "        self.length = len(data)\n",
    "\n",
    "    def __len__(self):\n",
    "        return self.length\n",
    "\n",
    "    def __getitem__(self, index):\n",
    "        data = torch.from_numpy(self.data)\n",
    "        label = torch.from_numpy(self.label)\n",
    "\n",
    "        return data[index], label[index]\n",
    "\n",
    "\n",
    "##Positional Encoder Layer\n",
    "class PositionalEncoding(nn.Module):\n",
    "    def __init__(self, d_model, dropout=0.1, max_len=5000):\n",
    "        super(PositionalEncoding, self).__init__()\n",
    "        self.dropout = nn.Dropout(p=dropout)\n",
    "        pe = torch.zeros(max_len, d_model)\n",
    "        position = torch.arange(0, max_len, dtype=torch.float).unsqueeze(1)\n",
    "        div_term = torch.exp(torch.arange(0, d_model, 2).float() * (-math.log(10000.0) / d_model))\n",
    "        \n",
    "        ##the sine function is used to represent the odd-numbered sub-vectors\n",
    "        pe[:, 0::2] = torch.sin(position * div_term)\n",
    "        ##the cosine function is used to represent the even-numbered sub-vectors\n",
    "        pe[:, 1::2] = torch.cos(position * div_term)\n",
    "        pe = pe.unsqueeze(0).transpose(0, 1)\n",
    "        self.register_buffer('pe', pe)\n",
    "\n",
    "    def forward(self, x):\n",
    "        x = x + self.pe[:x.size(0), :]\n",
    "        return self.dropout(x)\n",
    "    \n",
    "##CIForm\n",
    "##function\n",
    "    #annotating cell type identification of scRNA-seq data\n",
    "##parameters\n",
    "    #input_dim  :Default is equal to gap\n",
    "    #nhead      :Number of heads in the attention mechanism\n",
    "    #d_model    :Default is equal to gap\n",
    "    #num_classes:Number of cell types\n",
    "    #dropout    :dropout rate which is used to prevent model overfitting\n",
    "class CIForm(nn.Module):\n",
    "    def __init__(self, input_dim, nhead=2, d_model=80, num_classes=2, dropout=0.1):\n",
    "        super().__init__()\n",
    "        #TransformerEncoderLayer with self-attention \n",
    "        self.encoder_layer = nn.TransformerEncoderLayer(\n",
    "            d_model=d_model, dim_feedforward=1024, nhead=nhead, dropout=dropout\n",
    "        )\n",
    "        \n",
    "        #Positional Encoding with self-attention \n",
    "        self.positionalEncoding = PositionalEncoding(d_model=d_model, dropout=dropout)\n",
    "        \n",
    "        #Classification layer\n",
    "        self.pred_layer = nn.Sequential(\n",
    "            nn.Linear(d_model, d_model),\n",
    "            nn.ReLU(),\n",
    "            nn.Linear(d_model, num_classes)\n",
    "        )\n",
    "\n",
    "    def forward(self, mels):\n",
    "        out = mels.permute(1, 0, 2)\n",
    "        #Positional Encoding layer\n",
    "        out = self.positionalEncoding(out)\n",
    "        #Transformer Encoder layer layer\n",
    "        out = self.encoder_layer(out)\n",
    "        out = out.transpose(0, 1)\n",
    "        #Pooling layer\n",
    "        out = out.mean(dim=1)\n",
    "        #Classification layer\n",
    "        out = self.pred_layer(out)\n",
    "        return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69f627b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import KFold, StratifiedKFold\n",
    "from sklearn.metrics import accuracy_score,f1_score,recall_score,precision_score\n",
    "from tqdm.auto import tqdm\n",
    "\n",
    "##main\n",
    "##parameters\n",
    "    #s                  :the length of a sub-vector\n",
    "    #referece_datapath  :the path of referece dataset\n",
    "    #query_datapath     :the path pf test dataset\n",
    "    #label_path         :the path of cell type annotation file(.csv)\n",
    "    #T                 \n",
    "    #                  T==True：The rows represent the cells and the columns represent the genes\n",
    "    #                  T=False: The rows represent the genes and the columns represent the cells\n",
    "def main(referece_datapath,label_path,query_datapath,s,T = False):\n",
    "    gap = s                        #the length of a sub-vector\n",
    "    topgenes = 2000                #the number of HVGs\n",
    "    \n",
    "    #Getting the data which input into the CIForm and the cell type annotations\n",
    "    Data, labels, cell_types = getData(gap, referece_datapath,topgenes,label_path=label_path,T = False)\n",
    "    \n",
    "    #Set parameters of CIForm\n",
    "    d_models = s   \n",
    "    heads = 64                     #the number of heads in self-attention mechanism\n",
    "    num_classes = len(cell_types)  #the number of cell types\n",
    "\n",
    "    lr = 0.0001                    #learning rate\n",
    "    dp = 0.1                       #dropout rate\n",
    "    n_epochs = 20                  #the number of epoch\n",
    "    \n",
    "    #Constructing the CIForm model\n",
    "    model = Classifier(input_dim=d_models, nhead=heads, d_model=d_models,\n",
    "                       num_classes=num_classes,dropout=dp)\n",
    "    #Setting loss function\n",
    "    criterion = nn.CrossEntropyLoss()\n",
    "    #Setting optimization function\n",
    "    optimizer = torch.optim.Adam(model.parameters(), lr=lr, weight_decay=1e-5)\n",
    "\n",
    "    l = Data.shape[0]\n",
    "    #Setting the batch size varies according to the number of cells\n",
    "    if(l > 5000):\n",
    "        batch_sizes = 512   #the size of batch\n",
    "    else:\n",
    "        batch_sizes = 256   #the size of batch\n",
    "\n",
    "\n",
    "    #the number of cell types\n",
    "    num_classes = len(cell_types)\n",
    "    \n",
    "    #Model verification by 5-KFold\n",
    "    skf = StratifiedKFold(n_splits=5,random_state=2021, shuffle=True)\n",
    "    fold = 0\n",
    "    Indexs = []\n",
    "    for index in range(len(Data)):\n",
    "        Indexs.append(index)\n",
    "    Indexs = np.asarray(Indexs)\n",
    "    for train_index, test_index in skf.split(Data, labels):\n",
    "        fold = fold + 1\n",
    "        X_train, X_test = Data[train_index], Data[test_index]\n",
    "        X_train = np.asarray(X_train)\n",
    "        X_test = np.asarray(X_test)\n",
    "\n",
    "        y_train, y_test = labels[train_index], labels[test_index]\n",
    "\n",
    "        y_train = np.asarray(y_train)\n",
    "        y_test = np.asarray(y_test)\n",
    "        \n",
    "        #Setting the training dataset\n",
    "        train_dataset = myDataSet(data=X_train,label=y_train)\n",
    "        train_loader = DataLoader(train_dataset,batch_size=batch_sizes,shuffle=True,pin_memory=True)\n",
    "        \n",
    "        #Setting the test dataset\n",
    "        test_dataset = myDataSet(data=X_test,label=y_test)\n",
    "        test_loader = DataLoader(test_dataset,batch_size=batch_sizes,shuffle=False,pin_memory=True)\n",
    "\n",
    "        #startting training CIForm.Using training data to train CIForm\n",
    "        #n_epochs: the times of Training \n",
    "        model.train()\n",
    "        for epoch in range(n_epochs):\n",
    "            # model.train()\n",
    "            # These are used to record information in training.\n",
    "            train_loss = []\n",
    "            train_accs = []\n",
    "            train_f1s = []\n",
    "            for batch in tqdm(train_loader):\n",
    "                # A batch consists of scRNA-seq data and corresponding cell type annotations.\n",
    "                data, labels = batch\n",
    "                logits = model(data)\n",
    "                labels = torch.tensor(labels, dtype=torch.long)\n",
    "                loss = criterion(logits, labels)\n",
    "                optimizer.zero_grad()\n",
    "                loss.backward()\n",
    "                optimizer.step()\n",
    "                #Getting the predicted cell type\n",
    "                preds = logits.argmax(1)\n",
    "                preds = preds.cpu().numpy()\n",
    "                labels = labels.cpu().numpy()\n",
    "                #Metrics\n",
    "                acc = accuracy_score(labels, preds)\n",
    "                f1 = f1_score(labels,preds,average='macro')\n",
    "                train_loss.append(loss.item())\n",
    "                train_accs.append(acc)\n",
    "                train_f1s.append(f1)\n",
    "            train_loss = sum(train_loss) / len(train_loss)\n",
    "            train_acc = sum(train_accs) / len(train_accs)\n",
    "            train_f1 = sum(train_f1s) / len(train_f1s)\n",
    "\n",
    "            print(f\"[ Train | {epoch + 1:03d}/{n_epochs:03d} ] loss = {train_loss:.5f}, acc = {train_acc:.5f}, f1 = {train_f1:.5f}\")\n",
    "            \n",
    "            ##Start the validation model, which predicts the cell types in the test dataset\n",
    "            model.eval()\n",
    "            test_accs = []\n",
    "            test_f1s = []\n",
    "            y_predict = []\n",
    "            labelss = []\n",
    "            for batch in tqdm(test_loader):\n",
    "                # A batch consists of scRNA-seq data and corresponding cell type annotations.\n",
    "                data, labels = batch\n",
    "                with torch.no_grad():\n",
    "                    logits = model(data)\n",
    "                    \n",
    "                #Getting the predicted cell type\n",
    "                preds = logits.argmax(1)\n",
    "                preds = preds.cpu().numpy().tolist()\n",
    "                labels = labels.cpu().numpy().tolist()\n",
    "                \n",
    "                #Metrics\n",
    "                acc = accuracy_score(labels, preds)\n",
    "                f1 = f1_score(labels, preds, average='macro')\n",
    "                test_f1s.append(f1)\n",
    "                test_accs.append(acc)\n",
    "                \n",
    "                y_predict.extend(preds)\n",
    "                labelss.extend(labels)\n",
    "            test_acc = sum(test_accs) / len(test_accs)\n",
    "            test_f1 = sum(test_f1s) / len(test_f1s)\n",
    "            print(\"---------------------------------------------end test---------------------------------------------\")\n",
    "            #Metrics\n",
    "            all_acc = accuracy_score(labelss, y_predict)\n",
    "            all_f1 = f1_score(labelss, y_predict, average='macro')\n",
    "            print(\"all_acc:\", all_acc,\"all_f1:\", all_f1)\n",
    "\n",
    "            labelsss = []\n",
    "            y_predicts = []\n",
    "            for i in labelss:\n",
    "                labelsss.append(cell_types[i])\n",
    "            for i in y_predict:\n",
    "                y_predicts.append(cell_types[i])\n",
    "            \n",
    "            \n",
    "            #Storing predicted cell types and the CIForm\n",
    "            log_dir = \"log/\"\n",
    "            if (not os.path.isdir(log_dir)):\n",
    "                os.makedirs(log_dir)\n",
    "            \n",
    "            np.save(log_dir  + 'y_tests.npy', labelsss)\n",
    "            np.save(log_dir  + 'y_predicts.npy', y_predicts)\n",
    "            torch.save(model.state_dict(), log_dir + 'CIForm.tar')\n",
    "\n",
    "            with open(log_dir + \"resilt.txt\", \"a\") as f:\n",
    "                f.writelines(\"acc:\" + str(all_acc) + \"\\n\")\n",
    "                f.writelines('f1:' + str(all_f1) + \"\\n\")\n",
    "\n",
    "referece_datapath = \"Dataset\\scRNAseq_Benchmark_datasets\\Intra-dataset\\Pancreatic_data\\Human\\Human.csv\" \n",
    "label_path = \"Dataset\\scRNAseq_Benchmark_datasets\\Intra-dataset\\Pancreatic_data\\Human\\Labels.csv\"\n",
    "query_datapath = \"\"\n",
    "s = 1024  # the length of a sub-vector\n",
    "##parameters\n",
    "    #s                  :the length of a sub-vector\n",
    "    #referece_datapath  :the path of referece dataset\n",
    "    #query_datapath     :the path pf test dataset\n",
    "    #label_path         :the path of cell type annotation file(.csv)\n",
    "    #T                 \n",
    "    #                  T==True：The rows represent the cells and the columns represent the genes\n",
    "    #                  T=False: The rows represent the genes and the columns represent the cells\n",
    "main(referece_datapath,label_path,query_datapath,s,T = False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}