open_file.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"3\"\n",
    "\n",
    "GENOMES = { \"mm10\" : \"/users/kcochran/genomes/mm10_no_alt_analysis_set_ENCODE.fasta\",\n",
    "            \"hg38\" : \"/users/kcochran/genomes/GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta\" }\n",
    "\n",
    "ROOT = \"/users/kcochran/projects/cs197_cross_species_domain_adaptation/\"\n",
    "DATA_DIR = ROOT + \"data/\"\n",
    "\n",
    "SPECIES = [\"mouse\", \"human\"]\n",
    "\n",
    "TFS = [\"CTCF\", \"CEBPA\", \"HNF4A\", \"RXRA\"]\n",
    "\n",
    "TEST_DIR = \"/users/tkanell/team_covariates/data/human/accessibility/peaks.bed\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import gzip\n",
    "from collections import defaultdict\n",
    "import random\n",
    "import numpy as np\n",
    "from pyfaidx import Fasta\n",
    "from torch.utils.data import Dataset\n",
    "#from torch.utils.data import Dataset\n",
    "import pyBigWig\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Generator(Dataset):\n",
    "    letter_dict = {\n",
    "        'a':[1,0,0,0],'c':[0,1,0,0],'g':[0,0,1,0],'t':[0,0,0,1],\n",
    "        'n':[0,0,0,0],'A':[1,0,0,0],'C':[0,1,0,0],'G':[0,0,1,0],\n",
    "        'T':[0,0,0,1],'N':[0,0,0,0]}\n",
    "\n",
    "    def __init__(self, species, tf, train_val_test,\n",
    "                 seq_len = 2114, profile_len = 1000, return_labels = True):\n",
    "        \n",
    "        #assert train_val_test in [\"train\", \"val\", \"test\"]\n",
    "        ## note: kelly will give you these files, but they are basically the same as normal peak files\n",
    "        self.peakfile = TEST_DIR\n",
    "        if train_val_test == \"train\":\n",
    "            self.peakfile = PEAKS_DIR + species + \"/\" + tf + \"/filtered_peaks_chr3toY.bed\"\n",
    "        elif train_val_test == \"val\":\n",
    "            self.peakfile = PEAKS_DIR + species + \"/\" + tf + \"/filtered_peaks_chr1.bed\"\n",
    "        # else:\n",
    "        #     self.peakfile = PEAKS_DIR + species + \"/\" + tf + \"/filtered_peaks_chr2.bed\"\n",
    "            \n",
    "        #self.pos_bw = BIGWIGS_DIR + species + \"/\" + tf + \"/final.pos.bigWig\"\n",
    "        #self.neg_bw = BIGWIGS_DIR + species + \"/\" + tf + \"/final.neg.bigWig\"\n",
    "        self.prof_len = profile_len\n",
    "        self.max_jitter = 0\n",
    "        self.return_labels = return_labels\n",
    "        \n",
    "        #self.genome_file = GENOMES[species]\n",
    "        self.seq_len = seq_len\n",
    "\n",
    "        self.set_len()\n",
    "        self.coords = self.get_coords()\n",
    "        self.seqs_onehot = self.convert(self.coords)\n",
    "        self.profiles, self.logcounts = self.get_profiles_and_logcounts(self.coords)\n",
    "\n",
    "\n",
    "    def __len__(self):\n",
    "        return self.len\n",
    "    \n",
    "    \n",
    "    def set_len(self):\n",
    "        with open(self.peakfile) as f:\n",
    "            self.len = sum([1 for _ in f])\n",
    "\n",
    "\n",
    "    def get_coords(self):\n",
    "        with open(self.peakfile) as posf:\n",
    "            coords_tmp = [line.split()[:3] for line in posf]  # expecting bed file format\n",
    "        \n",
    "        coords = []\n",
    "        for coord in coords_tmp:\n",
    "            chrom, start, end = coord[0], int(coord[1]), int(coord[2])\n",
    "            window_start, window_end = expand_window(start, end,\n",
    "                                                     self.seq_len + 2 * self.max_jitter)\n",
    "            coords.append((coord[0], window_start, window_end))  # no strand consideration\n",
    "        return coords\n",
    "            \n",
    "\n",
    "    def get_profiles_and_logcounts(self, coords):\n",
    "        profiles = []\n",
    "        logcounts = []\n",
    "\n",
    "        with pyBigWig.open(self.pos_bw) as pos_bw_reader:\n",
    "            with pyBigWig.open(self.neg_bw) as neg_bw_reader:\n",
    "                for chrom, start, end in coords:\n",
    "                    # need to trim the profile length to match model output size\n",
    "                    # this is smaller than the input size bc of the receptive field\n",
    "                    # and deconv layer kernel width\n",
    "                    prof_start, prof_end = expand_window(start, end,\n",
    "                                                 self.prof_len + 2 * self.max_jitter)\n",
    "                    \n",
    "                    pos_profile = np.array(pos_bw_reader.values(chrom, prof_start, prof_end))\n",
    "                    pos_profile[np.isnan(pos_profile)] = 0\n",
    "                    neg_profile = np.array(neg_bw_reader.values(chrom, prof_start, prof_end))\n",
    "                    neg_profile[np.isnan(neg_profile)] = 0\n",
    "                    profile = np.array([pos_profile, neg_profile])\n",
    "                    \n",
    "                    pos_logcount = np.log(np.sum(pos_profile) + 1)\n",
    "                    neg_logcount = np.log(np.sum(neg_profile) + 1)\n",
    "                    logcount = np.array([pos_logcount, neg_logcount])\n",
    "\n",
    "                    profiles.append(profile)\n",
    "                    logcounts.append(logcount)\n",
    "                    \n",
    "        profiles = np.array(profiles)\n",
    "        logcounts = np.array(logcounts)\n",
    "        return profiles, logcounts\n",
    "                \n",
    "\n",
    "    def convert(self, coords):\n",
    "        seqs_onehot = []\n",
    "        with Fasta(self.genome_file) as converter:\n",
    "            for coord in coords:\n",
    "                chrom, start, stop = coord\n",
    "                assert chrom in converter\n",
    "                seq = converter[chrom][start:stop].seq\n",
    "                seq_onehot = np.array([self.letter_dict.get(x,[0,0,0,0]) for x in seq])\n",
    "                seqs_onehot.append(seq_onehot)\n",
    "\n",
    "        seqs_onehot = np.array(seqs_onehot)\n",
    "        return seqs_onehot\n",
    "\n",
    "\n",
    "    def __getitem__(self, batch_index):\t\n",
    "        # get coordinates\n",
    "        onehot = self.seqs_onehot[batch_index]\n",
    "        assert onehot.shape[0] > 0, onehot.shape\n",
    "\n",
    "        onehot = torch.tensor(onehot, dtype=torch.float).permute(1, 0)\n",
    "        \n",
    "        if not self.return_labels:\n",
    "            return onehot\n",
    "        else:\n",
    "            # get profiles and logcounts for the two strands\n",
    "            profiles = self.profiles[batch_index]\n",
    "            logcounts = self.logcounts[batch_index]\n",
    "\n",
    "            profiles = torch.tensor(profiles, dtype=torch.float)\n",
    "            logcounts = torch.tensor(logcounts, dtype=torch.float)\n",
    "            return onehot, profiles, logcounts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_profiles_and_logcounts():\n",
    "    profiles = []\n",
    "    logcounts = []\n",
    "\n",
    "    with pyBigWig.open(self.pos_bw) as pos_bw_reader:\n",
    "        with pyBigWig.open(self.neg_bw) as neg_bw_reader:\n",
    "            for chrom, start, end in coords:\n",
    "                # need to trim the profile length to match model output size\n",
    "                # this is smaller than the input size bc of the receptive field\n",
    "                # and deconv layer kernel width\n",
    "                prof_start, prof_end = expand_window(start, end,\n",
    "                                             self.prof_len + 2 * self.max_jitter)\n",
    "\n",
    "                pos_profile = np.array(pos_bw_reader.values(chrom, prof_start, prof_end))\n",
    "                pos_profile[np.isnan(pos_profile)] = 0\n",
    "                neg_profile = np.array(neg_bw_reader.values(chrom, prof_start, prof_end))\n",
    "                neg_profile[np.isnan(neg_profile)] = 0\n",
    "                profile = np.array([pos_profile, neg_profile])\n",
    "\n",
    "                pos_logcount = np.log(np.sum(pos_profile) + 1)\n",
    "                neg_logcount = np.log(np.sum(neg_profile) + 1)\n",
    "                logcount = np.array([pos_logcount, neg_logcount])\n",
    "\n",
    "                profiles.append(profile)\n",
    "                logcounts.append(logcount)\n",
    "\n",
    "    profiles = np.array(profiles)\n",
    "    logcounts = np.array(logcounts)\n",
    "    return profiles, logcounts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This working\n"
     ]
    },
    {
     "ename": "UnicodeDecodeError",
     "evalue": "'utf-8' codec can't decode byte 0x8b in position 1: invalid start byte",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mUnicodeDecodeError\u001b[0m                        Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_35211/3301993907.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"This working\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mtrain_gen\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mGenerator\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"human\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"CTCF\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"va\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"second\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;31m#source = Generator(train_species, tf, \"val\")\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/tmp/ipykernel_35211/1675135143.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, species, tf, train_val_test, seq_len, profile_len, return_labels)\u001b[0m\n\u001b[1;32m     27\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mseq_len\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mseq_len\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     28\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 29\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_len\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     30\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcoords\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_coords\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     31\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mseqs_onehot\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconvert\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcoords\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/tmp/ipykernel_35211/1675135143.py\u001b[0m in \u001b[0;36mset_len\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m     39\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mset_len\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m         \u001b[0;32mwith\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpeakfile\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 41\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlen\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0m_\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     42\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     43\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/tmp/ipykernel_35211/1675135143.py\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m     39\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mset_len\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m         \u001b[0;32mwith\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpeakfile\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 41\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlen\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0m_\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     42\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     43\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/197/lib/python3.7/codecs.py\u001b[0m in \u001b[0;36mdecode\u001b[0;34m(self, input, final)\u001b[0m\n\u001b[1;32m    320\u001b[0m         \u001b[0;31m# decode input (taking the buffer into account)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    321\u001b[0m         \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuffer\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0minput\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 322\u001b[0;31m         \u001b[0;34m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconsumed\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_buffer_decode\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0merrors\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfinal\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    323\u001b[0m         \u001b[0;31m# keep undecoded input until the next call\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    324\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuffer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mconsumed\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mUnicodeDecodeError\u001b[0m: 'utf-8' codec can't decode byte 0x8b in position 1: invalid start byte"
     ]
    }
   ],
   "source": [
    "print(\"This working\")\n",
    "train_gen = Generator(\"human\", \"CTCF\", \"va\")\n",
    "print(\"second\")\n",
    "#source = Generator(train_species, tf, \"val\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mykernel",
   "language": "python",
   "name": "mykernel"
  },
  "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}