4_deep_learning_for_sequences.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 特定の転写因子が結合するか否かの予測"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "転写因子の結合サイトに関連した実験が、ChIP-seq等を用いて、多数実施されるようになってきました。ここでは転写因子CTCFを対象として実施したChIP-seqの実験結果を用いて、特定の配列を与えた際に、CTCFが結合するか否かを予測するクラス分類問題を考えます。\n",
    "\n",
    "ここでは、CTCFの結合の有無を機械学習的に調べるために、以下の様なセッティングを考えます。\n",
    "\n",
    "* クラス分類問題を考える。つまり、何かをクラス0、何かをクラス1として、0と1を分類する問題とする。\n",
    "* クラス1はCTCFが結合する配列である。学習データは、ChIP-seqのCTCFに関する実験で、ピークの中心周辺の配列で作成する。\n",
    "    * 簡単のためピーク前後5０bp(合計100塩基）とします。\n",
    "* クラス0はCTCFが結合しない配列であり、同一の実験で、ピークが無い領域の配列とする\n",
    "    * 長さは揃えて１００塩基とします。\n",
    "* これらのクラスを分ける分類器が作成できたら、任意の配列を入力として、結合するかしないかの判定をする機械の作成ができる。\n",
    "\n",
    "<img src=\"img/DL_chipseq.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "今回、ChIP-seqのデータからピークコールを行った後のデータが格納されている、[ChIP-Atlas](https://chip-atlas.org/) のデータを利用します。CTCFのChIP-seqを行った実験として、19人分のCTCFの結合サイトを解析した[Kasowski et al.](http://science.sciencemag.org/content/342/6159/750.long)の中の一つ、SRX356455 ( [ChIP-AtlasのSRX356455のページ](https://chip-atlas.org/view?id=SRX356455))　を利用します。ピークコール後の結果は、ChIP-AtlasのDownloadsから、ダウンロードできます。\n",
    "\n",
    "ここでは、しきい値 p<1E-0.5 の結果から、ピークの前後１００塩基、および、ピークの無い位置からランダムに100塩基を取り出した前処理済みのデータをdata/sequences 以下に用意したので、これを用います。```SRX356455.05_peak_seq_100.txt```がピーク前後、```SRX356455.05_random_seq_100.txt```が、それ以外の領域から取ったものです。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ライブラリの読み込み"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "学習は、配列の読み込みとネットワーク構造を除けば、特徴量からの学習と同様です。深層学習で利用するライブラリを読み込みます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.optim as optim\n",
    "from torch.optim import lr_scheduler\n",
    "from torchvision import models, transforms\n",
    "import torch.utils.data as data\n",
    "import torch.nn.functional as F\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "import os\n",
    "import copy\n",
    "import argparse\n",
    "import time\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# データの読み込み"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2節で酵母画像の特徴量から学習をした際と、同じ流れです。\n",
    "\n",
    "配列解析に向けた大きな変更点は、塩基（ATGC）をそのまま学習に用いた場合、計算機にはただのアルファベットと思われてしまうので、これをベクトル情報に変換します。その方法は、酵母の芽の大きさの分類（no, medium,...） を0次元目、1次元目・・・と割当てた場合と同様で、Aを0次元目、Tを1次元目、としたベクトルを作成します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "SEQ_LENGTH = 100\n",
    "\n",
    "\n",
    "def make_dataset(datadir):\n",
    "    pos_seq = \"SRX356455.05_peak_seq_100.txt\"\n",
    "    neg_seq = \"SRX356455.05_random_seq_100.txt\"\n",
    "    # id      chr     start   end     seq\n",
    "    data = pd.read_csv(os.path.join(datadir, \"sequences\", pos_seq), sep=\"\\t\")\n",
    "    sequences = []\n",
    "    classes = []\n",
    "    for index, row in data[[\"id\", \"seq\"]].iterrows():\n",
    "        y = 1\n",
    "        seq_vector = seq2vector(row[\"seq\"])\n",
    "        if len(seq_vector) == 0:\n",
    "            continue\n",
    "        sequences.append(seq2vector(row[\"seq\"]))\n",
    "        classes.append(np.array(y))\n",
    "    data = pd.read_csv(os.path.join(datadir, \"sequences\", neg_seq), sep=\"\\t\")\n",
    "    for index, row in data[[\"id\", \"seq\"]].iterrows():\n",
    "        y = 0\n",
    "        seq_vector = seq2vector(row[\"seq\"])\n",
    "        if len(seq_vector) == 0:\n",
    "            continue\n",
    "        sequences.append(seq2vector(row[\"seq\"]))\n",
    "        classes.append(np.array(y))\n",
    "    return sequences, classes\n",
    "\n",
    "\n",
    "def seq2vector(seq):\n",
    "    if type(seq) is not str: # Case on Null sequence\n",
    "        return np.zeros((0,0))\n",
    "    seq_array = np.zeros((4, SEQ_LENGTH))\n",
    "    flag = 0\n",
    "    for i in range(SEQ_LENGTH):\n",
    "        s = seq[i]\n",
    "        if s == \"a\" or s == \"A\":\n",
    "            seq_array[0, i] = 1\n",
    "        elif s == \"c\" or s == \"C\":\n",
    "            seq_array[1, i] = 1\n",
    "        elif s == \"g\" or s == \"G\":\n",
    "            seq_array[2, i] = 1\n",
    "        elif s == \"t\" or s == \"T\":\n",
    "            seq_array[3, i] = 1\n",
    "        else:\n",
    "            flag += 1\n",
    "    if len(seq) == flag: # Case on N sequence\n",
    "        return np.zeros((0,0))\n",
    "    seq_array = seq_array.astype(np.float32)\n",
    "    return seq_array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "データの読み込みも、特徴量からの学習時と同様です。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "datadir = \"data\"\n",
    "\n",
    "class DatasetFolder(data.Dataset):\n",
    "    def __init__(self, X, y):\n",
    "        self.samples = X\n",
    "        self.targets = y\n",
    "        self.transforms = transforms.Compose([\n",
    "            ToTensorOfTarget()\n",
    "        ])\n",
    "\n",
    "    def __getitem__(self, index):\n",
    "        sample = self.samples[index]\n",
    "        sample = self.transforms(sample)\n",
    "        target = self.targets[index]\n",
    "        target = self.transforms(target)\n",
    "        return sample, target\n",
    "\n",
    "    def __len__(self):\n",
    "        return len(self.samples)\n",
    "\n",
    "\n",
    "class ToTensorOfTarget(object):\n",
    "    def __call__(self, target):\n",
    "        return torch.from_numpy(target)\n",
    "\n",
    "# 全体を、training, valid, testに分ける。ここでは、3:1:1 に分割。\n",
    "# training + valid が、機械学習の training data 相当。\n",
    "X, y = make_dataset(datadir)\n",
    "X_tmp, X_test, y_tmp, y_test = train_test_split(\n",
    "    X, y, test_size = 0.20)\n",
    "X_train, X_val, y_train, y_val = train_test_split(\n",
    "    X_tmp, y_tmp, test_size = 0.25\n",
    ")\n",
    "\n",
    "\n",
    "sequence_datasets = {\n",
    "    'train':DatasetFolder(X_train, y_train),\n",
    "    'val':DatasetFolder(X_val, y_val),\n",
    "    'test': DatasetFolder(X_test, y_test)\n",
    "}\n",
    "\n",
    "dataset_sizes = {x: len(sequence_datasets[x]) for x in ['train', 'val', 'test']}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ネットワークモデルの定義"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "配列解析では、画像解析の様なゴールドスタンダードのモデルが存在する訳ではありません。特徴量抽出で行ったようなDense層を重ねることも一つですし、ローカルな特徴量を取るために画像同様Convoludionを利用する方法、更には、配列を時系列の音声処理の様にみなして解析を実施していく方法(WaveNet)も考えられます。\n",
    "\n",
    "ここでは、前節の画像と同様のConvolutionを用いた方法で、シンプルにネットワークを構築します。\n",
    "ネットワーク構成も、Conv -> Pooling -> Conv -> Pooling -> Dense -> Output という前節同様の構成を利用しています。\n",
    "変更点は、画像が２次元情報だったのに対して、配列は1次元情報なので、1次元方法にのみConvolutionとPoolingを実施している点です。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Net(nn.Module):\n",
    "    def __init__(self):\n",
    "        super(Net, self).__init__()\n",
    "        self.conv1 = nn.Conv1d(4, 8, 3, padding=1) # 100\n",
    "        self.pool = nn.MaxPool1d(2) # 50\n",
    "        self.conv2 = nn.Conv1d(8, 16, 3, padding=1) # 50\n",
    "        self.fc1 = nn.Linear(16 * 25 * 1, 20) # channel_num * x * y\n",
    "        self.fc2 = nn.Linear(20, 2)\n",
    "\n",
    "    def forward(self, x):\n",
    "        x = self.pool(F.relu(self.conv1(x))) # 100 -> 100 -> 50\n",
    "        x = self.pool(F.relu(self.conv2(x))) # 50 -> 50 -> 25\n",
    "        x = x.view(-1, 16 * 25 * 1) # channel_num * x * y\n",
    "        x = F.relu(self.fc1(x))\n",
    "        x = self.fc2(x)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "次に、学習ステップの詳細を定義します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_test_accuracy(model, criterion, optimizer, phase):\n",
    "    running_loss = 0.0\n",
    "    running_corrects = 0\n",
    "    model.train(False)\n",
    "\n",
    "    for inputs, labels in dataloaders[phase]:\n",
    "        inputs = inputs.to(device)\n",
    "        labels = labels.to(device)\n",
    "        \n",
    "        # 訓練のときだけ履歴を保持する\n",
    "        with torch.set_grad_enabled(phase == 'train'):\n",
    "            outputs = model(inputs)\n",
    "            #_, classnums = torch.max(labels, 1)\n",
    "            _, preds = torch.max(outputs, 1)\n",
    "            loss = criterion(outputs, labels)\n",
    "\n",
    "        # 統計情報\n",
    "        running_loss += loss.item() * inputs.size(0)\n",
    "        running_corrects += torch.sum(preds == labels)\n",
    "\n",
    "    # サンプル数で割って平均を求める\n",
    "    epoch_loss = running_loss / dataset_sizes[phase]\n",
    "    epoch_acc = running_corrects.double() / dataset_sizes[phase]\n",
    "    print('On Test:\\tLoss: {:.4f} Acc: {:.4f}'.format(epoch_loss, epoch_acc))\n",
    "\n",
    "\n",
    "def train_model(model, criterion, optimizer, scheduler, num_epochs=25):\n",
    "    since = time.time()\n",
    "    # 途中経過でモデル保存するための初期化\n",
    "    best_model_wts = copy.deepcopy(model.state_dict())\n",
    "    best_acc = 0.0\n",
    "    # 時間計測用\n",
    "    end = time.time()\n",
    "\n",
    "    print(model)\n",
    "    print()\n",
    "\n",
    "    for epoch in range(num_epochs):\n",
    "        print('Epoch:{}/{}'.format(epoch, num_epochs - 1), end=\"\")\n",
    "\n",
    "        # 各エポックで訓練+バリデーションを実行\n",
    "        for phase in ['train', 'val']:\n",
    "            if phase == 'train':\n",
    "                scheduler.step()\n",
    "                model.train(True)  # training mode\n",
    "            else:\n",
    "                model.train(False)  # evaluate mode\n",
    "\n",
    "            running_loss = 0.0\n",
    "            running_corrects = 0\n",
    "\n",
    "            for inputs, labels in dataloaders[phase]:\n",
    "                inputs = inputs.to(device)\n",
    "                labels = labels.to(device)\n",
    "\n",
    "                optimizer.zero_grad()\n",
    "\n",
    "                # 訓練のときだけ履歴を保持する\n",
    "                with torch.set_grad_enabled(phase == 'train'):\n",
    "                    outputs = model(inputs)\n",
    "                    #_, classnums = torch.max(labels, 1)\n",
    "                    _, preds = torch.max(outputs, 1)\n",
    "                    loss = criterion(outputs, labels)\n",
    "                    # backward + optimize only if in training phase\n",
    "                    if phase == 'train':\n",
    "                        loss.backward()\n",
    "                        optimizer.step()\n",
    "\n",
    "                # 統計情報\n",
    "                running_loss += loss.item() * inputs.size(0)\n",
    "                running_corrects += torch.sum(preds == labels)\n",
    "\n",
    "            # サンプル数で割って平均を求める\n",
    "            epoch_loss = running_loss / dataset_sizes[phase]\n",
    "            epoch_acc = running_corrects.double() / dataset_sizes[phase]\n",
    "\n",
    "            print('\\t{} Loss: {:.4f} Acc: {:.4f} Time: {:.4f}'.format(phase, epoch_loss, epoch_acc, time.time()-end), end=\"\")\n",
    "\n",
    "            # 精度が改善したらモデルを保存する\n",
    "            if phase == 'val' and epoch_acc > best_acc:\n",
    "                best_acc = epoch_acc\n",
    "                best_model_wts = copy.deepcopy(model.state_dict())\n",
    "            end = time.time()\n",
    "\n",
    "        print()\n",
    "\n",
    "    time_elapsed = time.time() - since\n",
    "    print()\n",
    "    print('Training complete in {:.0f}m {:.0f}s'.format(time_elapsed // 60, time_elapsed % 60))\n",
    "    print('Best val acc: {:.4f}'.format(best_acc))\n",
    "\n",
    "    # load best model weights\n",
    "    model.load_state_dict(best_model_wts)\n",
    "    return model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "バッチごとの読み込み設定をしましょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# バッチサイズ分のデータを読み込む。\n",
    "# training はデータをシャッフルし、読み込み始める配列をランダムにする。\n",
    "# 他はシャッフルの必要なし。\n",
    "batch_size=64\n",
    "workers=0\n",
    "dataloaders = {\n",
    "    'train': torch.utils.data.DataLoader(\n",
    "        sequence_datasets['train'],\n",
    "        batch_size=batch_size,\n",
    "        shuffle=True,\n",
    "        num_workers=workers),\n",
    "    'val': torch.utils.data.DataLoader(\n",
    "        sequence_datasets['val'],\n",
    "        batch_size=batch_size,\n",
    "        shuffle=False,\n",
    "        num_workers=workers),\n",
    "    'test': torch.utils.data.DataLoader(\n",
    "        sequence_datasets['test'],\n",
    "        batch_size=batch_size,\n",
    "        shuffle=False,\n",
    "        num_workers=workers)\n",
    "}\n",
    "dataset_sizes = {x: len(sequence_datasets[x]) for x in ['train', 'val', 'test']}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "では、早速学習をしてみましょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Net(\n",
      "  (conv1): Conv1d(4, 8, kernel_size=(3,), stride=(1,), padding=(1,))\n",
      "  (pool): MaxPool1d(kernel_size=2, stride=2, padding=0, dilation=1, ceil_mode=False)\n",
      "  (conv2): Conv1d(8, 16, kernel_size=(3,), stride=(1,), padding=(1,))\n",
      "  (fc1): Linear(in_features=400, out_features=20, bias=True)\n",
      "  (fc2): Linear(in_features=20, out_features=2, bias=True)\n",
      ")\n",
      "\n",
      "Epoch:0/9\ttrain Loss: 0.5766 Acc: 0.6933 Time: 6.0771\tval Loss: 0.5411 Acc: 0.7212 Time: 0.8388\n",
      "Epoch:1/9\ttrain Loss: 0.5245 Acc: 0.7387 Time: 4.9355\tval Loss: 0.5114 Acc: 0.7481 Time: 0.8380\n",
      "Epoch:2/9\ttrain Loss: 0.4668 Acc: 0.7806 Time: 7.6239\tval Loss: 0.4533 Acc: 0.7926 Time: 0.8236\n",
      "Epoch:3/9\ttrain Loss: 0.4306 Acc: 0.8034 Time: 4.8213\tval Loss: 0.4166 Acc: 0.8133 Time: 0.8161\n",
      "Epoch:4/9\ttrain Loss: 0.4074 Acc: 0.8166 Time: 8.6646\tval Loss: 0.4018 Acc: 0.8243 Time: 0.8111\n",
      "Epoch:5/9\ttrain Loss: 0.3906 Acc: 0.8259 Time: 4.6591\tval Loss: 0.3900 Acc: 0.8303 Time: 0.8693\n",
      "Epoch:6/9\ttrain Loss: 0.3709 Acc: 0.8356 Time: 5.1666\tval Loss: 0.3798 Acc: 0.8345 Time: 1.8978\n",
      "Epoch:7/9\ttrain Loss: 0.3556 Acc: 0.8455 Time: 7.6154\tval Loss: 0.3700 Acc: 0.8414 Time: 0.8238\n",
      "Epoch:8/9\ttrain Loss: 0.3420 Acc: 0.8530 Time: 4.9556\tval Loss: 0.3594 Acc: 0.8439 Time: 1.6486\n",
      "Epoch:9/9\ttrain Loss: 0.3339 Acc: 0.8568 Time: 8.5369\tval Loss: 0.3598 Acc: 0.8468 Time: 0.8179\n",
      "\n",
      "Training complete in 1m 13s\n",
      "Best val acc: 0.8468\n",
      "On Test:\tLoss: 0.3540 Acc: 0.8501\n"
     ]
    }
   ],
   "source": [
    "device_name = \"cpu\"\n",
    "epochs = 10\n",
    "batch_size = 64\n",
    "lr = 0.01\n",
    "momentum = 0.9\n",
    "\n",
    "device = torch.device(device_name)\n",
    "model = Net()\n",
    "model = model.to(device_name)\n",
    "\n",
    "# 損失関数、最適化方法、学習率の更新方法を定義\n",
    "criterion = nn.CrossEntropyLoss()\n",
    "optimizer = optim.SGD(model.parameters(), lr=lr, momentum=momentum)\n",
    "exp_lr_scheduler = lr_scheduler.StepLR(optimizer, step_size=10, gamma=0.7)\n",
    "\n",
    "# 実際の学習を実施する\n",
    "# 結果出力用ファイルのprefix\n",
    "model = train_model(model, criterion, optimizer, exp_lr_scheduler, num_epochs=epochs)\n",
    "# テストデータでの精度を求める\n",
    "print_test_accuracy(model, criterion, optimizer, 'test')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ここまでで、100塩基の配列を与えれば、CTCFが結合しそうか否かを判別することができそうな感じが見えてきます。一方で、CTCFの結合サイトの配列が何であるかは、ここからは分からないことが、深層学習がブラックボックスだと言われる所以となります。\n",
    "\n",
    "機械学習を用いた配列解析は以上です。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "今回の内容をまとめたソースコードは、```src/chipseq_prediction_cnn.py``` になります。\n",
    "シェルから\n",
    "```bash\n",
    "$ python src/chipseq_prediction_cnn.py data -b 64 --lr 0.1 --epoch 20 -j 0\n",
    "```\n",
    "の様な形で実行可能です。\n",
    "\n",
    "です。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 参考：学習データの作成"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* ここから先は、上記で利用している配列ファイルの作成方法なので、実行の必要はありません。（結構時間もかかります）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ChIP-atlasからのデータのダウンロードは、以下の様に行いました。\n",
    "ChIP-atlasから、対象の実験(SRX356455)を検索します。\n",
    "\n",
    "<img src=\"img/DL_chip_atlas_search.png\">\n",
    "\n",
    "次に、Downloadから、BEDファイルをダウンロードします。\n",
    "\n",
    "<img src=\"img/DL_chip_atlas_dl.png\">\n",
    "\n",
    "ダウンロードしたBEDファイルは、```script```ディレクトリに入れておきます。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ChIP-AtlasからBED形式のピーク情報(SRX356455.05.bed など)をダウンロードして、```script```のディレクトリに入れます。\n",
    "また、ChIP-Atlasは、現在は今回はヒトのhg19を利用しているので、ヒトゲノムの配列をダウンロードします。\n",
    "```bash\n",
    "$ cd script\n",
    "$ mkdir ref_dir\n",
    "$ cd ref_dir\n",
    "$ wget -c http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz\n",
    "$ tar zxvf chromFa.tar.gz\n",
    "$ cd ..\n",
    "```\n",
    "\n",
    "まず、BED形式から、ピーク位置情報を抽出します。\n",
    "```bash\n",
    "$ python get_chipseq_peak_info.py SRX356455.05.bed 100\n",
    "```\n",
    "これで、SRX356455.05_peak_extract.txt, SRX356455.05_peak_info.txtの２つのファイルが作成されます。それぞれ\n",
    "\n",
    "* _peak_info.txt: 結合領域\n",
    "* _peak_extract.txt: ピークの中心周辺100bpの領域\n",
    "\n",
    "を表しています。\n",
    "```bash\n",
    "$ wc -l SRX356455.05_peak_extract.txt\n",
    "   54613 SRX356455.05_peak_extract.txt\n",
    "```\n",
    "5.5万箇所くらいあることがわかります。\n",
    "\n",
    "次に、ピーク以外の領域を取り出します。ここでは、5万配列を取り出します（ピークが5万箇所くらいあるので、その数に近い数を取り出します）。\n",
    "```bash\n",
    "$ python random_seq_except_listed.py ref_dir 100 50000 -l SRX356455.05_peak_info.txt -o SRX356455.05_random_100.txt\n",
    "```\n",
    "少々時間がかかります。選択された領域は、```SRX356455.05_random_100.txt```に格納されます。\n",
    "\n",
    "最後に、作成した _peak_extract と、_ranom_100.txt のファイルを\n",
    "\n",
    "```bash\n",
    "$ python extract_seq.py SRX356455.05_peak_extract.txt ref_dir --id-column name > SRX356455.05_peak_seq_100.txt\n",
    "$ python extract_seq.py SRX356455.05_random_100.txt ref_dir > SRX356455.05_random_seq_100.txt\n",
    "```\n",
    "これで、```SRX356455.05_peak_seq_100.txt```(ピーク周辺), ```SRX356455.05_random_seq.txt```(ピークの無いところ)の情報が得られます。前半は、このデータを利用しました。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}