{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Phylogenetic Invariants" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AAAA AAAC AAAG AAAT AACA AACC AACG AACT AAGA AAGC AAGG AAGT AATA AATC AATG AATT\n", "ACAA ACAC ACAG ACAT ACCA ACCC ACCG ACCT ACGA ACGC ACGG ACGT ACTA ACTC ACTG ACTT\n", "AGAA AGAC AGAG AGAT AGCA AGCC AGCG AGCT AGGA AGGC AGGG AGGT AGTA AGTC AGTG AGTT\n", "ATAA ATAC ATAG ATAT ATCA ATCC ATCG ATCT ATGA ATGC ATGG ATGT ATTA ATTC ATTG ATTT\n", "\n", "CAAA CAAC CAAG CAAT CACA CACC CACG CACT CAGA CAGC CAGG CAGT CATA CATC CATG CATT\n", "CCAA CCAC CCAG CCAT CCCA CCCC CCCG CCCT CCGA CCGC CCGG CCGT CCTA CCTC CCTG CCTT\n", "CGAA CGAC CGAG CGAT CGCA CGCC CGCG CGCT CGGA CGGC CGGG CGGT CGTA CGTC CGTG CGTT\n", "CTAA CTAC CTAG CTAT CTCA CTCC CTCG CTCT CTGA CTGC CTGG CTGT CTTA CTTC CTTG CTTT\n", "\n", "GAAA GAAC GAAG GAAT GACA GACC GACG GACT GAGA GAGC GAGG GAGT GATA GATC GATG GATT\n", "GCAA GCAC GCAG GCAT GCCA GCCC GCCG GCCT GCGA GCGC GCGG GCGT GCTA GCTC GCTG GCTT\n", "GGAA GGAC GGAG GGAT GGCA GGCC GGCG GGCT GGGA GGGC GGGG GGGT GGTA GGTC GGTG GGTT\n", "GTAA GTAC GTAG GTAT GTCA GTCC GTCG GTCT GTGA GTGC GTGG GTGT GTTA GTTC GTTG GTTT\n", "\n", "TAAA TAAC TAAG TAAT TACA TACC TACG TACT TAGA TAGC TAGG TAGT TATA TATC TATG TATT\n", "TCAA TCAC TCAG TCAT TCCA TCCC TCCG TCCT TCGA TCGC TCGG TCGT TCTA TCTC TCTG TCTT\n", "TGAA TGAC TGAG TGAT TGCA TGCC TGCG TGCT TGGA TGGC TGGG TGGT TGTA TGTC TGTG TGTT\n", "TTAA TTAC TTAG TTAT TTCA TTCC TTCG TTCT TTGA TTGC TTGG TTGT TTTA TTTC TTTG TTTT\n", "\n", "[['AAAA' 'AAAC' 'AAAG' 'AAAT']\n", " ['ACAA' 'ACAC' 'ACAG' 'ACAT']\n", " ['AGAA' 'AGAC' 'AGAG' 'AGAT']\n", " ['ATAA' 'ATAC' 'ATAG' 'ATAT']]\n" ] } ], "source": [ "pi_string = \"\"\"\\\n", "AAAA AAAC AAAG AAAT AACA AACC AACG AACT AAGA AAGC AAGG AAGT AATA AATC AATG AATT\n", "ACAA ACAC ACAG ACAT ACCA ACCC ACCG ACCT ACGA ACGC ACGG ACGT ACTA ACTC ACTG ACTT\n", "AGAA AGAC AGAG AGAT AGCA AGCC AGCG AGCT AGGA AGGC AGGG AGGT AGTA AGTC AGTG AGTT\n", "ATAA ATAC ATAG ATAT ATCA ATCC ATCG ATCT ATGA ATGC ATGG ATGT ATTA ATTC ATTG ATTT\n", "\n", "CAAA CAAC CAAG CAAT CACA CACC CACG CACT CAGA CAGC CAGG CAGT CATA CATC CATG CATT\n", "CCAA CCAC CCAG CCAT CCCA CCCC CCCG CCCT CCGA CCGC CCGG CCGT CCTA CCTC CCTG CCTT\n", "CGAA CGAC CGAG CGAT CGCA CGCC CGCG CGCT CGGA CGGC CGGG CGGT CGTA CGTC CGTG CGTT\n", "CTAA CTAC CTAG CTAT CTCA CTCC CTCG CTCT CTGA CTGC CTGG CTGT CTTA CTTC CTTG CTTT\n", "\n", "GAAA GAAC GAAG GAAT GACA GACC GACG GACT GAGA GAGC GAGG GAGT GATA GATC GATG GATT\n", "GCAA GCAC GCAG GCAT GCCA GCCC GCCG GCCT GCGA GCGC GCGG GCGT GCTA GCTC GCTG GCTT\n", "GGAA GGAC GGAG GGAT GGCA GGCC GGCG GGCT GGGA GGGC GGGG GGGT GGTA GGTC GGTG GGTT\n", "GTAA GTAC GTAG GTAT GTCA GTCC GTCG GTCT GTGA GTGC GTGG GTGT GTTA GTTC GTTG GTTT\n", "\n", "TAAA TAAC TAAG TAAT TACA TACC TACG TACT TAGA TAGC TAGG TAGT TATA TATC TATG TATT\n", "TCAA TCAC TCAG TCAT TCCA TCCC TCCG TCCT TCGA TCGC TCGG TCGT TCTA TCTC TCTG TCTT\n", "TGAA TGAC TGAG TGAT TGCA TGCC TGCG TGCT TGGA TGGC TGGG TGGT TGTA TGTC TGTG TGTT\n", "TTAA TTAC TTAG TTAT TTCA TTCC TTCG TTCT TTGA TTGC TTGG TTGT TTTA TTTC TTTG TTTT\n", "\"\"\"\n", "\n", "## print as string\n", "print pi_string\n", "\n", "## also save as indexable array\n", "pi_arr = np.array(pi_string.split()).reshape(16,16)\n", "print pi_arr[:4,:4]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create a simple alignment for four taxa\n", "Ordered so the first two rows have a more similar sequence than the latter two (e.g., see CATCAT vs. GATGAT). " ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[['C' 'A' 'T' 'C' 'A' 'T' 'C' 'A' 'T' 'C']\n", " ['C' 'A' 'T' 'C' 'A' 'T' 'G' 'A' 'T' 'C']\n", " ['G' 'A' 'T' 'G' 'A' 'T' 'G' 'A' 'T' 'C']\n", " ['G' 'A' 'T' 'G' 'A' 'T' 'G' 'A' 'T' 'C']] ...\n" ] } ], "source": [ "seqalign = np.array([\n", " list(\"CATCATCATCAT\")*10,\n", " list(\"CATCATGATCAT\")*10,\n", " list(\"GATGATGATCCC\")*10,\n", " list(\"GATGATGATCCC\")*10])\n", "\n", "print seqalign[:, :10], '...'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A trick for dealing with large seq arrays\n", "We convert it to integers because they are easier to deal with than characters. The 'view' argument is used here to convert chars to 8bit integers. I then reassign the ints as 0, 1, 2, or 3, which can be used to quickly fill the the count array using a set of rules designated below. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[['C' 'A' 'T' 'C' 'A' 'T' 'C' 'A' 'T' 'C']\n", " ['C' 'A' 'T' 'C' 'A' 'T' 'G' 'A' 'T' 'C']\n", " ['G' 'A' 'T' 'G' 'A' 'T' 'G' 'A' 'T' 'C']\n", " ['G' 'A' 'T' 'G' 'A' 'T' 'G' 'A' 'T' 'C']]\n", "[[67 65 84 67 65 84 67 65 84 67]\n", " [67 65 84 67 65 84 71 65 84 67]\n", " [71 65 84 71 65 84 71 65 84 67]\n", " [71 65 84 71 65 84 71 65 84 67]]\n", "[[1 0 3 1 0 3 1 0 3 1]\n", " [1 0 3 1 0 3 2 0 3 1]\n", " [2 0 3 2 0 3 2 0 3 1]\n", " [2 0 3 2 0 3 2 0 3 1]]\n" ] } ], "source": [ "## print as char\n", "print seqalign[:, :10]\n", "\n", "## get a copy with dtype=int8\n", "seqint = np.copy(seqalign.view(np.int8))\n", "print seqint[:, :10]\n", "\n", "## convert ints to indexes of matrix\n", "seqint[seqint == 65] = 0 ## As\n", "seqint[seqint == 67] = 1 ## Cs\n", "seqint[seqint == 71] = 2 ## Ts\n", "seqint[seqint == 84] = 3 ## Gs\n", "\n", "print seqint[:, :10]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fill the invariants counts array\n", "This records the number of times each of the patterns in this array are observed. We will also fill two \"alternate arrays\", representing the site patterns observed for other two possible resolutions of our quartet. This first array assumes taxa are ordered as 12|34. We can fill the two alternates (13|24 and 14|23) from the patterns observed here. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[30 0 0 0 0 10 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 10 0 0 0 0 20 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 10 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 10 0 0 0 0 0 0 0 0 0 30]]\n" ] } ], "source": [ "## create a 16x16 array of zeros\n", "invcounts = np.zeros((16,16), dtype=int)\n", "\n", "## fill the array according to a set of rules\n", "for idx in xrange(seqalign.shape[1]):\n", " i = seqint[:, idx] \n", " invcounts[(4*i[0])+i[1], (4*i[2])+i[3]] += 1\n", " \n", "print invcounts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fill the three alternative matrices with site counts" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[[30 0 0 0 0 10 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 10 0 0 0 0 20 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 10 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 10 0 0 0 0 0 0 0 0 0 30]]\n", "\n", " [[30 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 10 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 20 0 0 0 10 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 10 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 30]]\n", "\n", " [[30 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 10 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 20 0 0 0 10 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 10 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 30]]]\n" ] } ], "source": [ "## Let's create three arrays \n", "mats_ints = np.zeros((3, 16, 16), dtype=int)\n", "mats_ints[0] = invcounts\n", "\n", "## Fill the alternates\n", "x = 0\n", "for i in xrange(0, 16, 4):\n", " for j in xrange(0, 16, 4):\n", " mats_ints[1, i:i+4, j:j+4] = mats_ints[0, x].reshape(4, 4)\n", " mats_ints[2, i:i+4, j:j+4] = mats_ints[0, x].reshape(4, 4).T\n", " x += 1\n", "\n", "print mats_ints[:, :, :]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### What does this look like as strings?\n", "We print this just for our edification." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[['AAAA' 'AAAC' 'AAAG' 'AAAT' 'AACA' 'AACC' 'AACG' 'AACT' 'AAGA' 'AAGC']\n", " ['ACAA' 'ACAC' 'ACAG' 'ACAT' 'ACCA' 'ACCC' 'ACCG' 'ACCT' 'ACGA' 'ACGC']\n", " ['AGAA' 'AGAC' 'AGAG' 'AGAT' 'AGCA' 'AGCC' 'AGCG' 'AGCT' 'AGGA' 'AGGC']\n", " ['ATAA' 'ATAC' 'ATAG' 'ATAT' 'ATCA' 'ATCC' 'ATCG' 'ATCT' 'ATGA' 'ATGC']\n", " ['CAAA' 'CAAC' 'CAAG' 'CAAT' 'CACA' 'CACC' 'CACG' 'CACT' 'CAGA' 'CAGC']\n", " ['CCAA' 'CCAC' 'CCAG' 'CCAT' 'CCCA' 'CCCC' 'CCCG' 'CCCT' 'CCGA' 'CCGC']\n", " ['CGAA' 'CGAC' 'CGAG' 'CGAT' 'CGCA' 'CGCC' 'CGCG' 'CGCT' 'CGGA' 'CGGC']\n", " ['CTAA' 'CTAC' 'CTAG' 'CTAT' 'CTCA' 'CTCC' 'CTCG' 'CTCT' 'CTGA' 'CTGC']\n", " ['GAAA' 'GAAC' 'GAAG' 'GAAT' 'GACA' 'GACC' 'GACG' 'GACT' 'GAGA' 'GAGC']\n", " ['GCAA' 'GCAC' 'GCAG' 'GCAT' 'GCCA' 'GCCC' 'GCCG' 'GCCT' 'GCGA' 'GCGC']\n", " ['GGAA' 'GGAC' 'GGAG' 'GGAT' 'GGCA' 'GGCC' 'GGCG' 'GGCT' 'GGGA' 'GGGC']\n", " ['GTAA' 'GTAC' 'GTAG' 'GTAT' 'GTCA' 'GTCC' 'GTCG' 'GTCT' 'GTGA' 'GTGC']\n", " ['TAAA' 'TAAC' 'TAAG' 'TAAT' 'TACA' 'TACC' 'TACG' 'TACT' 'TAGA' 'TAGC']\n", " ['TCAA' 'TCAC' 'TCAG' 'TCAT' 'TCCA' 'TCCC' 'TCCG' 'TCCT' 'TCGA' 'TCGC']\n", " ['TGAA' 'TGAC' 'TGAG' 'TGAT' 'TGCA' 'TGCC' 'TGCG' 'TGCT' 'TGGA' 'TGGC']\n", " ['TTAA' 'TTAC' 'TTAG' 'TTAT' 'TTCA' 'TTCC' 'TTCG' 'TTCT' 'TTGA' 'TTGC']]\n", "\n", " [['AAAA' 'AAAC' 'AAAG' 'AAAT' 'ACAA' 'ACAC' 'ACAG' 'ACAT' 'AGAA' 'AGAC']\n", " ['AACA' 'AACC' 'AACG' 'AACT' 'ACCA' 'ACCC' 'ACCG' 'ACCT' 'AGCA' 'AGCC']\n", " ['AAGA' 'AAGC' 'AAGG' 'AAGT' 'ACGA' 'ACGC' 'ACGG' 'ACGT' 'AGGA' 'AGGC']\n", " ['AATA' 'AATC' 'AATG' 'AATT' 'ACTA' 'ACTC' 'ACTG' 'ACTT' 'AGTA' 'AGTC']\n", " ['CAAA' 'CAAC' 'CAAG' 'CAAT' 'CCAA' 'CCAC' 'CCAG' 'CCAT' 'CGAA' 'CGAC']\n", " ['CACA' 'CACC' 'CACG' 'CACT' 'CCCA' 'CCCC' 'CCCG' 'CCCT' 'CGCA' 'CGCC']\n", " ['CAGA' 'CAGC' 'CAGG' 'CAGT' 'CCGA' 'CCGC' 'CCGG' 'CCGT' 'CGGA' 'CGGC']\n", " ['CATA' 'CATC' 'CATG' 'CATT' 'CCTA' 'CCTC' 'CCTG' 'CCTT' 'CGTA' 'CGTC']\n", " ['GAAA' 'GAAC' 'GAAG' 'GAAT' 'GCAA' 'GCAC' 'GCAG' 'GCAT' 'GGAA' 'GGAC']\n", " ['GACA' 'GACC' 'GACG' 'GACT' 'GCCA' 'GCCC' 'GCCG' 'GCCT' 'GGCA' 'GGCC']\n", " ['GAGA' 'GAGC' 'GAGG' 'GAGT' 'GCGA' 'GCGC' 'GCGG' 'GCGT' 'GGGA' 'GGGC']\n", " ['GATA' 'GATC' 'GATG' 'GATT' 'GCTA' 'GCTC' 'GCTG' 'GCTT' 'GGTA' 'GGTC']\n", " ['TAAA' 'TAAC' 'TAAG' 'TAAT' 'TCAA' 'TCAC' 'TCAG' 'TCAT' 'TGAA' 'TGAC']\n", " ['TACA' 'TACC' 'TACG' 'TACT' 'TCCA' 'TCCC' 'TCCG' 'TCCT' 'TGCA' 'TGCC']\n", " ['TAGA' 'TAGC' 'TAGG' 'TAGT' 'TCGA' 'TCGC' 'TCGG' 'TCGT' 'TGGA' 'TGGC']\n", " ['TATA' 'TATC' 'TATG' 'TATT' 'TCTA' 'TCTC' 'TCTG' 'TCTT' 'TGTA' 'TGTC']]\n", "\n", " [['AAAA' 'AACA' 'AAGA' 'AATA' 'ACAA' 'ACCA' 'ACGA' 'ACTA' 'AGAA' 'AGCA']\n", " ['AAAC' 'AACC' 'AAGC' 'AATC' 'ACAC' 'ACCC' 'ACGC' 'ACTC' 'AGAC' 'AGCC']\n", " ['AAAG' 'AACG' 'AAGG' 'AATG' 'ACAG' 'ACCG' 'ACGG' 'ACTG' 'AGAG' 'AGCG']\n", " ['AAAT' 'AACT' 'AAGT' 'AATT' 'ACAT' 'ACCT' 'ACGT' 'ACTT' 'AGAT' 'AGCT']\n", " ['CAAA' 'CACA' 'CAGA' 'CATA' 'CCAA' 'CCCA' 'CCGA' 'CCTA' 'CGAA' 'CGCA']\n", " ['CAAC' 'CACC' 'CAGC' 'CATC' 'CCAC' 'CCCC' 'CCGC' 'CCTC' 'CGAC' 'CGCC']\n", " ['CAAG' 'CACG' 'CAGG' 'CATG' 'CCAG' 'CCCG' 'CCGG' 'CCTG' 'CGAG' 'CGCG']\n", " ['CAAT' 'CACT' 'CAGT' 'CATT' 'CCAT' 'CCCT' 'CCGT' 'CCTT' 'CGAT' 'CGCT']\n", " ['GAAA' 'GACA' 'GAGA' 'GATA' 'GCAA' 'GCCA' 'GCGA' 'GCTA' 'GGAA' 'GGCA']\n", " ['GAAC' 'GACC' 'GAGC' 'GATC' 'GCAC' 'GCCC' 'GCGC' 'GCTC' 'GGAC' 'GGCC']\n", " ['GAAG' 'GACG' 'GAGG' 'GATG' 'GCAG' 'GCCG' 'GCGG' 'GCTG' 'GGAG' 'GGCG']\n", " ['GAAT' 'GACT' 'GAGT' 'GATT' 'GCAT' 'GCCT' 'GCGT' 'GCTT' 'GGAT' 'GGCT']\n", " ['TAAA' 'TACA' 'TAGA' 'TATA' 'TCAA' 'TCCA' 'TCGA' 'TCTA' 'TGAA' 'TGCA']\n", " ['TAAC' 'TACC' 'TAGC' 'TATC' 'TCAC' 'TCCC' 'TCGC' 'TCTC' 'TGAC' 'TGCC']\n", " ['TAAG' 'TACG' 'TAGG' 'TATG' 'TCAG' 'TCCG' 'TCGG' 'TCTG' 'TGAG' 'TGCG']\n", " ['TAAT' 'TACT' 'TAGT' 'TATT' 'TCAT' 'TCCT' 'TCGT' 'TCTT' 'TGAT' 'TGCT']]]\n" ] } ], "source": [ "## What does it look like as character strings\n", "mats_strs = np.zeros((3, 16, 16), dtype=\"S4\")\n", "mats_strs[0] = pi_arr\n", "\n", "## fill alternates\n", "x = 0\n", "for i in xrange(0, 16, 4):\n", " for j in xrange(0, 16, 4):\n", " mats_strs[1, i:i+4, j:j+4] = mats_strs[0, x].reshape(4, 4)\n", " mats_strs[2, i:i+4, j:j+4] = mats_strs[0, x].reshape(4, 4).T\n", " x += 1\n", "\n", "## print just first 10 columns so that it's more readable\n", "print mats_strs[:, :, :10]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Calculate the matrix rank\n", "In this case you can clearly see that the first matrix is the better fit, it has two values that *vanish*, as we expect invariants to do, and thus the rank is 4, where it is 6 for the other two topologies. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mat=0, rank=4\n", "mat=1, rank=6\n", "mat=2, rank=6\n" ] } ], "source": [ "### Well, first we can just calculate the matrix rank\n", "### b/c if they have different rank then we just choose \n", "### the lowest one\n", "for mat in xrange(3):\n", " print \"mat={}, rank={}\".format(mat, np.linalg.matrix_rank(mats_ints[mat, :, :]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Why no SVD? \n", "We didn't need it because there were so many empty columns/rows. However, most of the time with real data, and especially large data sets, the matrices will have little or no empty columns, and so the matrices will all be of equal rank. In that case we need to use SVD to calculate which invariants approach most closely to zero from the SVD scores. \n", "\n", "Here we calculate the squared sum of SVD scores above matrix rank 4, the minimum rank observed in the data set... (not sure why this is, exactly). " ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mat=0, rank=4, sum=2.08152798678e-15\n", "mat=1, rank=6, sum=14.1421356237\n", "mat=2, rank=6, sum=14.1421356237\n" ] } ], "source": [ "for idx in xrange(3):\n", " ## if all \n", " np.linalg.matrix_rank(mats_ints[idx])\n", "\n", " ## calculate SVD\n", " scores = np.linalg.svd(mats_ints[idx])[1]\n", " \n", " ## sum square of X smallest scores\n", " print \"mat={}, rank={}, sum={}\"\\\n", " .format(idx, \n", " np.linalg.matrix_rank(mats_ints[idx]),\n", " np.sqrt(np.sum(scores[4:]**2)))\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### There we have it, the first topology (12|34) is best." ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 0 }