{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dbbe425c-4a8b-439c-97a1-cee37a070f41",
   "metadata": {},
   "source": [
    "# TPCAV detailed usage tutorial\n",
    "\n",
    "This tutorial goes through the detailed steps of TPCAV, make sure you installed tpcav package first:\n",
    "\n",
    "`pip install tpcav`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3a7dc41c-124f-448a-a9fb-5625a54d4729",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2026-02-25 15:38:56--  https://hgdownload.gi.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz\n",
      "Resolving hgdownload.gi.ucsc.edu (hgdownload.gi.ucsc.edu)... 128.114.119.163\n",
      "Connecting to hgdownload.gi.ucsc.edu (hgdownload.gi.ucsc.edu)|128.114.119.163|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 983659424 (938M) [application/x-gzip]\n",
      "Saving to: ‘data/hg38.fa.gz’\n",
      "\n",
      "hg38.fa.gz          100%[===================>] 938.09M  39.0MB/s    in 26s     \n",
      "\n",
      "2026-02-25 15:39:23 (36.4 MB/s) - ‘data/hg38.fa.gz’ saved [983659424/983659424]\n",
      "\n",
      "--2026-02-25 15:40:00--  https://raw.githubusercontent.com/seqcode/TPCAV/main/data/motif-clustering-v2.1beta_consensus_pwms.test.meme\n",
      "Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.111.133, 185.199.109.133, 185.199.108.133, ...\n",
      "Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.111.133|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 3504 (3.4K) [text/plain]\n",
      "Saving to: ‘data/motif-clustering-v2.1beta_consensus_pwms.test.meme’\n",
      "\n",
      "motif-clustering-v2 100%[===================>]   3.42K  --.-KB/s    in 0s      \n",
      "\n",
      "2026-02-25 15:40:00 (25.5 MB/s) - ‘data/motif-clustering-v2.1beta_consensus_pwms.test.meme’ saved [3504/3504]\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Fasta(\"data/hg38.fa\")"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# data download\n",
    "!mkdir data/\n",
    "!wget https://hgdownload.gi.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz -P data/\n",
    "!gunzip data/hg38.fa.gz\n",
    "!wget https://raw.githubusercontent.com/seqcode/TPCAV/main/data/motif-clustering-v2.1beta_consensus_pwms.test.meme -P data/\n",
    "\n",
    "import pyfaidx\n",
    "pyfaidx.Fasta('data/hg38.fa')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fbed531-7834-457b-a9fe-c878537a026c",
   "metadata": {},
   "source": [
    "## Prepare the model\n",
    "\n",
    "Since TPCAV is a concept based attribution method, the first step is to construct candidate concepts for testing, each concept is a set of input examples that share a similar pattern (e.g. motif).\n",
    "\n",
    "Assume we start with a dummy model that has two linear layers taking 1024bp long one-hot coded DNA input for predicting a scalar value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bf92bb3c-89de-472e-8872-3a54ce835476",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "\n",
    "class DummyModelSeq(torch.nn.Module):\n",
    "    def __init__(self):\n",
    "        super().__init__()\n",
    "        self.layer1 = torch.nn.Linear(1024, 1)\n",
    "        self.layer2 = torch.nn.Linear(4, 1)\n",
    "\n",
    "    def forward(self, seq):\n",
    "        y_hat = self.layer1(seq)\n",
    "        y_hat = y_hat.squeeze(-1)\n",
    "        y_hat = self.layer2(y_hat)\n",
    "        return y_hat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80d1891d-ab65-4299-aa6a-723b7150f233",
   "metadata": {},
   "source": [
    "## Construct concepts\n",
    "\n",
    "Now we start constructing concepts, a concept is basically an iterator of the inputs that represent a prominent pattern. TPCAV provides a class `ConceptBuilder` to construct common concepts in genomics field, here we use it to construct motif concepts, defined as sets of random genomic sequences inserted by motif instances.\n",
    "\n",
    "There are 10 motif concepts built in this case given the test meme file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "16ba9d74-a736-4ffb-8c9c-724ab68c7699",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jmy5455/miniforge3/envs/test/lib/python3.14/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[(Concept(1, 'AC0001:GATA-PROP:GATA'),\n",
       "  Concept(2, 'AC0001:GATA-PROP:GATA_perm')),\n",
       " (Concept(3, 'AC0002:PROP-ALX:Homeodomain'),\n",
       "  Concept(4, 'AC0002:PROP-ALX:Homeodomain_perm')),\n",
       " (Concept(5, 'AC0003:HNF1A-HNF1B:Homeodomain'),\n",
       "  Concept(6, 'AC0003:HNF1A-HNF1B:Homeodomain_perm')),\n",
       " (Concept(7, 'AC0004:ZSCAN:C2H2_ZF'), Concept(8, 'AC0004:ZSCAN:C2H2_ZF_perm')),\n",
       " (Concept(9, 'AC0005:POU3F-POU1F:Homeodomain,POU'),\n",
       "  Concept(10, 'AC0005:POU3F-POU1F:Homeodomain,POU_perm'))]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from pathlib import Path\n",
    "from tpcav.concepts import ConceptBuilder\n",
    "\n",
    "motif_path = Path(\"data\") / \"motif-clustering-v2.1beta_consensus_pwms.test.meme\"\n",
    "\n",
    "# create concept builder to generate concepts\n",
    "builder = ConceptBuilder(\n",
    "    genome_fasta=\"data/hg38.fa\",\n",
    "    input_window_length=1024,\n",
    "    bws=None,\n",
    "    num_motifs=12,\n",
    "    include_reverse_complement=True,\n",
    "    min_samples=1000,\n",
    "    batch_size=8,\n",
    ")\n",
    "# use random regions as control  \n",
    "builder.build_control()\n",
    "# use meme motif PWMs to build motif concepts, one concept per motif\n",
    "builder.add_meme_motif_concepts(str(motif_path))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71e96db7-d9a8-47b5-9d7e-d4eae2b631dc",
   "metadata": {},
   "source": [
    "Concept constructed by `ConceptBuilder` contains an iterator of two things: fasta sequences strings and array of bigwig signal tracks. If your model takes different formats of inputs, you can provide a transformation function to fit your model. Make sure your transformation function returns a tuple and your model forward function can be called by `model(*inputs)`\n",
    "\n",
    "Below is the example of obtaining one hot coded DNA sequence inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2ef8bb2a-c6bc-49b6-a369-3e672e6ef5af",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jmy5455/miniforge3/envs/test/lib/python3.14/site-packages/torch/utils/data/dataloader.py:668: UserWarning: 'pin_memory' argument is set as true but no accelerator is found, then device pinned memory won't be used.\n",
      "  warnings.warn(warn_msg)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(['ATTTTACTAACTGCCACACTGCATAGCTTTTATAAGCATATTAATAGATTTACGGGAGAGAGGCACCAGAGTGGATTACAACACAGGACTGACCTTTTACCACCCGGCTGAAGCTCGAAAAACAAAACCAAATCCCATCTTACTAACCATACCATAACTAGGTTATGGAAATACCCCAGAATACAAGTTTGTTTGTTTTTTAAGATGGTTTTCCTTCCCTGAGAGAACAATCCTGAAACTTAAGTCAATAAAGATACTGCCATTATTAAGCATTTTGCCTCCTCGGGTCAGACTTGAGGCCTTTTCCCCCCTTTCTTCTCGGTAATAAGCTTACTGCCTTAGATCACGGTGCATCCTAGGCAGGTTATGCCATTTTAATCTCATAAATATGAAAATAAAAGACATAATCAGAGCAATGTCGTTATTTTTTATTTTTAAACCAAAGGAAGCCAAAATAGCAACAGAGAAGAAAGGAAATATAATATCATTAACTCTGGGGCCAAGGTGCATTTGCCCTAGTTCTTTTTTTAATAATGCGCTTATTTTTCTTCCTCCCCATTCCACAGCACTCAGCATCTTCTATCTTGTCTCCCTTGGTGTAGCTCTAGTGTTCTAATTCTATCCCATCTTCTTCTACTGAAATAATGAACTTTTTAAAGAGAAAATTAAGTGAATGTCTTTGTAGCAATGAACCTCCCCCTCTCGAGAATAAAAGGGAGCTTTGGAGATAAGGCAATGCTTAACCCAGTGAACGAATTGCTTACTGATGGGGTCAGGTAAAAGGATTACAGGGTGGCAGAGCTGACTCTGCCTACTGATATCGGCCATGTTAAATAATTCCCGTATTATATACTCAACATCTATATAACAATGTTTTTGTTAGCTTGTTTCTATTTTATAACAGCCCTGAAATGTAAAATGGTTTTATTTTATTTATTTATTATAATGCTGTAAGTTTTTTGAGATAATGACATAACCTGTCGCCCAGGCTGGAGTGCAGTGGCGCGACCTCGGCTCACTGCAA',\n",
       "  'ATTTTACTAACTGCCACACTGCATAGCTTTTATAAGAATATTAATAGATTTACGGGAGAGAGGCACCAGAGTGGATTACAACACAGGACTGACCTTTTACCACCCGGCTGAAGCTCGAAAAACAAAACCAAATCCCATCTTACTAACCATACCATAACTACGTTATGGAAATACCCCAGAATACAAGTTTGTTTGTTTTTTAAGATGGTTTTCCTTCCCTGAGAGAACAATCCTGAAACTTAAGTCAATAAAGATACTCGCATTATTAAGCATTTTGCCTCCTCGGGTCAGACTTGAGGCCTTTTCCCCCCTTTCTTCTCGGTAATAACCTTACTGCCTTAGATCACGGTGCATCCTAGGCAGGTTATGCCATTTTAATCTCATAAATATGAAAATAAAAGACATAATCAGAGACATGTCGTTATTTTTTATTTTTAAACCAAAGGAAGCCAAAATAGCAACAGAGAAGAAAGGAAATATAATATCCTTAACTCTGGGGCCAAGGTGCATTTGCCCTAGTTCTTTTTTTAATAATGGCGTTATTTTTCTTCCTCCCCATTCCACAGCACTCAGCATCTTCTATCTTGTCTCCCTTGGTGTAGCTCTAGTGTTCTAATTCTATCCCATCTTCTTCTACTGAAATAATGAACTTTTTAAAGAGAAAATTAAGTGAATGTCTTTGTAGCAATGAACCTCCCCCTCTCGAGAATAAAAGGGAGCTTTGGAGATAAGGCAATGCTTAACCCAGTGAACGAATTGCTTACTGATGGGGTCAGGTAAAATCATTACAGGGTGGCAGAGCTGACTCTGCCTACTGATATCGGCCATGTTAAATAATTGCCGTATTATATACTCAACATCTATATAACAATGTTGTTGTTAGCTTGTTTCTATTTTATAACAGCCCTGAAATGTAAAATGGTTTTATTTTATTTATTTATTATAATTGTGTTAGTTTTTTGAGATAATGACATTACCTGTCGCCCAGGCTGGAGTGCAGTGGCGCGACCTCGGCTCACTGCAA',\n",
       "  'CTCTTATTTCCTGCTATTATGCCTAATAGCGTTATACCTGTAGTAACCTATCTGTAGTTACCTGTAAGTAACTCTAATAAACTGATCACTAAAATAGAGTGTCTCAAAAGGATATTTTAAGCTGATAAGATCCAAATATAAAGAAAGACGGAAATACATGGAATTGTCTGGGGCCAAGAAATACTTGATAATTTTTGCCAGTAAAGTAAGAGTTGCCATAGTTCCCTTGGCTTTAGGGCATTGGCTATAAGGTCATAATCTTGCCCCTAGAGATAAATAAGGTTATTATAACTATGAAACATTATTTTATTAAAAGTGAATATACAAAATAACTCAGAAATACATGATCATATTTTGTCTGTGGAGTACAAGCATTTATTTGTCTGCTGATTACTTCTTTGGGGGAACTACCCATCCTTCATTCTTGGTCAGTATAAGGGAGTTATCTAAAGTAACCATGGCTCCAGGGGTGGGCATGTGATCCAAATTTTGCTAATTATTATTAGACATCCTTGGAGGGATATAACGACGTTATGGTGAGCATATTACTCAAGCTTGCCAATTAATGTCGTTATTTGAGAATTTAGCAAAAATTACTGGGAAAGAGGCATTTCCTTTCCATTTGAGTTACTGAACTGATAGAAAGTAACTGATAAACTCTATAATGCCATTATGGCCGTTCCTACCCTCTGCATCGGGAAGGCATATTTGACAATAAAGCTGATACAGAAGAAAACAGATTTAAGACGTGGAGAAACAGTGTTCTGATAAGGCAATTATAGCAAATTGATTCAACACCCAATTTAGTCATGTGAGGAAATGAGTTCTTTTTTGTTTTATGATAATGCTATTATGGTTTTTGTGTAATGGTATAATGGAAATGAGTTCTTTTTTGTTTTATGCATAACGATATTTTGTTTTTGTCCCCTGTAACTGAAAGGAGCTTGGCTAATTAAGTCAGGAATTTGCATATAGCTTTTAGATATAAATTAAATAGTTTTGTAATTGATATGCTTACTTTTCT',\n",
       "  'TAAGTAGTATTCTATCGAATATATATATACCACAGTTTCTTTATGCACTCATTGATTGTTGGGCATTTGGGTTGGCTCTGTTTTTGCAATTGTGAATTGTGCTGCTATATACATACTATGATTAATATCTTTATAAGCTGTTTAGTTCTTTTACTGCGATTACATACCCAGTAGTGGGATTGCTGGATCAAATGGTAGTTTTACTTTTAGTCCTTTAAGCAATCTCTCCACGCTGTTTTCCACAGTGGTTGTTTGTTTACATTCCCATTAGCAGTGCAGAAGTGTTCCCTATTATAAAGCTCTTTTCCAACATCTACTATTTTTACAATAATCGTACTATGGCCATTCTAACAGAAGTAAGGTGATATCACACTGTGCTTTTGATTTGTGTTTCCCTGATCATTAGTGAATAACTCCATTATTTCATATGTTTGTTGGCATAAGATCGTTAATCTTTTGACTAATGTCCTTATATGTCCTTAGCCCACTTTTTGATGGGATTATAACACCATTATTTGTTGATTTGTTTGAGTTCATTGTAGAATCTGGATATTAGTCCTCTGTCAGGTGTATAGATTGTAAAGATTTTTTCCCATTCTGTGGGTTGTCTGTTTACTTTGCTGACTGTCCCTTTTCCCACTCAAAAGTTTCTAGCAACCGCAGTTAAAAAAGACAAAGTGGAACATTATGTAATGATAAAAGGACTAGTCCAACAAGAAAATATCACAATCCTAGCCATACATGCACCTAACACTGGAGCTCCCAGATTTGTAAAACGATTACTAATATAATGGCATTATGAGATAGACAGCAACACAGTAATAGCTTAATACAGTTATACTCCACTGACAGCACTAGGCAGGTCATTAAGACAGAAACTCAAGAAAGAACAATGGATTTAAACTACACTCCGGAACAAATGGACTTAACACCTATATACAGAGCATTTTATCCAATAACCGCAGAATATACATTCTTTTCAACAGCATATGTTAAGGCTGTTATGATAGGTCATATGACAGGC',\n",
       "  'CTCTTATTTCCTGCTATTATGCCCAATAGCGTTATACCTGTAGTAACCTATCTGTAGTTACCTGTAAGTAACTCTAATAAACTTATCACTAAAATAGAGTGTCTCAAAAGGATATTTTAAGCTGATAAGATCCAAATATAAAGAAAGACGGAAATACATGGAATTGTCTGGGGCCAAGAAATACTTGATAATTTTTGCCAGTAAAGTAAGAGTTGCCATAGTTCCCTTGGCTTTAGGGCATTGGCTATAAGGTCATTATCTTGCCCCTAGAGATAAATAAGCTTATTATAACTATGAAACATTATTTTATTAAAAGTGAATATACAAAATAACTCAGAAATACATGATCATATTTTGTCTGTGGAGTACAAGCATTTATTTGTCTGCTGATTACTTCTTTGGGGGAACTACCCATCCTTCATTCTTGGTCAGTATCATGGCGTTATCTAAAGTAACCATGGCTCCAGGGGTGGGCATGTGATCCAAATTTTGCTAATTATTATTAGACATCCTTGGAGGGATATAACGACGTTATGGTGAGCATATTACTCAAGCTTGCCAATTAATGTCGTTATTTGAGAATTTAGCAAAAATTACTGGGAAAGAGGCATTTCCTTTCCATTTGAGTTACTGAACTGATAGAAAGTAACTGATAAACTCTATAATGGCATTATGGCCGTTCCTACCCTCTGCATCGGGAAGGCATATTTGACAATAAAGCTGATACAGAAGAAAACAGATTTAAGACGTGGAGAAACAGTGTTCTGATAAGGCACTTATAGCAAATTGATTCAACACCCAATTTAGTCATGTGAGGAAATGAGTTCTTTTTTGTTTTATGATTATGCTATTATGGTTTTTGTGTAATGGTATTATGGAAATGAGTTCTTTTTTGTTTTATGCATAACGAGATTGTGTTTTTGTCCCCTGTAACTGAAAGGAGCTTGGCTAATTAAGTCAGGAATTTGCATATAGCTTTTAGATATAAATTAAATAGTTTTGTAATTGATATGCTTACTTTTCT',\n",
       "  'TAAGTAGTATTCTATCGAATATATATATACCACAGTTTCTTTATGCACTCATTGATTGTTGGGCATTTGGGTTGGCTCTGTTTTTGCAATTGTGAATTGTGCTGCTATATACATACTATCATTAATATCTTTCTAAGATCTTTATTTCTTTTACTTGCCTTACATACCCAGTAGTGGGATTGCTGGATCAAATGGTAGTTTTACTTTTAGTCCTTTAAGCAATCTCTCCACGCTGTTTTCCACAGTGGTTGTTTGTTTACATTCCCATTAGCAGTGCAGAAGTGTTCCCTATTATAAACCGCTTGTCCAACATCTACTATTTTTAAAATAATGGTACTATGGCCATTCTAACAGAAGTAAGGTGATATCACACTGTGCTTTTGATTTGTGTTTCCCTGATCATTAGTGAATAACTCCATTATTTCATATGTTTGTTGGCATAAGATCGTTAATCTTTTGACTAATGTCCTTATATGTCCTTAGCCCACTTTTTGATGGGATTATAACACCATTATTTGTTGATTTGTTTGAGTTCATTGTAGAATCTGGATATTAGTCCTCTGTCAGGTGTATAGATTGTAAAGATTTTTTCCCATTCTGTGGGTTGTCTGTTTACTTTGCTGACTGTCCCTTTTCCCACTCAAAAGTTTCTAGCAACCGCAGTTAAAAAAGACAAAGTGGAACATTATGTAATGATAAAAGGACTAGTCCAACAAGAAAATATCACAATCCTAGCCATACATGCACCTAACACTGGAGCTCCCAGATTTGTAAAACGATTACTAATATAATTGCATTATGAGATAGACAGCAACACAGTAATAGCTTAATACCGTTATACTCCACTGACAGCACTAGGCAGGTCATTAAGACAGAAACTCAAGAAAGAACAATGGATTTAAACTACACTCCGGAACAAATGGACTTAACACCTATATACAGAGCATTTTATCCAATAACCGCAGAATATACATTCTTTTCAACAGCATATGTTAATGGTGTTATGATAGGTCATATGACAGGC',\n",
       "  'TTATCAAGTATCTGCAGCTAATGGACTTATCTCAAACCTAAAAATAAGAAATGAAAAAAATGCCTAAGAAAATTGTAATATAACGAAGTTTTTCATCACACATGCTAAAACTTTTCTGGATGAAAGATACAAACTGGCCATTAATGCGGTTATTGAATGGAACTTAACCGAGGAGTTTGGCAGAGGAATGGCTGGGTTATGTAATGCCCGTATTTTCTCTGCATTTAGCTGAAAGCATGGTTGTCCAGAATTAAGGCTATTTCTCAGTCTATAATTATATTATGCAATGGTCATGTACCTCTGCATTTAGCTGAAAGCATGGTTGTCCAGAATTAAATAACGAAATTATTCTCTCTTAACAGCTGGCAATGGTCATGTAACTAAGTTCTGGGTAATGATATGTTGAGTGGGAGTGGTGTGTGCAGTATTCTGATTGAGCCTGTGGAGGGAAGCATGTGCTCTTCACTTCCTGTTTCCTCCCTTCTAGCTGCCTGGAACATAAACTTAGTAATGAAATTACTTAGGAAAAATCGTAACAACGTTATCTTCCTTGCTCTTAATCTTATTAAATGAATAATAGAACTTAAACATAGGTATTTTGACACAAATTATCCACTAAGCATATAATTCAAAATATTCGATTGAATTCTTTACTCCCACAATAATTTAGCATTCTAAAACATTAATTTGCTGAAGAGAAATTGGCTTTTGGATGAGTTAGGGCATTTACTAAAAGGTTCCTTTATTACTGAAGAGTGGTTAGTAATAATGAGATTAAAGATAATTTGAAATTATTAATTTTATGAGATGATCTTTAAAAATACTTTTTGAAAATCAAATTATGCTTTCAGATTTTCCAGGTCAAAATTTGTATTGTAATTTTATTACAATTGTGTTAGAATTAAACACACATACAATCAAATTATGAAAGACAGACACAGTTGGAAGAGTCATACTCCCTCATTATAAAAGTTACTACAAAATTATAGAATTCAAATAATCTACTTATTGGCAAAAGGATAAA',\n",
       "  'TTATCAAGTATCTGCAGCCTATGTACTTATCTCAAACCTAAAAATAAGAAATGAAAAAAATGCCTAAGAAAATTGTAATATAACGACGTTGTTCATCACACATGCTAAAACTTTTCTGGATGAAAGATACAAACTGGCCATTAATGCCGTTATTGAATGGAACTTAACCGAGGAGTTTGGCAGAGGAATGGCTGGGTTATGTAATTCCCGTATTTTCTCTGCATTTAGCTGAAAGCATGGTTGTCCAGAATTAAGGCTATTTCTCAGTCTATAATTATCTTATGCAATGGTCATGTACCTCTGCATTTAGCTGAAAGCATGGTTGTCCAGAATTAAATAACCAAATTATTCTCTCTTAACAGCTGGCAATGGTCATGTAACTAAGTTCTGGGTAATGATATGTTGAGTGGGAGTGGTGTGTGCAGTATTCTGATTGAGCCTGTGGAGGGAAGCATGTGCTCTTCACTTCCTGTTTCCTCCCTTCTAGCTGCCTGGAACATAAACTTAGTAATGAAATTACTTAGGAAAAATCGTAACACCGTTATCTTCCTTGCTCTTAATCTTATTAAATGAATAATAGAACTTAAACATAGGTATTTTGACACAAATTATCCACTAAGCATATAATTCAAAATATTCGATTGAATTCTTTACTCCCACAATAATTTAGCATTCTAAAACATTAATTTGCTGAAGAGAAATTGGCTTTTGGATGAGTTAGGGCATTTACTAAAAGGTTCCTTTATTACTGAAGAGTGGTTAGTAATAATGACATTAAAGATAATTTGAAATTATTAATTTTATGAGATGATCTTTAAAAATACTTTTTGAAAATCAAATTATGCTTTCAGATTTTCCAGGTCAAAATTTGTATTGTAATTTTATTACAATTGTGTTAGAATTAAACACACATAAAATCAAATTATGAAAGACAGACACAGTTGGAAGAGTCATACTCCCTCATTATAAAAGTTACTACAAAATTATAGAATTCAAATAATCTACTTATTGGCAAAAGGATAAA'],\n",
       " tensor([[nan],\n",
       "         [nan],\n",
       "         [nan],\n",
       "         [nan],\n",
       "         [nan],\n",
       "         [nan],\n",
       "         [nan],\n",
       "         [nan]]))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# each batch is a tuple of a fasta sequence list and an array of bigwig signals\n",
    "next(iter(builder.concepts[0].data_iter))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0504fa0f-f05b-4878-80d7-58c697cb5344",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tpcav import helper\n",
    "\n",
    "# apply transformation function to obtain one-hot encoded sequences\n",
    "def transform_fasta_to_one_hot_seq(seq, chrom):\n",
    "    # `seq` is a list of fasta sequences\n",
    "    # `chrom` is a numpy array of bigwig signals of shape [-1, # bigwigs, len]\n",
    "    return (helper.fasta_to_one_hot_sequences(seq),) # it has to return a tuple of inputs, even if there is only one input\n",
    "builder.apply_transform(transform_fasta_to_one_hot_seq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "689f99f3-ebda-4fb2-8661-44577017c198",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(tensor([[[1., 0., 0.,  ..., 0., 1., 1.],\n",
       "          [0., 0., 0.,  ..., 1., 0., 0.],\n",
       "          [0., 0., 0.,  ..., 0., 0., 0.],\n",
       "          [0., 1., 1.,  ..., 0., 0., 0.]],\n",
       " \n",
       "         [[0., 0., 0.,  ..., 0., 0., 0.],\n",
       "          [1., 0., 1.,  ..., 0., 1., 0.],\n",
       "          [0., 0., 0.,  ..., 0., 0., 0.],\n",
       "          [0., 1., 0.,  ..., 1., 0., 1.]],\n",
       " \n",
       "         [[1., 0., 0.,  ..., 0., 1., 1.],\n",
       "          [0., 0., 0.,  ..., 1., 0., 0.],\n",
       "          [0., 0., 0.,  ..., 0., 0., 0.],\n",
       "          [0., 1., 1.,  ..., 0., 0., 0.]],\n",
       " \n",
       "         ...,\n",
       " \n",
       "         [[0., 0., 1.,  ..., 1., 1., 1.],\n",
       "          [0., 0., 0.,  ..., 0., 0., 0.],\n",
       "          [0., 0., 0.,  ..., 0., 0., 0.],\n",
       "          [1., 1., 0.,  ..., 0., 0., 0.]],\n",
       " \n",
       "         [[0., 0., 0.,  ..., 0., 1., 0.],\n",
       "          [0., 0., 0.,  ..., 0., 0., 0.],\n",
       "          [1., 1., 1.,  ..., 0., 0., 0.],\n",
       "          [0., 0., 0.,  ..., 1., 0., 1.]],\n",
       " \n",
       "         [[1., 0., 0.,  ..., 0., 0., 0.],\n",
       "          [0., 0., 0.,  ..., 1., 1., 0.],\n",
       "          [0., 0., 0.,  ..., 0., 0., 0.],\n",
       "          [0., 1., 1.,  ..., 0., 0., 1.]]]),)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# after transformation each batch is one-hot coded DNA arrays wrapped in tuple\n",
    "next(iter(builder.concepts[0].data_iter))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f476a6dd-c2b8-42a4-8580-6cafba3dd43f",
   "metadata": {},
   "source": [
    "## Train concept activaton vector (CAV) classifier\n",
    "\n",
    "After concept construction, now it's time to train the linear classifier for every concept in your model. You need to wrap your model instance by `TPCAV` class and the layer for interpretation provided, then we apply PCA transformation given the concepts we just built to decorrelate the embedding space to solve the correlated redundant feature issue.\n",
    "\n",
    "Note: You can either provide a layer name (should be able to retrieve the corresponding layer by `getattr(model, layer_name)`) or the exact layer variable (make sure it's from the instantiated model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5d492620-87d0-45e4-9764-eb7943cc099e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:tpcav.tpcav_model:Sampled 10 activations from concept AC0001:GATA-PROP:GATA\n",
      "INFO:tpcav.tpcav_model:Sampled 10 activations from concept AC0002:PROP-ALX:Homeodomain\n",
      "INFO:tpcav.tpcav_model:Sampled 10 activations from concept AC0003:HNF1A-HNF1B:Homeodomain\n",
      "INFO:tpcav.tpcav_model:Sampled 10 activations from concept AC0004:ZSCAN:C2H2_ZF\n",
      "INFO:tpcav.tpcav_model:Sampled 10 activations from concept AC0005:POU3F-POU1F:Homeodomain,POU\n",
      "INFO:tpcav.tpcav_model:Sampled 10 activations from concept random_regions\n",
      "INFO:tpcav.tpcav_model:Sampled 10 activations from concept AC0001:GATA-PROP:GATA_perm\n",
      "INFO:tpcav.tpcav_model:Sampled 10 activations from concept AC0002:PROP-ALX:Homeodomain_perm\n",
      "INFO:tpcav.tpcav_model:Sampled 10 activations from concept AC0003:HNF1A-HNF1B:Homeodomain_perm\n",
      "INFO:tpcav.tpcav_model:Sampled 10 activations from concept AC0004:ZSCAN:C2H2_ZF_perm\n",
      "INFO:tpcav.tpcav_model:Sampled 10 activations from concept AC0005:POU3F-POU1F:Homeodomain,POU_perm\n"
     ]
    }
   ],
   "source": [
    "from tpcav.tpcav_model import TPCAV\n",
    "# create TPCAV model on top of your model\n",
    "tpcav_model = TPCAV(DummyModelSeq(), layer_name=\"layer1\")\n",
    "# fit PCA on sampled all concept activations\n",
    "tpcav_model.fit_pca(\n",
    "    concepts=builder.all_concepts(),\n",
    "    num_samples_per_concept=10,\n",
    "    num_pc=\"full\",\n",
    ")\n",
    "# you can save the tpcav model in case of future use\n",
    "torch.save(tpcav_model, \"data/tmp_tpcav_model.pt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7177e0eb-6e56-4141-b5ab-94fed823f282",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[-0.1796,  0.0593, -0.6600,  0.7271],\n",
       "        [ 0.7280, -0.6324, -0.2646, -0.0088],\n",
       "        [-0.4566, -0.7081,  0.4251,  0.3308],\n",
       "        [-0.4789, -0.3085, -0.5601, -0.6015]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# now tpcav_model contains the necessary parameters for PCA transformation\n",
    "tpcav_model.Vh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3d8dbd7-1acb-4811-9cbd-581781f68126",
   "metadata": {},
   "source": [
    "Using the fitted TPCAV model, we can train the concept activation vectors (CAVs) on the concepts we built, here we use the random genomic background concept as control"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6c07c332-ee51-4925-80f2-ab06e397715f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:tpcav.cavs:Submitted CAV training for concept AC0001:GATA-PROP:GATA\n",
      "INFO:tpcav.cavs:Submitted CAV training for concept AC0002:PROP-ALX:Homeodomain\n",
      "INFO:tpcav.cavs:Submitted CAV training for concept AC0003:HNF1A-HNF1B:Homeodomain\n",
      "INFO:tpcav.cavs:Submitted CAV training for concept AC0004:ZSCAN:C2H2_ZF\n",
      "INFO:tpcav.cavs:Submitted CAV training for concept AC0005:POU3F-POU1F:Homeodomain,POU\n",
      "INFO:tpcav.cavs:Best Params: {'alpha': 1e-06} | Iterations: 20\n",
      "INFO:tpcav.cavs:Best Params: {'alpha': 0.01} | Iterations: 18\n",
      "INFO:tpcav.cavs:[train] Accuracy: 0.5536\n",
      "INFO:tpcav.cavs:[val] Accuracy: 0.5500\n",
      "INFO:tpcav.cavs:[train] Accuracy: 0.5417\n",
      "INFO:tpcav.cavs:[val] Accuracy: 0.6000\n",
      "INFO:tpcav.cavs:[test] Accuracy: 0.5000\n",
      "INFO:tpcav.cavs:[test] Accuracy: 0.4500\n",
      "INFO:tpcav.cavs:Best Params: {'alpha': 1e-08} | Iterations: 14\n",
      "INFO:tpcav.cavs:Best Params: {'alpha': 0.0001} | Iterations: 18\n",
      "INFO:tpcav.cavs:[train] Accuracy: 0.5893\n",
      "INFO:tpcav.cavs:[train] Accuracy: 0.5833\n",
      "INFO:tpcav.cavs:[val] Accuracy: 0.4000\n",
      "INFO:tpcav.cavs:[val] Accuracy: 0.3000\n",
      "INFO:tpcav.cavs:[test] Accuracy: 0.3500\n",
      "INFO:tpcav.cavs:[test] Accuracy: 0.4500\n",
      "INFO:tpcav.cavs:Best Params: {'alpha': 1e-06} | Iterations: 13\n",
      "INFO:tpcav.cavs:[train] Accuracy: 0.5893\n",
      "INFO:tpcav.cavs:[val] Accuracy: 0.5500\n",
      "INFO:tpcav.cavs:[test] Accuracy: 0.6000\n"
     ]
    }
   ],
   "source": [
    "from tpcav.cavs import CavTrainer\n",
    "\n",
    "# create trainer for computing CAVs\n",
    "cav_trainer = CavTrainer(tpcav_model, penalty=\"l2\")\n",
    "# set control concept for CAV training\n",
    "cav_trainer.set_control(builder.control_concepts[0], num_samples=100)\n",
    "# train CAVs for all concepts\n",
    "cav_trainer.train_concepts(\n",
    "    builder.concepts, 100, output_dir=\"data/cavs/\", num_processes=2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "684ebb1f-a0b6-4c87-9142-a4fce408206d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# after training, the cavs are stored\n",
    "cav_trainer.cav_weights\n",
    "# you can save the cav_trainer for future use\n",
    "torch.save(cav_trainer, \"cav_trainer.pt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9f92e1d-7af4-4b75-bc3f-13a498764a01",
   "metadata": {},
   "source": [
    "## Layer attribution & TPCAV score calculation\n",
    "\n",
    "Finally, we can compute the layer attribution scores on the test regions against control regions, and obtain the final TPCAV score for each concept\n",
    "\n",
    "Note the default layer attribution in TPCAV uses DeepLIFT multiplier x abs(input - baselines), which is different from default Captum implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "81403bb9-5649-4307-94d8-d44a6f523241",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jmy5455/miniforge3/envs/test/lib/python3.14/site-packages/captum/attr/_core/deep_lift.py:294: UserWarning: Input Tensor 1 did not already require gradients, required_grads has been set automatically.\n",
      "  gradient_mask = apply_gradient_requirements(inputs_tuple)\n",
      "/home/jmy5455/miniforge3/envs/test/lib/python3.14/site-packages/captum/log/dummy_log.py:39: UserWarning: Setting forward, backward hooks and attributes on non-linear\n",
      "               activations. The hooks and attributes will be removed\n",
      "            after the attribution is finished\n",
      "  return func(*args, **kwargs)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'AC0001:GATA-PROP:GATA': 3.912023005428146,\n",
       " 'AC0002:PROP-ALX:Homeodomain': 4.61512051684126,\n",
       " 'AC0003:HNF1A-HNF1B:Homeodomain': 1.0216512475319812,\n",
       " 'AC0004:ZSCAN:C2H2_ZF': 2.9652730660692823,\n",
       " 'AC0005:POU3F-POU1F:Homeodomain,POU': 3.4965075614664802}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create input regions and baseline regions for attribution\n",
    "random_regions_1 = helper.random_regions_dataframe(\n",
    "    \"data/hg38.fa.fai\", 1024, 100, seed=1\n",
    ")\n",
    "random_regions_2 = helper.random_regions_dataframe(\n",
    "    \"data/hg38.fa.fai\", 1024, 100, seed=2\n",
    ")\n",
    "# create iterators to yield one-hot encoded sequences from the region dataframes\n",
    "def pack_data_iters(df):\n",
    "    seq_fasta_iter = helper.dataframe_to_fasta_iter(\n",
    "        df, \"data/hg38.fa\", batch_size=8\n",
    "    )\n",
    "    seq_one_hot_iter = (\n",
    "        helper.fasta_to_one_hot_sequences(seq_fasta)\n",
    "        for seq_fasta in seq_fasta_iter\n",
    "    )\n",
    "    return zip(seq_one_hot_iter, )\n",
    "# compute layer attributions given the iterators of testing regions and control regions\n",
    "attributions = tpcav_model.layer_attributions(\n",
    "    pack_data_iters(random_regions_1), pack_data_iters(random_regions_2)\n",
    ")\n",
    "# compute TPCAV scores for the concept\n",
    "cav_trainer.tpcav_score_all_concepts_log_ratio(attributions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11f64f51-dd1d-48c3-b75f-70d5fa8fee95",
   "metadata": {},
   "source": [
    "You can visualize the cosine similarity matrix between concepts vectors combined with the TPCAV scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "37a3ff15-bfc3-436a-9866-37530e4aac25",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 50428 characters]",
      "text/plain": [
       "<Figure size 1000x1000 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cav_trainer.plot_cavs_similaritiy_heatmap(attributions, fscore_thresh=0.0)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "test",
   "language": "python",
   "name": "test"
  },
  "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.14.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
