hclu_dev.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 187,
   "id": "762ab4f1-1325-44fa-b762-7cc116dd1e8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyranges as pr\n",
    "import pyfaidx\n",
    "import re\n",
    "import textwrap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "id": "2e9b0263-3c75-4caa-9bbb-bb7d63e30da6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_chr(seq, ofile, chr_name, line_lim=None):\n",
    "    with open(ofile, 'w') as ofile:\n",
    "        ofile.write(f'>{chr_name}\\n')\n",
    "        if line_lim:\n",
    "            wrapped_seq = textwrap.wrap(seq, width=line_lim)\n",
    "            for seq in wrapped_seq:\n",
    "                ofile.write(seq+'\\n')\n",
    "        else:\n",
    "            ofile.write(str(seq))\n",
    "        \n",
    "def get_gene_seq(fa_file,\n",
    "                gtf_file, \n",
    "                gene,\n",
    "                ofile=None,\n",
    "                slack=0):\n",
    "        \n",
    "    gtf_df = pr.read_gtf(gtf_file, as_df=True)\n",
    "\n",
    "    gtf_df = gtf_df.loc[(gtf_df.gene_name==gene)&(gtf_df.Feature=='gene')]\n",
    "    assert len(gtf_df.index) == 1\n",
    "\n",
    "    # get start and end of gene\n",
    "    start = min(gtf_df['Start'].values[0], gtf_df['End'].values[0])-slack\n",
    "    end = max(gtf_df['Start'].values[0], gtf_df['End'].values[0])+slack\n",
    "    ch = gtf_df['Chromosome'].values[0]\n",
    "    print(ch)\n",
    "    print(start)\n",
    "    print(end)\n",
    "    strand = gtf_df['Strand'].values[0]\n",
    "\n",
    "    fa = pyfaidx.Fasta(fa_file)\n",
    "    if strand == '+':\n",
    "        gene_seq = fa[ch][start:end]\n",
    "    else: \n",
    "        gene_seq = fa[ch][start:end].complement\n",
    "\n",
    "    if ofile:\n",
    "        write_chr(gene_seq.seq, ofile, gene)\n",
    "    \n",
    "    return gene_seq.seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "e9d255a7-b841-4691-86d6-c3d0f31b0fd0",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr14\n",
      "65966982\n",
      "65983047\n"
     ]
    }
   ],
   "source": [
    "# mouse\n",
    "gene = 'Clu'\n",
    "gtf_file = 'ref/annot.gtf'\n",
    "fa_file = 'ref/genome.fa'\n",
    "ofile = 'ref/clu.fa'\n",
    "seq = get_gene_seq(fa_file,\n",
    "             gtf_file, \n",
    "             gene,\n",
    "             ofile,\n",
    "             slack=1500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "id": "0dc8de5c-939e-4dd3-8ed1-14676ab97e5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def replace_seq(mod_fa,\n",
    "                seq):\n",
    "    \"\"\"\n",
    "    Replace sequence in seq with sequence in mod_fa. Mod_fa must contain some overlap\n",
    "    with the sequence in seq to find the insertion sites.\n",
    "    \n",
    "    Parameters:\n",
    "        mod_fa (str): Path to fasta file with sequence to insert\n",
    "        seq (str): Sequence to add mod_fa into\n",
    "    \"\"\"\n",
    "    \n",
    "    # get anchors at 5' and 3' ends of sequence\n",
    "    fa = pyfaidx.Fasta(modified_fa)\n",
    "    mod_seq = fa[list(fa.keys())[0]][0:-1].seq\n",
    "    \n",
    "    anchor_5 = mod_seq[:20]\n",
    "    anchor_3 = mod_seq[-20:]\n",
    "    \n",
    "    start_inds = [(m.start(), m.end()) for m in re.finditer(anchor_5, seq, re.IGNORECASE)]\n",
    "    assert len(start_inds) == 1\n",
    "\n",
    "    end_inds = [(m.start(), m.end()) for m in re.finditer(anchor_3, seq, re.IGNORECASE)]\n",
    "    assert len(end_inds) == 1\n",
    "\n",
    "    assert seq[start_inds[0][0]:start_inds[0][1]].lower() == anchor_5.lower()\n",
    "    assert seq[end_inds[0][0]:end_inds[0][1]].lower() == anchor_3.lower()\n",
    "\n",
    "    # replace_seq = seq[start_inds[0][1]:end_inds[0][0]] # this will replace excluding anchors\n",
    "    replace_seq = seq[start_inds[0][0]:end_inds[0][1]] # this will replace inluding anchors\n",
    "    with open('ref/replace_seq.fa', 'w') as ofile:\n",
    "        ofile.write(replace_seq)\n",
    "        \n",
    "    assert replace_seq[:20].lower() == anchor_5.lower()\n",
    "    assert replace_seq[-20:].lower() == anchor_3.lower()\n",
    "    \n",
    "    # print(replace_seq[0])\n",
    "    # print(replace_seq[-1])\n",
    "\n",
    "    inds = [(m.start(), m.end()) for m in re.finditer(replace_seq, seq)]\n",
    "    assert len(inds) == 1  \n",
    "    \n",
    "    new_len = len(seq)-len(replace_seq)+len(mod_seq)\n",
    "    new_seq = seq.replace(replace_seq, mod_seq)\n",
    "    assert len(new_seq) == new_len\n",
    "    return new_seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "id": "2ab7f349-1bda-4f38-a999-186e9ec25d2e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr14\n",
      "65965482\n",
      "65984547\n"
     ]
    }
   ],
   "source": [
    "gene = 'Clu'\n",
    "gtf_file = 'ref/annot.gtf'\n",
    "fa_file = 'ref/genome.fa'\n",
    "seq = get_gene_seq(fa_file,\n",
    "             gtf_file, \n",
    "             gene,\n",
    "             ofile = 'ref/clu.fa',\n",
    "             slack=3000)\n",
    "\n",
    "mod_fa = 'ref/Clu-2kb_info/ID16848_consensus_040622.fas'\n",
    "ofile = 'ref/hClu.fa'\n",
    "\n",
    "new_seq = replace_seq(mod_fa,\n",
    "                      seq)\n",
    "write_chr(new_seq, ofile, 'hClu', line_lim=80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "id": "6be019b4-736f-4b15-a2fa-5d236741be1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # human\n",
    "# gene = 'CLU'\n",
    "# gtf_file = 'ref/human_annot.gtf'\n",
    "# fa_file = 'ref/human_genome.fa'\n",
    "# ofile = 'ref/human_clu.fa'\n",
    "# get_gene_seq(fa_file,\n",
    "#              gtf_file, \n",
    "#              gene,\n",
    "#              ofile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "id": "1ca36d64-8eda-4f86-a447-c0d48225e150",
   "metadata": {},
   "outputs": [],
   "source": [
    "modified_fa = 'ref/Clu-2kb_info/ID16848_consensus_040622.fas'\n",
    "fa = pyfaidx.Fasta(modified_fa)\n",
    "mod_seq = fa[list(fa.keys())[0]][0:-1].seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "9b0010fc-57e9-4eb6-af11-d326b242062c",
   "metadata": {},
   "outputs": [],
   "source": [
    "mouse_5_anchor='ATAACCCGCAGATCCAT' # hand checked if these exist in the Clu mm10\n",
    "mouse_3_anchor='TGTACAGGCACTCAA'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "id": "43a044d9-8c94-4ab4-a8b6-91a6996dbdb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "start_inds = [(m.start(), m.end()) for m in re.finditer(mouse_5_anchor, seq)]\n",
    "assert len(start_inds) == 1\n",
    "\n",
    "end_inds = [(m.start(), m.end()) for m in re.finditer(mouse_3_anchor, seq)]\n",
    "assert len(end_inds) == 1\n",
    "\n",
    "assert seq[start_inds[0][0]:start_inds[0][1]] == mouse_5_anchor\n",
    "assert seq[end_inds[0][0]:end_inds[0][1]] == mouse_3_anchor\n",
    "\n",
    "replace_seq = seq[start_inds[0][1]:end_inds[0][0]]\n",
    "# print(replace_seq[0])\n",
    "# print(replace_seq[-1])\n",
    "\n",
    "inds = [(m.start(), m.end()) for m in re.finditer(replace_seq, seq)]\n",
    "assert len(inds) == 1\n",
    "\n",
    "# seq = seq.replace(replace_seq, mod_seq)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "1415c2de-98b0-439c-8cfe-39eafb1b79b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19065"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(seq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "f304af02-581e-476c-904f-d5c78abfc27a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATAACCCGCAGATCCATA\n"
     ]
    }
   ],
   "source": [
    "# temp_ind = seq.find(mouse_5_anchor)\n",
    "# start_ind = seq.find(mouse_5_anchor)+len(mouse_5_anchor)\n",
    "# print(seq[temp_ind:start_ind+1])\n",
    "\n",
    "# end_ind = seq.find(mouse_3_anchor)\n",
    "# temp_ind = seq.find(mouse_3_anchor)+len(mouse_3_anchor)\n",
    "# # print(seq[end_ind:temp_ind])\n",
    "\n",
    "# # replace_seq = seq[start_ind+2:end_ind-1]\n",
    "# # # replace_seq = seq[start_ind:end_ind]\n",
    "\n",
    "# # print(replace_seq[0])\n",
    "# # print(replace_seq[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "de95cfa7-0131-44eb-8737-6c2f9539dccd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "end_ind"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b019ee0f-74a3-4b4c-b493-f48fa9a828db",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6caa1bc3-94cc-4a27-9205-f5b5a1c4b164",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "8bd07498-0897-45b6-b1be-789c081912e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # mouse\n",
    "\n",
    "\n",
    "# gtf_df = pr.read_gtf(gtf_file, as_df=True)\n",
    "\n",
    "# gtf_df = gtf_df.loc[(gtf_df.gene_name==gene)&(gtf_df.Feature=='gene')]\n",
    "# assert len(gtf_df.index) == 1\n",
    "\n",
    "# # get start and end of gene\n",
    "# start = min(gtf_df['Start'].values[0], gtf_df['End'].values[0])\n",
    "# end = max(gtf_df['Start'].values[0], gtf_df['End'].values[0])\n",
    "# ch = gtf_df['Chromosome'].values[0]\n",
    "\n",
    "# fa = pyfaidx.Fasta(fa_file)\n",
    "# gene_seq = fa[ch][start:end]\n",
    "\n",
    "# with open(mouse_ofile, 'w') as ofile:\n",
    "#     ofile.write('>clu\\n')\n",
    "#     ofile.write(str(gene_seq.seq))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "0e9cf7fc-9c0e-4e96-bdd6-39fadb29e9f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # human\n",
    "# gene = 'CLU'\n",
    "# gtf_file = 'ref/human_annot.gtf'\n",
    "# fa_file = 'ref/human_genome.fa'\n",
    "# human_ofile = 'ref/human_clu.fa'\n",
    "\n",
    "# gtf_df = pr.read_gtf(gtf_file, as_df=True)\n",
    "\n",
    "# gtf_df = gtf_df.loc[(gtf_df.gene_name==gene)&(gtf_df.Feature=='gene')]\n",
    "# assert len(gtf_df.index) == 1\n",
    "\n",
    "# # get start and end of gene\n",
    "# start = min(gtf_df['Start'].values[0], gtf_df['End'].values[0])\n",
    "# end = max(gtf_df['Start'].values[0], gtf_df['End'].values[0])\n",
    "# ch = gtf_df['Chromosome'].values[0]\n",
    "# strand = gtf_df['Strand'].values[0]\n",
    "\n",
    "# fa = pyfaidx.Fasta(fa_file)\n",
    "# if strand == '+':\n",
    "#     gene_seq = fa[ch][start:end]\n",
    "# else: \n",
    "#     gene_seq = fa[ch][start:end].complement\n",
    "\n",
    "# with open(human_ofile, 'w') as ofile:\n",
    "#     ofile.write('>clu\\n')\n",
    "#     ofile.write(str(gene_seq.seq))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "0522e453-137a-4685-95c0-b55c4d14f114",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "fe652a8e-ac40-41ef-8c02-5395d68e1ebb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "60b77e6a-b7a6-42a0-b085-35f32e5848cf",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "04f41a1a-1ce2-4745-a8ce-7ea530a949c1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Chromosome</th>\n",
       "      <th>Source</th>\n",
       "      <th>Feature</th>\n",
       "      <th>Start</th>\n",
       "      <th>End</th>\n",
       "      <th>Score</th>\n",
       "      <th>Strand</th>\n",
       "      <th>Frame</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_type</th>\n",
       "      <th>...</th>\n",
       "      <th>transcript_type</th>\n",
       "      <th>transcript_name</th>\n",
       "      <th>transcript_support_level</th>\n",
       "      <th>tag</th>\n",
       "      <th>havana_transcript</th>\n",
       "      <th>exon_number</th>\n",
       "      <th>exon_id</th>\n",
       "      <th>ont</th>\n",
       "      <th>protein_id</th>\n",
       "      <th>ccdsid</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1133585</th>\n",
       "      <td>chr8</td>\n",
       "      <td>HAVANA</td>\n",
       "      <td>gene</td>\n",
       "      <td>27596916</td>\n",
       "      <td>27615031</td>\n",
       "      <td>.</td>\n",
       "      <td>-</td>\n",
       "      <td>.</td>\n",
       "      <td>ENSG00000120885.21</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>ncRNA_host</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 24 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        Chromosome  Source Feature     Start       End Score Strand Frame  \\\n",
       "1133585       chr8  HAVANA    gene  27596916  27615031     .      -     .   \n",
       "\n",
       "                    gene_id       gene_type  ... transcript_type  \\\n",
       "1133585  ENSG00000120885.21  protein_coding  ...             NaN   \n",
       "\n",
       "        transcript_name transcript_support_level         tag  \\\n",
       "1133585             NaN                      NaN  ncRNA_host   \n",
       "\n",
       "        havana_transcript exon_number exon_id  ont protein_id ccdsid  \n",
       "1133585               NaN         NaN     NaN  NaN        NaN    NaN  \n",
       "\n",
       "[1 rows x 24 columns]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gtf_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "547d6076-d857-46ad-b5d5-1ad9714f1ef5",
   "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.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}