segment_convergence_outlier_detection.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a1fbbcb5-b920-4c42-a0fe-bbfae9863188",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "import pickle,sys,os\n",
    "import antropy as ant\n",
    "from numpy import *\n",
    "import heapq, qnorm\n",
    "import ruptures as rpt\n",
    "from scipy.stats import norm, pearsonr, zscore, rankdata\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.covariance import EllipticEnvelope\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "dc27512a-b707-4303-a2df-93098867cf0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('../../Human_hg38_chrom_ends.pkl','rb') as a:\n",
    "    ENDS = pickle.load(a)\n",
    "\n",
    "with open('../../Human_hg38_auto_10kb_windows.pkl','rb') as b:\n",
    "    W = pickle.load(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "83da0f94-9ae8-4244-be1a-1d50c8c3e81e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_numpy_arr(infile,columns=None):\n",
    "    arr_file = load(infile)\n",
    "    if columns is None:\n",
    "        arr = arr_file['arr']\n",
    "    else:\n",
    "        arr = arr_file['arr'][:,columns]\n",
    "    for cols in range(arr.shape[1]):\n",
    "        csum = sum(arr[:,cols])\n",
    "        arr[:,cols] /= csum\n",
    "    return arr\n",
    "\n",
    "def get_full_table(arr,win=10000):\n",
    "    rows = 0\n",
    "    full_table = []\n",
    "    for ch in range(22):\n",
    "        chrom = 'chr'+str(ch+1)\n",
    "        bound = ENDS[chrom]\n",
    "        windows = bound // win\n",
    "        full_table.append(arr[rows:rows+windows,:])\n",
    "        rows += windows\n",
    "    if rows!=arr.shape[0]:\n",
    "        print('error')\n",
    "    return full_table\n",
    "\n",
    "def get_full_table_k(full_table,win=10000,k=1):\n",
    "    if k == 1:\n",
    "        return full_table\n",
    "    else:\n",
    "        ft = []\n",
    "        for ch in range(22):\n",
    "            chrom = 'chr'+str(ch+1)\n",
    "            bound = ENDS[chrom]\n",
    "            win2 = win*k\n",
    "            windows2 = bound // win2\n",
    "            arr_i = full_table[ch]\n",
    "            arr_k = zeros((windows2,arr_i.shape[1]))\n",
    "            for j in range(windows2-1):\n",
    "                arr_k[j] = arr_i[j*k:(j+1)*k,:].sum(axis=0)\n",
    "            arr_k[-1] = arr_i[(windows2-1)*k:,:].sum(axis=0)\n",
    "            ft.append(arr_k)\n",
    "        return ft\n",
    "\n",
    "def get_full_table_array_k(full_table,win=10000,k=1):\n",
    "    if k == 1:\n",
    "        for ch in range(22):\n",
    "            if ch == 0:\n",
    "                arr = full_table[ch]\n",
    "            else:\n",
    "                arr = concatenate((arr,full_table[ch]),axis=0)\n",
    "    else:\n",
    "        ft = []\n",
    "        for ch in range(22):\n",
    "            chrom = 'chr'+str(ch+1)\n",
    "            bound = ENDS[chrom]\n",
    "            win2 = win*k\n",
    "            windows2 = bound // win2\n",
    "            arr_i = full_table[ch]\n",
    "            arr_k = zeros((windows2,arr_i.shape[1]))\n",
    "            for j in range(windows2-1):\n",
    "                arr_k[j] = arr_i[j*k:(j+1)*k,:].sum(axis=0)\n",
    "            arr_k[-1] = arr_i[(windows2-1)*k:,:].sum(axis=0)\n",
    "            if ch==0:\n",
    "                arr = arr_k\n",
    "            else:\n",
    "                arr = concatenate((arr,arr_k),axis=0)\n",
    "            #ft.append(arr_k)\n",
    "    return arr\n",
    "\n",
    "def get_full_table_array_k_chr(full_table,win=10000,k=1):\n",
    "    ends = []\n",
    "    if k == 1:\n",
    "        for ch in range(22):\n",
    "            if ch == 0:\n",
    "                arr = full_table[ch]\n",
    "            else:\n",
    "                arr = concatenate((arr,full_table[ch]),axis=0)\n",
    "            ends.append(arr.shape[0])\n",
    "    else:\n",
    "        ft = []\n",
    "        for ch in range(22):\n",
    "            chrom = 'chr'+str(ch+1)\n",
    "            bound = ENDS[chrom]\n",
    "            win2 = win*k\n",
    "            windows2 = bound // win2\n",
    "            arr_i = full_table[ch]\n",
    "            arr_k = zeros((windows2,arr_i.shape[1]))\n",
    "            for j in range(windows2-1):\n",
    "                arr_k[j] = arr_i[j*k:(j+1)*k,:].sum(axis=0)\n",
    "            arr_k[-1] = arr_i[(windows2-1)*k:,:].sum(axis=0)\n",
    "            if ch==0:\n",
    "                arr = arr_k\n",
    "            else:\n",
    "                arr = concatenate((arr,arr_k),axis=0)\n",
    "            ends.append(arr.shape[0])\n",
    "    return arr, ends\n",
    "\n",
    "def get_rate_of_change(table):\n",
    "    for ch in range(22):\n",
    "        arr = table[ch]\n",
    "        if ch == 0:\n",
    "            new_arr = arr[1:,:] - arr[:-1,:]\n",
    "        else:\n",
    "            arr = arr[1:,:] - arr[:-1,:]\n",
    "            new_arr = concatenate((new_arr,arr),axis=0)\n",
    "    for cols in range(new_arr.shape[1]):\n",
    "        sigma = sqrt(sum([x**2 for x in new_arr[:,cols]])/new_arr.shape[0])\n",
    "        if sigma == 0:\n",
    "            print('samples with all zero values found')\n",
    "        new_arr[:,cols] /= sigma\n",
    "    print(new_arr.shape)\n",
    "    return new_arr\n",
    "\n",
    "def get_rate_of_change_chr(table,ch):\n",
    "    arr = table[ch]\n",
    "    new_arr = arr[1:,:] - arr[:-1,:]\n",
    "    for cols in range(new_arr.shape[1]):\n",
    "        sigma = sqrt(sum([x**2 for x in new_arr[:,cols]])/new_arr.shape[0])\n",
    "        new_arr[:,cols] /= sigma\n",
    "    #print(new_arr.shape)\n",
    "    return new_arr\n",
    "\n",
    "def get_arr_for_chrs(infile,columns=None,k=1):\n",
    "    df_arr = get_numpy_arr(infile,columns)\n",
    "    table = get_full_table(df_arr)\n",
    "    ftk = get_full_table_k(table,k=k)\n",
    "    tables = []\n",
    "    for ch in range(22):\n",
    "        arr = get_rate_of_change_chr(ftk,ch)\n",
    "        tables.append(arr)\n",
    "    return tables\n",
    "\n",
    "def get_arr(infile,columns=None,k=1):\n",
    "    df_arr = get_numpy_arr(infile,columns)\n",
    "    table = get_full_table(df_arr)\n",
    "    ftk, ends = get_full_table_array_k_chr(table,k=k)\n",
    "    #print(ends)\n",
    "    return ftk #arrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ac160fcc-3adc-4c20-a71a-a023a92fbef9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(287487, 3202)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k = 1\n",
    "thgs = get_arr('merged_table_no_correction_1000g_10kb.npz',k=k)\n",
    "thgs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9c95e9fe-2aab-4c8b-8caf-ac10944c1de8",
   "metadata": {},
   "outputs": [],
   "source": [
    "ctrl38s = get_arr('../../merged_table_no_correction_ctrl38_10kb.npz',k=k)\n",
    "ctdna36s = get_arr('../../merged_table_no_correction_ctdna36luad_10kb.npz',k=k)\n",
    "brcas = get_arr('../../merged_table_no_correction_brca_10kb.npz',k=k)\n",
    "crcs = get_arr('../../merged_table_no_correction_crc_10kb.npz',k=k)\n",
    "ctrl29s = get_arr('../../merged_table_no_correction_ctrl29_10kb.npz',k=k)\n",
    "ega8460 = get_arr('../../merged_table_no_correction_EGA8460_10kb.npz',k=k)\n",
    "ega5343_ctrl = get_arr('../../merged_table_no_correction_EGA5343_ctrl_10kb.npz',k=k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1cceb408-069f-4396-bf3f-32ac97886bb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "ega5343_prostate = get_arr('../../merged_table_no_correction_EGA5343_prostate_10kb.npz',k=k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1f498441-8715-44ad-addd-38b9e237b588",
   "metadata": {},
   "outputs": [],
   "source": [
    "arr_file = load('../../merged_table_no_correction_EGA7080_10kb.npz')\n",
    "samples = arr_file['samples']\n",
    "\n",
    "ctrl = [i for i in range(len(samples)) if samples[i].startswith('Ctrl')]\n",
    "ctdna = [i for i in range(len(samples)) if i not in ctrl] # and samples[i].endswith('_1')\n",
    "ega7080_ctrl = get_arr('../../merged_table_no_correction_EGA7080_10kb.npz',columns=ctrl,k=k)\n",
    "ega7080_ctdna = get_arr('../../merged_table_no_correction_EGA7080_10kb.npz',columns=ctdna,k=k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4c31ceb3-0d68-4ac9-afe8-63491fe69f0d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "51248\n",
      "[0, 21838, 43007, 61151, 77732, 93356, 108799, 122508, 135682, 146854, 158478, 170121, 181337, 190291, 197701, 205137, 211316, 216730, 222093, 226049, 231029, 233445, 236239]\n"
     ]
    }
   ],
   "source": [
    "ends = [24895, 49114, 68943, 87964, 106117, 123197,139131, 153644, 167483, 180862, \n",
    "        194370, 207697, 219133, 229837, 240036, 249069, 257394, 265431, 271292, 277736, 282406, 287487]\n",
    "feats = [i for i in range(thgs.shape[0]) if sum(thgs[i,:])>0]\n",
    "print(ends[-1]-len(feats))\n",
    "#zf = [i for i in range(thgs.shape[0]) if i not in feats]\n",
    "#print(len(zf))\n",
    "chr_splits = [0]\n",
    "j = 0\n",
    "for e in ends:\n",
    "    for i in range(j,len(feats)):\n",
    "        if feats[i] > e:\n",
    "            chr_splits.append(i)\n",
    "            j = i\n",
    "            break\n",
    "chr_splits.append(len(feats))\n",
    "print(chr_splits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "092d1455-ae19-4365-a8f9-735cc6704521",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_rank(arr):\n",
    "    ranked_arr = []\n",
    "    for i in range(arr.shape[1]):\n",
    "        ranked = rankdata(arr[:,i],'min')\n",
    "        ranked_arr.append(ranked)\n",
    "    return array(ranked_arr).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1a95c7d2-3931-482f-b11e-8b73a5596795",
   "metadata": {},
   "outputs": [],
   "source": [
    "rctrl38s = get_rank(ctrl38s[feats,:])\n",
    "rctdna36s = get_rank(ctdna36s[feats,:])\n",
    "rpost22s = get_rank(post22s[feats,:])\n",
    "rbrcas = get_rank(brcas[feats,:])\n",
    "rcrcs = get_rank(crcs[feats,:])\n",
    "rctrl29s = get_rank(ctrl29s[feats,:])\n",
    "rega5343_ctrl = get_rank(ega5343_ctrl[feats,:])\n",
    "rega7080_ctrl = get_rank(ega7080_ctrl[feats,:])\n",
    "rega8460 = get_rank(ega8460[feats,:])\n",
    "rega7080_ctdna = get_rank(ega7080_ctdna[feats,:])\n",
    "rthgs = get_rank(thgs[feats,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f1511658-ba06-49fa-8622-9c91e2233b98",
   "metadata": {},
   "outputs": [],
   "source": [
    "rega5343_prostate = get_rank(ega5343_prostate[feats,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "cd2a48e4-3335-4677-b8e7-04b64632adf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "rctrls = hstack((rega5343_ctrl,rega7080_ctrl,rctrl38s,rctrl29s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "691efc58-65dc-4954-a994-05bf7fa42cd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "rctdnas = hstack((rega7080_ctdna,rctdna36s,rbrcas,rcrcs,rega8460))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "a95b2e91-ad72-457c-912f-7611c6ae17f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# segmentation\n",
    "def get_breakpoints(arr=rthgs,splits=chr_splits):\n",
    "    bkps = {}\n",
    "    breaks = [0]\n",
    "    offset = 0\n",
    "    for ch in range(22):\n",
    "        ctrl_signals = mean(arr[splits[ch]:splits[ch+1],:],axis=1)\n",
    "        sigma = std(ctrl_signals)\n",
    "        n = len(ctrl_signals)\n",
    "        if n < 100:\n",
    "            offset+=n\n",
    "            continue\n",
    "        algo = rpt.Binseg(model=\"l1\",min_size=100).fit(ctrl_signals) # min size is 1Mbp # \"l1\", \"rbf\", \"linear\", \"normal\", \"ar\",...\n",
    "        my_bkps = algo.predict(pen=log(n)*sigma)\n",
    "        bkps['chr'+str(ch+1)] = [offset]+[x+offset for x in my_bkps]\n",
    "        breaks += [x+offset for x in my_bkps]\n",
    "        offset+= n\n",
    "    return bkps, breaks\n",
    "bkps, breaks = get_breakpoints()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "618b905c-e2b9-4c3e-aafc-8c03c891b3da",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_segment_mean(arr,breaks=breaks):\n",
    "    n = len(breaks)-1\n",
    "    arr_seg = []\n",
    "    for i in range(n):\n",
    "        arr_seg.append(mean(arr[breaks[i]:breaks[i+1],:],axis=0))\n",
    "    return array(arr_seg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "f4747093-49df-49f5-8119-f7062f42fb5f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14766.613891687324 9091.224919564898 8910.747007094713 5086.423075596531\n"
     ]
    }
   ],
   "source": [
    "sctrl38s = get_segment_mean(rctrl38s)\n",
    "sctrl29s = get_segment_mean(rctrl29s) \n",
    "sega7080_ctrl = get_segment_mean(rega7080_ctrl)\n",
    "sega5343_ctrl = get_segment_mean(rega5343_ctrl)\n",
    "print(median(std(sctrl38s,axis=1)),median(std(sctrl29s,axis=1)),median(std(sega7080_ctrl,axis=1)),median(std(sega5343_ctrl,axis=1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "380cbcc8-472c-4fa6-a606-caaedecb807a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('chr1', 98), ('chr2', 100), ('chr3', 78), ('chr4', 81), ('chr5', 75), ('chr6', 71), ('chr7', 67), ('chr8', 67), ('chr9', 48), ('chr10', 56), ('chr11', 56), ('chr12', 41), ('chr13', 32), ('chr14', 35), ('chr15', 36), ('chr16', 28), ('chr17', 20), ('chr18', 30), ('chr19', 21), ('chr20', 21), ('chr21', 12), ('chr22', 17)]\n"
     ]
    }
   ],
   "source": [
    "nsegments = []\n",
    "for j in range(1,23):\n",
    "    chrom = 'chr'+str(j)\n",
    "    nsegments.append((chrom,len(bkps[chrom])-1))\n",
    "print(nsegments)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc21e20c-fb36-4050-b7b8-c4d6d5c6b004",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "dffab232-5f84-4542-b030-5ba08720a1b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_chrom_mean(arr,splits=chr_splits):\n",
    "    n = len(splits)-1\n",
    "    arr_seg = []\n",
    "    for i in range(n):\n",
    "        arr_seg.append(mean(arr[splits[i]:splits[i+1],:],axis=0))\n",
    "    return array(arr_seg)\n",
    "\n",
    "def sort_arr_by_chr_mean(arr,chrom='chr1',splits=chr_splits):\n",
    "    i = int(chrom.split('hr')[1])-1\n",
    "    m = mean(arr[splits[i]:splits[i+1],:],axis=0)\n",
    "    sorted_idx = argsort(m)[::-1]\n",
    "    return sorted_idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "564a1b63-abd5-4cb5-b166-d7d7593998d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_rank_hist(arr, nbins=8):\n",
    "    hist = []\n",
    "    for s in range(arr.shape[1]):\n",
    "        h,_ = histogram(arr[:,s], bins=linspace(0,240000,nbins+1))\n",
    "        hist.append(h)\n",
    "    return array(hist).T\n",
    "\n",
    "def test_gaussian(arr,ctrl,error=1e-2): # feats x n\n",
    "    estimator = EllipticEnvelope(random_state=12345,contamination=error,support_fraction=1.0).fit(ctrl.T)\n",
    "    return estimator.predict(arr.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "9c50f436-1a7d-4ef5-ae80-3a45e33f8541",
   "metadata": {},
   "outputs": [],
   "source": [
    "def chr_hist_segment_classifier(arr,chrom,ctrls=rega7080_ctrl,error=1e-2,nbins=8,bps=bkps,splits=chr_splits):\n",
    "    sorted_idx = sort_arr_by_chr_mean(arr,chrom)\n",
    "    sorted_ctrl_idx = sort_arr_by_chr_mean(ctrls,chrom)\n",
    "    sarr = get_segment_mean(arr[:,sorted_idx],breaks=bps[chrom])\n",
    "    sctrl = get_segment_mean(ctrls[:,sorted_ctrl_idx],breaks=bps[chrom])\n",
    "    harr = get_rank_hist(sarr,nbins)#;print(mean(harr,axis=1).shape)\n",
    "    hctrl = get_rank_hist(sctrl,nbins)\n",
    "    arr_labels = test_gaussian(harr,hctrl,error=error)\n",
    "    return arr_labels[sorted_idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb4c7168-5a7f-4499-b36e-5eb1c878e162",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "51352d6d-daf7-4960-a6ba-0463fa1848e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ensemble_classifier_sample(arr,ctrls,incl=['chr1','chr2','chr3','chr4','chr5','chr6','chr7','chr12','chr13','chr15'],alpha=0.05,error=1e-64,nbins=8):#\n",
    "    arr_labels = zeros(arr.shape[1])\n",
    "    count = 0\n",
    "    for chrom in incl:#bkps:\n",
    "        nseg = len(bkps[chrom])-1\n",
    "        labels = array([max(x,0) for x in chr_hist_segment_classifier(arr,chrom,ctrls,error,nbins)])\n",
    "        labels *= nseg\n",
    "        arr_labels += labels\n",
    "        count += nseg\n",
    "    arr_labels/=count\n",
    "    return len([x for x in arr_labels if x<1-alpha])/len(arr_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "83acab07-5c33-4835-8a57-3b85345dc414",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ensemble_classifier(arr,ctrls,incl=['chr1','chr2','chr3','chr4','chr5','chr6','chr7','chr12','chr13','chr15'],alpha=0.05,error=1e-64,nbins=8,verbose=False):#\n",
    "    arr_labels = zeros(arr.shape[1])\n",
    "    count = 0\n",
    "    for chrom in incl:#bkps:\n",
    "        nseg = len(bkps[chrom])-1\n",
    "        labels = array([max(x,0) for x in chr_hist_segment_classifier(arr,chrom,ctrls,error,nbins)])\n",
    "        labels *= nseg\n",
    "        arr_labels += labels\n",
    "        count += nseg\n",
    "    arr_labels/=count\n",
    "    print(\"total number of segments:\",count)\n",
    "    print(\"proportion of outliers in arr:\",len([x for x in arr_labels if x<1-alpha])/len(arr_labels))\n",
    "\n",
    "def tuning(ctrls,a0=0.05,a1=0.25):\n",
    "    alpha = a0; alphas = []; spes = []\n",
    "    while alpha < a1:\n",
    "        count = 0\n",
    "        for i in range(ctrls.shape[1]):\n",
    "            count += ensemble_classifier_sample(ctrls[:,i].reshape(-1,1),ctrls[:,[j for j in range(ctrls.shape[1]) if j!=i]],nbins=8,alpha=alpha)\n",
    "        spe = 1-count/ctrls.shape[1]\n",
    "        alphas.append(alpha); spes.append(spe)\n",
    "        alpha += 0.01\n",
    "    plt.scatter(alphas,spes)\n",
    "\n",
    "def autotuning(ctrls,specificity=0.9,a0=0.05,a1=0.25):\n",
    "    alpha = a0; alphas = []; spes = []\n",
    "    while alpha < a1:\n",
    "        count = 0\n",
    "        for i in range(ctrls.shape[1]):\n",
    "            count += ensemble_classifier_sample(ctrls[:,i].reshape(-1,1),ctrls[:,[j for j in range(ctrls.shape[1]) if j!=i]],nbins=8,alpha=alpha)\n",
    "        spe = 1-count/ctrls.shape[1]\n",
    "        alphas.append(alpha); spes.append(spe)\n",
    "        if spe>=specificity:\n",
    "            print('alpha:',alpha)\n",
    "            print('leave-one-out-specificity:',spe)\n",
    "            return alpha\n",
    "        else:\n",
    "            alpha += 0.01\n",
    "    best_spe = max(spes)\n",
    "    for i,a in enumerate(alphas):\n",
    "        if spes[i]==best_spe:\n",
    "            print('alpha:',a)\n",
    "            print('leave-one-out-specificity:',best_spe)\n",
    "            return a\n",
    "\n",
    "def outlier_detection(arr,ctrls,specificity=0.9,condition_matched=True,a0=None,a1=None):\n",
    "    if condition_matched:\n",
    "        if a0 is None or a1 is None:\n",
    "            a0 = 0.05; a1=0.5\n",
    "    else: # if condition unmatched, consider larger alpha (0.1,0.15,0.2) and/or larger leave-one-out specificity to tolerate batch effects\n",
    "        if a0 is None or a1 is None:\n",
    "            a0 = 0.15; a1 = 0.5\n",
    "    alpha = autotuning(ctrls,specificity=specificity,a0=a0,a1=a1)\n",
    "    ensemble_classifier(arr,ctrls,alpha=alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "143110b2-c01c-4e30-a302-981347be43c6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha: 0.15\n",
      "leave-one-out-specificity: 0.9090909090909091\n",
      "total number of segments: 679\n",
      "proportion of outliers in arr: 0.970954356846473\n"
     ]
    }
   ],
   "source": [
    "outlier_detection(rega7080_ctdna,rega7080_ctrl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "0410e30b-476a-49df-b891-cf90140d99d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha: 0.1\n",
      "leave-one-out-specificity: 0.9560439560439561\n",
      "total number of segments: 679\n",
      "proportion of outliers in arr: 0.6846473029045643\n",
      "alpha: 0.1\n",
      "leave-one-out-specificity: 0.9560439560439561\n",
      "total number of segments: 679\n",
      "proportion of outliers in arr: 0.045454545454545456\n"
     ]
    }
   ],
   "source": [
    "outlier_detection(rega7080_ctdna,hstack((rega5343_ctrl,rctrl29s,rctrl38s)),specificity=0.95,condition_matched=False,a0=0.1,a1=0.5)\n",
    "outlier_detection(rega7080_ctrl,hstack((rega5343_ctrl,rctrl29s,rctrl38s)),specificity=0.95,condition_matched=False,a0=0.1,a1=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "e89e78fc-957a-429a-8469-dab85dbec64c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha: 0.05\n",
      "leave-one-out-specificity: 0.9210526315789473\n",
      "total number of segments: 679\n",
      "proportion of outliers in arr: 0.4444444444444444\n"
     ]
    }
   ],
   "source": [
    "outlier_detection(rctdna36s,rctrl38s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "d25b7bf8-ee13-442e-83bd-a5522a9ab014",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha: 0.15\n",
      "leave-one-out-specificity: 0.96\n",
      "total number of segments: 679\n",
      "proportion of outliers in arr: 0.25\n",
      "alpha: 0.15\n",
      "leave-one-out-specificity: 0.96\n",
      "total number of segments: 679\n",
      "proportion of outliers in arr: 0.2894736842105263\n"
     ]
    }
   ],
   "source": [
    "outlier_detection(rctdna36s,hstack((rega5343_ctrl,rctrl29s,rega7080_ctrl)),specificity=0.95,condition_matched=False)\n",
    "outlier_detection(rctrl38s,hstack((rega5343_ctrl,rctrl29s,rega7080_ctrl)),specificity=0.95,condition_matched=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "efc02f3b-315e-4c04-a299-e73feb6af01d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha: 0.47000000000000025\n",
      "leave-one-out-specificity: 0.9166666666666666\n",
      "total number of segments: 679\n",
      "proportion of outliers in arr: 1.0\n"
     ]
    }
   ],
   "source": [
    "outlier_detection(rega5343_prostate,rega5343_ctrl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "97793810-ce2e-4136-ab31-702dc42d1392",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha: 0.15\n",
      "leave-one-out-specificity: 0.9662921348314607\n",
      "total number of segments: 679\n",
      "proportion of outliers in arr: 1.0\n",
      "alpha: 0.15\n",
      "leave-one-out-specificity: 0.9662921348314607\n",
      "total number of segments: 679\n",
      "proportion of outliers in arr: 0.0\n"
     ]
    }
   ],
   "source": [
    "outlier_detection(rega5343_prostate,hstack((rctrl38s,rctrl29s,rega7080_ctrl)),specificity=0.95,condition_matched=False)\n",
    "outlier_detection(rega5343_ctrl,hstack((rctrl38s,rctrl29s,rega7080_ctrl)),specificity=0.95,condition_matched=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "1fb2014d-83a9-47ce-bd8c-94bbcfa334cd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha: 0.15\n",
      "leave-one-out-specificity: 0.9642857142857143\n",
      "total number of segments: 679\n",
      "proportion of outliers in arr: 0.0\n"
     ]
    }
   ],
   "source": [
    "outlier_detection(rctrl29s,hstack((rctrl38s,rega5343_ctrl,rega7080_ctrl)),specificity=0.95,condition_matched=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "2fbf3d7a-8c9f-4dde-9a10-50fab753f13f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha: 0.15\n",
      "leave-one-out-specificity: 0.9734513274336283\n",
      "total number of segments: 679\n",
      "proportion of outliers in arr: 1.0\n"
     ]
    }
   ],
   "source": [
    "outlier_detection(rbrcas,rctrls,specificity=0.95,condition_matched=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "55923e18-d1b5-4897-b233-26c0370e4714",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha: 0.15\n",
      "leave-one-out-specificity: 0.9734513274336283\n",
      "total number of segments: 679\n",
      "proportion of outliers in arr: 1.0\n"
     ]
    }
   ],
   "source": [
    "outlier_detection(rcrcs,rctrls,specificity=0.95,condition_matched=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "51d33ad5-3193-4ee0-87c5-457721b5af21",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha: 0.15\n",
      "leave-one-out-specificity: 0.9734513274336283\n",
      "total number of segments: 679\n",
      "proportion of outliers in arr: 1.0\n"
     ]
    }
   ],
   "source": [
    "outlier_detection(rega8460,rctrls,specificity=0.95,condition_matched=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c47fe13f-e07f-462d-b142-f444f6347349",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "id": "fe3da8d3-d749-49a2-a407-cdce4a7bea97",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1kGP contamination test\n",
    "df = pd.read_csv('3202_sequence_index.tsv',sep='\\t')\n",
    "runs = df['RUN_NAME'].tolist()\n",
    "files = df['LIBRARY_NAME'].tolist()\n",
    "sorted_files = sorted(files)\n",
    "batches = [x.split('-')[-2] for x in runs]\n",
    "bs = unique(batches)\n",
    "BATCHES = {}\n",
    "for batch in bs:\n",
    "    idxs = [i for i in range(len(batches)) if batches[i]==batch]\n",
    "    samples = array(files)[idxs]\n",
    "    BATCHES[batch] = list(samples)\n",
    "# ten unselected batches, batch1 for training, batch2-9 for testing \n",
    "batch1 = rthgs[:,[i for i in range(len(batches)) if batches[i]==bs[1]]]\n",
    "batch2 = rthgs[:,[i for i in range(len(batches)) if batches[i]==bs[2]]]\n",
    "batch3 = rthgs[:,[i for i in range(len(batches)) if batches[i]==bs[3]]]\n",
    "batch4 = rthgs[:,[i for i in range(len(batches)) if batches[i]==bs[4]]]\n",
    "batch5 = rthgs[:,[i for i in range(len(batches)) if batches[i]==bs[5]]]\n",
    "batch6 = rthgs[:,[i for i in range(len(batches)) if batches[i]==bs[6]]]\n",
    "batch7 = rthgs[:,[i for i in range(len(batches)) if batches[i]==bs[7]]]\n",
    "batch8 = rthgs[:,[i for i in range(len(batches)) if batches[i]==bs[8]]]\n",
    "batch9 = rthgs[:,[i for i in range(len(batches)) if batches[i]==bs[9]]]\n",
    "batch10 = rthgs[:,[i for i in range(len(batches)) if batches[i]==bs[10]]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "6e6f210d-7906-4c2b-b77f-5b2f6dea2151",
   "metadata": {},
   "outputs": [],
   "source": [
    "def contamination_test(chrom,test_set=hstack((batch2,batch3,batch4,batch5,batch6,batch7,batch8,batch9,batch10)),\n",
    "                       train_set=batch1,error=1e-32,bps=bkps,splits=chr_splits,nbins=8):\n",
    "    sorted_idx = sort_arr_by_chr_mean(test_set,chrom)\n",
    "    sorted_ctrl_idx = sort_arr_by_chr_mean(train_set,chrom)\n",
    "    sarr = get_segment_mean(test_set[:,sorted_idx],breaks=bps[chrom])\n",
    "    sctrl = get_segment_mean(train_set[:,sorted_ctrl_idx],breaks=bps[chrom])\n",
    "    harr = get_rank_hist(sarr,nbins)#;print(mean(harr,axis=1).shape)\n",
    "    hctrl = get_rank_hist(sctrl,nbins)\n",
    "    ctrl_labels = test_gaussian(hctrl,hctrl,error=error)\n",
    "    arr_labels = test_gaussian(harr,hctrl,error=error)\n",
    "    print(\"contamination rate for \"+chrom+\":\",len([x for x in arr_labels if x==-1])/len(arr_labels))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "7556b8a0-0b2d-488e-936d-74c951158086",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "contamination rate for chr1: 0.008583690987124463\n",
      "contamination rate for chr2: 0.008583690987124463\n",
      "contamination rate for chr3: 0.0\n",
      "contamination rate for chr4: 0.004291845493562232\n",
      "contamination rate for chr5: 0.017167381974248927\n",
      "contamination rate for chr6: 0.0\n",
      "contamination rate for chr7: 0.0\n",
      "contamination rate for chr8: 0.07725321888412018\n",
      "contamination rate for chr9: 0.27467811158798283\n",
      "contamination rate for chr10: 0.2575107296137339\n",
      "contamination rate for chr11: 0.34763948497854075\n",
      "contamination rate for chr12: 0.02145922746781116\n",
      "contamination rate for chr13: 0.012875536480686695\n",
      "contamination rate for chr14: 0.21888412017167383\n",
      "contamination rate for chr15: 0.017167381974248927\n",
      "contamination rate for chr16: 0.004291845493562232\n",
      "contamination rate for chr17: 0.0\n",
      "contamination rate for chr18: 0.03862660944206009\n",
      "contamination rate for chr19: 0.012875536480686695\n",
      "contamination rate for chr20: 0.0\n",
      "contamination rate for chr21: 0.08583690987124463\n",
      "contamination rate for chr22: 0.0\n"
     ]
    }
   ],
   "source": [
    "for ch in range(22):\n",
    "    chrom = 'chr'+str(ch+1)\n",
    "    contamination_test(chrom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "4c8ee147-02b6-45cc-928f-6a5438ecea86",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26 26 26 27 26 26 25 25 26 26\n"
     ]
    }
   ],
   "source": [
    "print(batch1.shape[1],batch2.shape[1],batch3.shape[1],batch4.shape[1],batch5.shape[1],\n",
    "      batch6.shape[1],batch7.shape[1],batch8.shape[1],batch9.shape[1],batch10.shape[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5ea7367-899b-4342-b1db-7325275fe178",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}