04_LGBM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn import datasets\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import roc_auc_score, accuracy_score\n",
    "from sklearn.inspection import permutation_importance\n",
    "from sklearn.metrics import r2_score\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "import lightgbm as lgb\n",
    "import math\n",
    "import scipy\n",
    "import conorm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_fasta(fasta_file):\n",
    "    content = {}\n",
    "    fas = open(fasta_file,'r')\n",
    "    seq = ''\n",
    "    last_entry = ''\n",
    "    while True:\n",
    "        line = fas.readline()\n",
    "        if len(line)==0:\n",
    "            content[last_entry]=seq\n",
    "            break\n",
    "        if '>' in line:\n",
    "            if last_entry=='':\n",
    "                line = line.rstrip()\n",
    "                last_entry=line\n",
    "            else:\n",
    "                line = line.rstrip()\n",
    "                content[last_entry]=seq\n",
    "                last_entry=line\n",
    "                seq=''\n",
    "        else:\n",
    "            line = line.rstrip()\n",
    "            seq += line\n",
    "    return content\n",
    "\n",
    "def read_gff_complete(gff_file):\n",
    "    import pandas as pd\n",
    "    columns=['SeqID','Source','type','Start',\n",
    "             'End','Score','Strand','Phase']\n",
    "    df_raw = []\n",
    "    retrieve_keys = []\n",
    "    additional_info = []\n",
    "    revised_info = []\n",
    "    with open(gff_file, 'r') as f:\n",
    "        while True:\n",
    "            line = f.readline().rstrip()\n",
    "            if line == '':\n",
    "                f.close()\n",
    "                break\n",
    "            elif not line.startswith('#'):\n",
    "                line = line.split('\\t')\n",
    "                df_raw.append(line[:-1])\n",
    "                info_dict = {}\n",
    "                for x in line[-1].split(';'):\n",
    "                    if x.count('=') == 1:\n",
    "                        k, v = x.split('=')\n",
    "                        info_dict[k] = v\n",
    "                additional_info.append(info_dict)\n",
    "                retrieve_keys += list(info_dict.keys())\n",
    "    retrieve_keys = list(np.unique(retrieve_keys))\n",
    "    for t in additional_info:\n",
    "        entry_list = []\n",
    "        for k in retrieve_keys:\n",
    "            if k in t:\n",
    "                entry_list.append(t[k])\n",
    "            else:\n",
    "                entry_list.append('-')\n",
    "        revised_info.append(entry_list)\n",
    "    df = pd.DataFrame(df_raw,columns=columns)\n",
    "    df[['Start','End']] = df[['Start','End']].astype(int)\n",
    "    df[retrieve_keys] = revised_info\n",
    "    return df.sort_values(by='Start').reset_index(drop=True)\n",
    "\n",
    "def reverse_complement(s):\n",
    "    rc_dict = {'A':'T','T':'A','C':'G','G':'C','a':'t','t':'a','c':'g','g':'c'}\n",
    "    comp = [rc_dict[x] for x in list(s)]\n",
    "    comp.reverse()\n",
    "    return ''.join(comp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature compilation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "tss = pd.read_csv('./data/LGBM/ref/MTB TSS+-35bp sequence.csv')\n",
    "data = pd.read_csv('./data/RPKM Mtb.csv').set_index('Unnamed: 0')\n",
    "sample = pd.read_excel('./data/sample/Sample Mtb.xlsx').set_index('Run')\n",
    "gff = read_gff_complete('./ref/Mtb_GCF_000195955.2_ASM19595v2_genomic.gff')\n",
    "mtb_g = read_fasta('./data/LGBM/ref/Mycobacterium_tuberculosis_H37Rv_genome_v4.fasta')\n",
    "gseq = mtb_g[list(mtb_g.keys())[0]]\n",
    "info_df = pd.read_csv(\"./data/LGBM/ref/TP and feature Mtb_LGBM format.csv\").set_index('gene')\n",
    "tf = pd.read_excel('./data/LGBM/ref/TF_Chipseq_ncomms6829_MOESM46_ESM.xlsx',sheet_name='same_as_sheet1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "TF_cutoff = 20\n",
    "filtered = []\n",
    "counts = tf.groupby('Regulator').count()\n",
    "for x,c in zip(counts.index, counts['Gene']):\n",
    "    if c>=TF_cutoff:\n",
    "        filtered.append(tf[tf['Regulator']==x])\n",
    "filtered = pd.concat(filtered)\n",
    "regulators = filtered['Regulator'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "gsize=len(gseq)\n",
    "tf_table = []\n",
    "tss_plus = np.sort(tss[tss['strand']=='+']['TSS.identifier'].values)\n",
    "tss_minus = np.sort(tss[tss['strand']=='-']['TSS.identifier'].values)\n",
    "c = 0\n",
    "tp = []\n",
    "nearest_tss = []\n",
    "start_codon = []\n",
    "genes = []\n",
    "strand_sign = []\n",
    "exp = []\n",
    "gene_center = []\n",
    "strand_list = []\n",
    "for x in info_df.index:\n",
    "    if x in gff['locus_tag'].values:\n",
    "        start,end,strand = gff[gff['locus_tag']==x][['Start','End','Strand']].values[0]\n",
    "        if strand == '+':\n",
    "            start_codon.append(start)\n",
    "            strand_sign.append(1)\n",
    "            if start <100:\n",
    "                gtss = tss_plus[-1]\n",
    "            else:\n",
    "                gtss = tss_plus[np.where((tss_plus-start)<10)[0][-1]]\n",
    "        else:\n",
    "            start_codon.append(end)\n",
    "            gtss = tss_minus[np.where((tss_minus-end)>-10)[0][0]]\n",
    "            strand_sign.append(-1)\n",
    "        strand_list.append(strand)\n",
    "        nearest_tss.append(gtss)\n",
    "        genes.append(x)\n",
    "        tp.append(info_df.loc[x,'TP'])\n",
    "        exp.append(info_df.loc[x,'Median_logRPKM'])\n",
    "        gene_center.append((start+end)/2)\n",
    "    else:\n",
    "        print(x)\n",
    "tss_df = pd.DataFrame({'TP':tp,'TSS':nearest_tss,\n",
    "                       'Start_codon':start_codon,\n",
    "                       'Strand_sign':strand_sign,\n",
    "                       'Expression':exp,'Center':gene_center,'Strand':strand_list},index=genes)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Retieve promoter sequencing feature from MathFeature server"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "up = 60\n",
    "down = 20\n",
    "\n",
    "f = open('./data/LGBM/ref/tss_seq.fa','w')\n",
    "pro_sequences = []\n",
    "\n",
    "tsites = []\n",
    "for i,(tsite,s) in enumerate(tss_df[['TSS','Strand']].values):\n",
    "    if tsite not in tsites:\n",
    "        tsites.append(tsite)\n",
    "        g = tss_df.index[i]\n",
    "        if s=='+':\n",
    "            seq = gseq[(tsite-up):(tsite+down)]\n",
    "        else:\n",
    "            seq = reverse_complement(gseq[(tsite-down-1):(tsite+up-1)])\n",
    "        pro_sequences.append(seq)\n",
    "        f.write('>TSS_{}\\n'.format(tsite))\n",
    "        f.write('{}\\n\\n'.format(seq))\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "f1 = pd.read_csv('./data/LGBM/ref/tss80_fourier_ext.csv')\n",
    "f1['label'] = [int(x.split('_')[1]) for x in f1['nameseq'].values]\n",
    "f1 = f1.set_index('label')\n",
    "f1_cols = f1.columns[1:]\n",
    "tss_df[['MF_{}'.format(x) for x in f1_cols]] = f1.loc[tss_df['TSS'].values][f1_cols].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3891, 26)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tss_df.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TF binding feature compilation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n",
      "/var/folders/1y/_2nhqmfd56s2k1hktst6p9b80000gn/T/ipykernel_70874/4258671928.py:21: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
      "  tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;\n"
     ]
    }
   ],
   "source": [
    "basic_features_cols = info_df.columns[17:]\n",
    "tss_df[basic_features_cols] = info_df.loc[tss_df.index,basic_features_cols].values;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "101"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(basic_features_cols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3891, 127)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tss_df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['MF_median', 'MF_maximum', 'MF_minimum', 'MF_peak',\n",
       "       'MF_none_levated_peak', 'MF_sample_standard_deviation',\n",
       "       'MF_population_standard_deviation', 'MF_percentile15',\n",
       "       'MF_percentile25', 'MF_percentile50',\n",
       "       ...\n",
       "       'l_beta', 'u_beta', 'beta_max_span', 'M', 'l_M', 'u_M', 'M_span',\n",
       "       'Vulnerability.Index', 'VI.Lower.Bound', 'VI.Upper.Bound'],\n",
       "      dtype='object', length=119)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tss_df.columns[8:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "X = tss_df[tss_df.columns[8:]].values\n",
    "X[np.isnan(X)]=0\n",
    "Y = np.log2(tss_df['TP'].values).reshape(-1,1) # use log2-transformed TP instead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "def min_max(data):\n",
    "    normd = (data-data.min(axis=0)[np.newaxis,:])/(data.max(axis=0)-data.min(axis=0))[np.newaxis,:]\n",
    "    return normd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_norm = min_max(X)\n",
    "Y_norm = min_max(Y).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# initiate importance dict\n",
    "importance_dict = pd.DataFrame({'feature_name':tss_df.columns[8:]})\n",
    "\n",
    "# record fit stats\n",
    "R = []\n",
    "MSE = []\n",
    "\n",
    "# iterative training\n",
    "iter_n = 100\n",
    "test_size = 0.4\n",
    "num_leaves = 30\n",
    "learning_rate = 0.015\n",
    "feature_fraction=0.8\n",
    "n_estimators = 200\n",
    "max_depth=15\n",
    "max_bin=10\n",
    "bagging_fraction=1\n",
    "\n",
    "\n",
    "for i in range(0,iter_n):\n",
    "    X_train, X_test, y_train, y_test = train_test_split(X_norm, Y_norm, \n",
    "                                                        test_size=test_size, \n",
    "                                                        random_state=i)\n",
    "    lgb_train = lgb.Dataset(X_train, y_train)\n",
    "    lgb_eval = lgb.Dataset(X_test, y_test, reference=lgb_train)\n",
    "    gbm = lgb.LGBMRegressor(objective='regression', reg_lambda=1,\n",
    "                            max_bin=max_bin,\n",
    "                            num_leaves=num_leaves,\n",
    "                            learning_rate=learning_rate, \n",
    "                            verbose=0,\n",
    "                            n_estimators=n_estimators,\n",
    "                            max_depth=max_depth,\n",
    "                            n_jobs=16)\n",
    "    \n",
    "    gbm.fit(X_train, y_train, \n",
    "            eval_set=[(X_test, y_test)], \n",
    "            eval_metric='l2')\n",
    "    \n",
    "    y_pred = gbm.predict(X_test, num_iteration=gbm.best_iteration_)\n",
    "    y_train_eval = gbm.predict(X_train, num_iteration=gbm.best_iteration_)\n",
    "    feature_importance = gbm.feature_importances_\n",
    "    \n",
    "    importance_dict_iter = pd.DataFrame({'feature_name':tss_df.columns[8:],\n",
    "                                         i:feature_importance})\n",
    "    importance_dict = pd.merge(importance_dict,importance_dict_iter,how='left',on='feature_name')\n",
    "\n",
    "    R.append([scipy.stats.pearsonr(y_test, y_pred)[0],\n",
    "              scipy.stats.pearsonr(y_train, y_train_eval)[0]])\n",
    "    MSE.append(mean_squared_error(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "TP                     0.635053\n",
       "TSS                       52810\n",
       "Start_codon               52831\n",
       "Strand_sign                   1\n",
       "Expression             8.756984\n",
       "                         ...   \n",
       "u_M                       0.582\n",
       "M_span                    0.417\n",
       "Vulnerability.Index      -0.696\n",
       "VI.Lower.Bound           -2.858\n",
       "VI.Upper.Bound            1.664\n",
       "Name: Rv0049, Length: 127, dtype: object"
      ]
     },
     "execution_count": 159,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tss_df.loc['Rv0049']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['operon_length', 'basePercent_GC', 'Activator.number', 'width',\n",
       "       'basePercent_G', 'MF_maximum', 'A3s', 'AA_C', 'AA_A',\n",
       "       'maxCov.TSS.site', 'distance_to_ori', 'AA_R', 'MF_amplitude',\n",
       "       'X2.5..quantile.κ', 'u_beta', 'l_gamma', 'AA_S', 'AA_Y',\n",
       "       'basePercent_T', 'basePercent_CT', 'AA_L', 'Repressor.number',\n",
       "       'X97.5..quantile.κ', 'Nc', 'AA_H', 'Rep.Direction', 'MF_minimum',\n",
       "       'basePercent_C', 'MF_kurtosis', 'MF_percentile25', 'MF_peak',\n",
       "       'MF_sample_standard_deviation', 'AA_D', 'MF_none_levated_peak',\n",
       "       'AA_E', 'C3s', 'maxCov.TSS.MeanCoverage', 'MF_percentile15',\n",
       "       'MF_skewness', 'AA_G', 'Geometric.mean.ω', 'AA_P', 'AA_F',\n",
       "       'maxCov.TSS.distance2gene', 'G3s', 'MF_median', 'Aromo', 'GC3s',\n",
       "       'MF_percentile75', 'AA_I', 'AA_V', 'AA_K', 'AA_M', 'AA_Q',\n",
       "       'n_guides', 'T3s', 'Mean.Pr.ω.1.', 'l_M', 'X97.5..quantile.θ',\n",
       "       'u_H', 'AA_N', 'X2.5..quantile.geometric.mean.ω', 'basePercent_A',\n",
       "       'MF_interquartile_range', 'AA_T', 'CAI', 'AA_W', 'Gravy',\n",
       "       'X97.5..quantile.geometric.mean.ω', 'Fop', 'X..sites.Pr.ω.1..90.',\n",
       "       'l_H', 'l_beta', 'str_span', 'CBI', 'Geometric.mean.θ', 'beta_max',\n",
       "       'beta_max_span', 'K', 'operon_innerLoc', 'X2.5..quantile.θ',\n",
       "       'Vulnerability.Index', 'Geometric.mean.κ', 'u_K', 'gene_gamma',\n",
       "       'VI.Lower.Bound', 'VI.Upper.Bound', 'l_K', 'u_M', 'M', 'u_gamma',\n",
       "       'H', 'Leaderless.TSS.number', 'M_span', 'maxCov.TSS.base_G',\n",
       "       'maxCov.TSS.base_A', 'strand', 'stopcodon_TGA', 'startcodon_GTG',\n",
       "       'startcodon_ATG', 'Leadered.TSS.number', 'maxCov.TSS.leaderseq',\n",
       "       'stopcodon_TAG', 'Rep.Transcript.Direction', 'stopcodon_TAA',\n",
       "       'startcodon_TTG', 'MF_coefficient_of_variation',\n",
       "       'ActivatingTarget.number', 'RepressingTarget.number',\n",
       "       'MF_variance', 'MF_population_standard_deviation',\n",
       "       'maxCov.TSS.base_C', 'startcodon_ATC', 'startcodon_ATA',\n",
       "       'MF_percentile50', 'MF_semi_interquartile_range', 'startcodon_ATT',\n",
       "       'startcodon_CTG', 'maxCov.TSS.base_T'], dtype=object)"
      ]
     },
     "execution_count": 154,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = importance_dict['feature_name'].values[np.flip(np.argsort(np.mean(importance_dict.values[:,1:],axis=1)))]\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [],
   "source": [
    "importance_dict.to_csv('./data/LGBM/ref/LGBM feature importance.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(\"./data/LGBM/ref/Prediction_LGBMR_TP_r.txt\", np.array(R), fmt='%f')\n",
    "np.savetxt(\"./data/LGBM/ref/Prediction_LGBMR_TP_mse.txt\", np.array(MSE), fmt='%f')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "tss_df[['MF_{}'.format(x) for x in f1_cols]].to_excel('./data/LGBM/ref/MathFeatureBinaryFourier.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0                 MF_median\n",
       "1                MF_maximum\n",
       "2                MF_minimum\n",
       "3                   MF_peak\n",
       "4      MF_none_levated_peak\n",
       "               ...         \n",
       "114                     u_M\n",
       "115                  M_span\n",
       "116     Vulnerability.Index\n",
       "117          VI.Lower.Bound\n",
       "118          VI.Upper.Bound\n",
       "Name: feature_name, Length: 119, dtype: object"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "importance_dict['feature_name']"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}