analyses/20191030_DMR_functional_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# intersect TE gff and gene gff with DMR files generated from 20191030 analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/Shelly/Documents/GitHub/Shelly_Pgenerosa/analyses\n"
     ]
    }
   ],
   "source": [
    "cd ~/Documents/GitHub/Shelly_Pgenerosa/analyses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "!mkdir 20191030_anno"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/Shelly/Documents/GitHub/Shelly_Pgenerosa/analyses/20191030_anno\n"
     ]
    }
   ],
   "source": [
    "cd 20191030_anno"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#sed to remove chr prefix from chromosome names in DMR bed files\n",
    "#bedtools intersect to find DMRs overlapping with genomic features\n",
    "# use -wb option to output gff info along with DMR info\n",
    "\n",
    "#for ambient samples at all time points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "FILES=/Users/Shelly/Documents/GitHub/Shelly_Pgenerosa/analyses/DMRs_heatmap/*_DMR250bp_MCmax30_cov5x_rms_results_collapsed_AOV0.1.DMR.bed\n",
    "for f in $FILES\n",
    "do\n",
    "sed 's/chrPGA/PGA/g' ${f} | \\\n",
    "intersectBed \\\n",
    "-a - \\\n",
    "-b /Volumes/web/Atumefaciens/20190928_Pgenerosa_v074.a4_gensas_annotation/Panopea-generosa-vv0.74.a4-merged-2019-10-07-4-46-46.gff3 \\\n",
    "-wb \\\n",
    "> $(basename ${f%_DMR250bp_MCmax30_cov5x_rms_results_collapsed_AOV0.1.DMR.bed})_0928GFF3.txt\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      94 amb_AllTimes_0928GFF3.txt\r\n",
      "      65 day10_AllpH_0928GFF3.txt\r\n",
      "      41 day135_AllpH_0928GFF3.txt\r\n",
      "      18 day145_AllpH_0928GFF3.txt\r\n",
      "     218 total\r\n"
     ]
    }
   ],
   "source": [
    "!wc -l *_0928GFF3.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#see what features DMRs are in "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  16 CDS\n",
      "  16 exon\n",
      "  24 gene\n",
      "  27 mRNA\n",
      "  11 repeat_region\n",
      "   9 CDS\n",
      "   9 exon\n",
      "  17 gene\n",
      "  20 mRNA\n",
      "  10 repeat_region\n",
      "   8 CDS\n",
      "   8 exon\n",
      "   9 gene\n",
      "  11 mRNA\n",
      "   5 repeat_region\n",
      "   2 CDS\n",
      "   2 exon\n",
      "   4 gene\n",
      "   5 mRNA\n",
      "   5 repeat_region\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "FILES=*_0928GFF3.txt\n",
    "for f in $FILES\n",
    "do\n",
    "awk -F\"\\t\" '{print $7}' ${f}|\\\n",
    "sort | \\\n",
    "uniq -c\n",
    "done\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PGA_scaffold10__49_contigs__length_53961475 11112810 11112940 6 gene PGEN_.00g208420\r\n",
      "PGA_scaffold10__49_contigs__length_53961475 11112810 11112940 6 mRNA PGEN_.00g208420.m01\r\n",
      "PGA_scaffold10__49_contigs__length_53961475 11112810 11112940 6 exon PGEN_.00g208420.m01.exon01\r\n",
      "PGA_scaffold10__49_contigs__length_53961475 11112810 11112940 6 CDS PGEN_.00g208420.m01.CDS01\r\n",
      "PGA_scaffold11__79_contigs__length_51449921 34242236 34242266 4 gene PGEN_.00g236130\r\n",
      "PGA_scaffold11__79_contigs__length_51449921 34242236 34242266 4 mRNA PGEN_.00g236130.m01\r\n",
      "PGA_scaffold11__79_contigs__length_51449921 46672468 46672516 4 gene PGEN_.00g240600\r\n",
      "PGA_scaffold11__79_contigs__length_51449921 46672468 46672516 4 mRNA PGEN_.00g240600.m01\r\n",
      "PGA_scaffold12__71_contigs__length_50438331 7000227 7000231 3 repeat_region 19535.GS22252516.PGEN_.repeat01127701\r\n",
      "PGA_scaffold12__71_contigs__length_50438331 7000227 7000253 3 gene PGEN_.00g245080\r\n"
     ]
    }
   ],
   "source": [
    "!head ~/Documents/GitHub/Shelly_Pgenerosa/analyses/20191030_anno/amb_AllTimes_0928GFF3.txt | \\\n",
    "awk -F\"\\t\" '{sub(/;.*/,\"\",$13); sub(\"ID=\",\"\",$13);print $1,$2,$3,$4,$7,$13}'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#create file with DMR and gene ID for each set of DMRs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "for f in *_0928GFF3.txt\n",
    "do\n",
    "cat ${f} | \\\n",
    "awk -F\"\\t\" '{sub(/;.*/,\"\",$13); sub(\"ID=\",\"\",$13);print $1,$2,$3,$4,$7,$13}' \\\n",
    "> $(basename ${f} _0928GFF3.txt).geneID.txt\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PGA_scaffold10__49_contigs__length_53961475 11112810 11112940 6 PGEN_.00g208420\r\n",
      "PGA_scaffold10__49_contigs__length_53961475 11112810 11112940 6 PGEN_.00g208420.m01\r\n",
      "PGA_scaffold10__49_contigs__length_53961475 11112810 11112940 6 PGEN_.00g208420.m01.exon01\r\n",
      "PGA_scaffold10__49_contigs__length_53961475 11112810 11112940 6 PGEN_.00g208420.m01.CDS01\r\n",
      "PGA_scaffold11__79_contigs__length_51449921 34242236 34242266 4 PGEN_.00g236130\r\n",
      "PGA_scaffold11__79_contigs__length_51449921 34242236 34242266 4 PGEN_.00g236130.m01\r\n",
      "PGA_scaffold11__79_contigs__length_51449921 46672468 46672516 4 PGEN_.00g240600\r\n",
      "PGA_scaffold11__79_contigs__length_51449921 46672468 46672516 4 PGEN_.00g240600.m01\r\n",
      "PGA_scaffold12__71_contigs__length_50438331 7000227 7000231 3 19535.GS22252516.PGEN_.repeat01127701\r\n",
      "PGA_scaffold12__71_contigs__length_50438331 7000227 7000253 3 PGEN_.00g245080\r\n"
     ]
    }
   ],
   "source": [
    "!head amb_AllTimes.geneID.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#create file with DMR and sp ID for each set of DMRs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "for f in *_0928GFF3.txt\n",
    "do\n",
    "cat ${f} | \\\n",
    "awk -F\"\\t\" '{sub(/.*;Notes=/,\"\",$13);sub(/sp\\|/,\"sp_\",$13);sub(/\\|.*/,\"\",$13);print $1,$2,$3,$4,$7,$13}'| \\\n",
    "grep \"sp_\" | sed 's/sp_//g' \\\n",
    "> $(basename ${f} _0928GFF3.txt).spID.txt\n",
    "done\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PGA_scaffold11__79_contigs__length_51449921 34242236 34242266 4 gene P78333\r\n",
      "PGA_scaffold11__79_contigs__length_51449921 46672468 46672516 4 gene Q3MHL6\r\n",
      "PGA_scaffold12__71_contigs__length_50438331 7000227 7000253 3 gene Q6PCR7\r\n",
      "PGA_scaffold12__71_contigs__length_50438331 44762795 44762874 6 gene Q5F3G0\r\n",
      "PGA_scaffold14__91_contigs__length_45393038 10741127 10741221 6 gene R9UTQ8\r\n",
      "PGA_scaffold14__91_contigs__length_45393038 15698493 15698502 4 gene E9Q309\r\n",
      "PGA_scaffold17__51_contigs__length_34923512 1185754 1185791 5 gene Q92851\r\n",
      "PGA_scaffold17__51_contigs__length_34923512 10219299 10219326 3 gene O60229\r\n",
      "PGA_scaffold17__51_contigs__length_34923512 11292320 11292361 4 gene A6QP79\r\n",
      "PGA_scaffold18__69_contigs__length_27737463 18578165 18578210 7 gene Q6A098\r\n"
     ]
    }
   ],
   "source": [
    "!head amb_AllTimes.spID.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#join gene ID and sp ID files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "for f in *.geneID.txt\n",
    "do\n",
    "cat ${f} | \\\n",
    "awk 'NR==FNR{a[$1,$2,$3,$4,$5]=$6;next}{if(($1,$2,$3,$4,$5)in a)print $0,a[$1,$2,$3,$4,$5];else print $0,\"NA\"}' \\\n",
    "$(basename ${f} .geneID.txt).spID.txt - \\\n",
    "> $(basename ${f} .geneID.txt).geneANDspID.txt\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PGA_scaffold10__49_contigs__length_53961475 11112810 11112940 6 gene PGEN_.00g208420 NA\r\n",
      "PGA_scaffold10__49_contigs__length_53961475 11112810 11112940 6 mRNA PGEN_.00g208420.m01 NA\r\n",
      "PGA_scaffold10__49_contigs__length_53961475 11112810 11112940 6 exon PGEN_.00g208420.m01.exon01 NA\r\n",
      "PGA_scaffold10__49_contigs__length_53961475 11112810 11112940 6 CDS PGEN_.00g208420.m01.CDS01 NA\r\n",
      "PGA_scaffold11__79_contigs__length_51449921 34242236 34242266 4 gene PGEN_.00g236130 P78333\r\n",
      "PGA_scaffold11__79_contigs__length_51449921 34242236 34242266 4 mRNA PGEN_.00g236130.m01 NA\r\n",
      "PGA_scaffold11__79_contigs__length_51449921 46672468 46672516 4 gene PGEN_.00g240600 Q3MHL6\r\n",
      "PGA_scaffold11__79_contigs__length_51449921 46672468 46672516 4 mRNA PGEN_.00g240600.m01 NA\r\n",
      "PGA_scaffold12__71_contigs__length_50438331 7000227 7000231 3 repeat_region 19535.GS22252516.PGEN_.repeat01127701 NA\r\n",
      "PGA_scaffold12__71_contigs__length_50438331 7000227 7000253 3 gene PGEN_.00g245080 Q6PCR7\r\n"
     ]
    }
   ],
   "source": [
    "!head amb_AllTimes.geneANDspID.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#count number of lines in all .geneANDspID.txt files to \n",
    "# verify they match  number of lines in _0928GFF3.txt files above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      94 amb_AllTimes.geneANDspID.txt\r\n",
      "      65 day10_AllpH.geneANDspID.txt\r\n",
      "      41 day135_AllpH.geneANDspID.txt\r\n",
      "      18 day145_AllpH.geneANDspID.txt\r\n",
      "     218 total\r\n"
     ]
    }
   ],
   "source": [
    "!wc -l *.geneANDspID.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#now to get the GO IDs for each SP ID I will use the Uniprot mapping function here:\n",
    "# https://www.uniprot.org/uploadlists/\n",
    "#saved file as \"uniprot-20191001_GO.tab\"\n",
    "\n",
    "#first take a look at the DMR/GFF column for matching to the uniprot GO data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A0A0R4IBK5\r\n",
      "A6H5Z3\r\n",
      "A6QP79\r\n",
      "B5X1P9\r\n",
      "E9Q309\r\n",
      "O35274\r\n",
      "O43374\r\n",
      "O60229\r\n",
      "P04052\r\n",
      "P23780\r\n",
      "P23919\r\n",
      "P46023\r\n",
      "P78333\r\n",
      "P9WQE8\r\n",
      "Q3MHL6\r\n",
      "Q5F3G0\r\n",
      "Q6A098\r\n",
      "Q6GNV7\r\n",
      "Q6PCR7\r\n",
      "Q80VI1\r\n",
      "Q8WUA2\r\n",
      "Q92851\r\n",
      "Q9NZJ4\r\n",
      "Q9UJS0\r\n",
      "R9UTQ8\r\n"
     ]
    }
   ],
   "source": [
    "!cat *.geneANDspID.txt | awk '{if($7!=\"NA\")print $7}' | sort | uniq "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      25\r\n"
     ]
    }
   ],
   "source": [
    "!cat *.geneANDspID.txt | awk '{if($7!=\"NA\")print $7}' | sort | uniq | wc -l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   1 A0A0R4IBK5\r\n",
      "   1 A6H5Z3\r\n",
      "   2 A6QP79\r\n",
      "   1 B5X1P9\r\n",
      "   1 E9Q309\r\n",
      "   1 O35274\r\n",
      "   2 O43374\r\n",
      "   1 O60229\r\n",
      "   1 P04052\r\n",
      "   3 P23780\r\n",
      "   1 P23919\r\n",
      "   1 P46023\r\n",
      "   1 P78333\r\n",
      "   1 P9WQE8\r\n",
      "   1 Q3MHL6\r\n",
      "   1 Q5F3G0\r\n",
      "   1 Q6A098\r\n",
      "   2 Q6GNV7\r\n",
      "   1 Q6PCR7\r\n",
      "   1 Q80VI1\r\n",
      "   1 Q8WUA2\r\n",
      "   1 Q92851\r\n",
      "   1 Q9NZJ4\r\n",
      "   1 Q9UJS0\r\n",
      "   1 R9UTQ8\r\n"
     ]
    }
   ],
   "source": [
    "!cat *.geneANDspID.txt | awk '{if($7!=\"NA\")print $7}' | sort | uniq -c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#next take a look at the Uniprot GO table I downloaded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yourlist:M201910316746803381A1F0E0DB47453E0216320D0AA5BAI\tEntry\tEntry name\tProtein names\tGene ontology IDs\tGene ontology (GO)\r\n",
      "A0A0R4IBK5\tA0A0R4IBK5\tR213A_DANRE\tE3 ubiquitin-protein ligase rnf213-alpha (EC 2.3.2.27) (EC 3.6.4.-) (Mysterin-A) (Mysterin-alpha) (RING finger protein 213-A) (RING finger protein 213-alpha) (RING-type E3 ubiquitin transferase rnf213-alpha)\tGO:0002040; GO:0004842; GO:0005829; GO:0006511; GO:0008015; GO:0016887; GO:0046872\tcytosol [GO:0005829]; ATPase activity [GO:0016887]; metal ion binding [GO:0046872]; ubiquitin-protein transferase activity [GO:0004842]; blood circulation [GO:0008015]; sprouting angiogenesis [GO:0002040]; ubiquitin-dependent protein catabolic process [GO:0006511]\r\n",
      "A6H5Z3\tA6H5Z3\tEXC6B_MOUSE\tExocyst complex component 6B (Exocyst complex component Sec15B) (SEC15-like protein 2)\tGO:0000145; GO:0006887; GO:0006893; GO:0006904\texocyst [GO:0000145]; exocytosis [GO:0006887]; Golgi to plasma membrane transport [GO:0006893]; vesicle docking involved in exocytosis [GO:0006904]\r\n",
      "A6QP79\tA6QP79\tCOL12_BOVIN\tCollectin-12 (Collectin placenta protein 1) (CL-P1)\tGO:0004888; GO:0005581; GO:0006910; GO:0006955; GO:0008329; GO:0016021; GO:0030169; GO:0030246; GO:0034138; GO:0046872; GO:0060355; GO:0071360\tcollagen trimer [GO:0005581]; integral component of membrane [GO:0016021]; carbohydrate binding [GO:0030246]; low-density lipoprotein particle binding [GO:0030169]; metal ion binding [GO:0046872]; signaling pattern recognition receptor activity [GO:0008329]; transmembrane signaling receptor activity [GO:0004888]; cellular response to exogenous dsRNA [GO:0071360]; immune response [GO:0006955]; phagocytosis, recognition [GO:0006910]; positive regulation of cell adhesion molecule production [GO:0060355]; toll-like receptor 3 signaling pathway [GO:0034138]\r\n",
      "B5X1P9\tB5X1P9\tABRX1_SALSA\tBRCA1-A complex subunit Abraxas 1 (Coiled-coil domain-containing protein 98) (Protein FAM175A)\tGO:0005634; GO:0006302; GO:0006325; GO:0010212; GO:0031593; GO:0045739; GO:0070531; GO:0072425\tBRCA1-A complex [GO:0070531]; nucleus [GO:0005634]; polyubiquitin modification-dependent protein binding [GO:0031593]; chromatin organization [GO:0006325]; double-strand break repair [GO:0006302]; positive regulation of DNA repair [GO:0045739]; response to ionizing radiation [GO:0010212]; signal transduction involved in G2 DNA damage checkpoint [GO:0072425]\r\n",
      "E9Q309\tE9Q309\tCE350_MOUSE\tCentrosome-associated protein 350 (Cep350)\tGO:0005634; GO:0005813; GO:0005814; GO:0005819; GO:0008017; GO:0034453; GO:0042995\tcell projection [GO:0042995]; centriole [GO:0005814]; centrosome [GO:0005813]; nucleus [GO:0005634]; spindle [GO:0005819]; microtubule binding [GO:0008017]; microtubule anchoring [GO:0034453]\r\n",
      "O35274\tO35274\tNEB2_RAT\tNeurabin-2 (Neurabin-II) (Neural tissue-specific F-actin-binding protein II) (PP1bp134) (Protein phosphatase 1 regulatory subunit 9B) (Spinophilin) (p130)\tGO:0001932; GO:0001975; GO:0003006; GO:0003779; GO:0004672; GO:0004864; GO:0005654; GO:0005737; GO:0005886; GO:0005912; GO:0007015; GO:0007568; GO:0007612; GO:0008022; GO:0008157; GO:0010243; GO:0014069; GO:0014070; GO:0015629; GO:0016358; GO:0016477; GO:0019722; GO:0019900; GO:0021766; GO:0021987; GO:0030027; GO:0030042; GO:0030175; GO:0030308; GO:0030425; GO:0030426; GO:0030864; GO:0031175; GO:0031749; GO:0032355; GO:0032515; GO:0032587; GO:0032591; GO:0034695; GO:0035094; GO:0035690; GO:0035902; GO:0042127; GO:0042493; GO:0043005; GO:0043025; GO:0043197; GO:0043200; GO:0044325; GO:0044326; GO:0044327; GO:0044877; GO:0046847; GO:0048545; GO:0050804; GO:0051015; GO:0060179; GO:0061458; GO:0071315; GO:0071364; GO:0071392; GO:0071407; GO:0097338; GO:1901653; GO:1903078; GO:1903119; GO:1903829; GO:1904372; GO:1904373; GO:1904386; GO:1990778; GO:1990780; GO:2000474\tactin cytoskeleton [GO:0015629]; adherens junction [GO:0005912]; cortical actin cytoskeleton [GO:0030864]; cytoplasm [GO:0005737]; cytoplasmic side of dendritic spine plasma membrane [GO:1990780]; dendrite [GO:0030425]; dendritic spine [GO:0043197]; dendritic spine head [GO:0044327]; dendritic spine membrane [GO:0032591]; dendritic spine neck [GO:0044326]; filopodium [GO:0030175]; growth cone [GO:0030426]; lamellipodium [GO:0030027]; neuron projection [GO:0043005]; neuronal cell body [GO:0043025]; nucleoplasm [GO:0005654]; plasma membrane [GO:0005886]; postsynaptic density [GO:0014069]; ruffle membrane [GO:0032587]; actin binding [GO:0003779]; actin filament binding [GO:0051015]; D2 dopamine receptor binding [GO:0031749]; ion channel binding [GO:0044325]; kinase binding [GO:0019900]; protein C-terminus binding [GO:0008022]; protein kinase activity [GO:0004672]; protein phosphatase 1 binding [GO:0008157]; protein phosphatase inhibitor activity [GO:0004864]; protein-containing complex binding [GO:0044877]; actin filament depolymerization [GO:0030042]; actin filament organization [GO:0007015]; aging [GO:0007568]; calcium-mediated signaling [GO:0019722]; cell migration [GO:0016477]; cellular response to drug [GO:0035690]; cellular response to epidermal growth factor stimulus [GO:0071364]; cellular response to estradiol stimulus [GO:0071392]; cellular response to morphine [GO:0071315]; cellular response to organic cyclic compound [GO:0071407]; cellular response to peptide [GO:1901653]; cerebral cortex development [GO:0021987]; dendrite development [GO:0016358]; developmental process involved in reproduction [GO:0003006]; filopodium assembly [GO:0046847]; hippocampus development [GO:0021766]; learning [GO:0007612]; male mating behavior [GO:0060179]; modulation of chemical synaptic transmission [GO:0050804]; negative regulation of cell growth [GO:0030308]; negative regulation of phosphoprotein phosphatase activity [GO:0032515]; neuron projection development [GO:0031175]; positive regulation of cellular protein localization [GO:1903829]; positive regulation of protein localization to actin cortical patch [GO:1904372]; positive regulation of protein localization to plasma membrane [GO:1903078]; protein localization to actin cytoskeleton [GO:1903119]; protein localization to cell periphery [GO:1990778]; regulation of cell population proliferation [GO:0042127]; regulation of opioid receptor signaling pathway [GO:2000474]; regulation of protein phosphorylation [GO:0001932]; reproductive system development [GO:0061458]; response to amino acid [GO:0043200]; response to amphetamine [GO:0001975]; response to clozapine [GO:0097338]; response to drug [GO:0042493]; response to estradiol [GO:0032355]; response to immobilization stress [GO:0035902]; response to kainic acid [GO:1904373]; response to L-phenylalanine derivative [GO:1904386]; response to nicotine [GO:0035094]; response to organic cyclic compound [GO:0014070]; response to organonitrogen compound [GO:0010243]; response to prostaglandin E [GO:0034695]; response to steroid hormone [GO:0048545]\r\n",
      "O43374\tO43374\tRASL2_HUMAN\tRas GTPase-activating protein 4 (Calcium-promoted Ras inactivator) (Ras p21 protein activator 4) (RasGAP-activating-like protein 2)\tGO:0000165; GO:0005096; GO:0005543; GO:0005829; GO:0005886; GO:0034260; GO:0046580; GO:0046872; GO:0071277\tcytosol [GO:0005829]; plasma membrane [GO:0005886]; GTPase activator activity [GO:0005096]; metal ion binding [GO:0046872]; phospholipid binding [GO:0005543]; cellular response to calcium ion [GO:0071277]; MAPK cascade [GO:0000165]; negative regulation of GTPase activity [GO:0034260]; negative regulation of Ras protein signal transduction [GO:0046580]\r\n",
      "O60229\tO60229\tKALRN_HUMAN\tKalirin (EC 2.7.11.1) (Huntingtin-associated protein-interacting protein) (Protein Duo) (Serine/threonine-protein kinase with Dbl- and pleckstrin homology domain)\tGO:0004674; GO:0005085; GO:0005089; GO:0005524; GO:0005654; GO:0005829; GO:0006468; GO:0007165; GO:0007186; GO:0007399; GO:0007528; GO:0007595; GO:0007613; GO:0008344; GO:0015629; GO:0016192; GO:0035023; GO:0035176; GO:0035556; GO:0042711; GO:0043065; GO:0043547; GO:0046872; GO:0046959; GO:0048013; GO:0051056; GO:0060125; GO:0060137; GO:0061003; GO:0070062\tactin cytoskeleton [GO:0015629]; cytosol [GO:0005829]; extracellular exosome [GO:0070062]; nucleoplasm [GO:0005654]; ATP binding [GO:0005524]; guanyl-nucleotide exchange factor activity [GO:0005085]; metal ion binding [GO:0046872]; protein serine/threonine kinase activity [GO:0004674]; Rho guanyl-nucleotide exchange factor activity [GO:0005089]; adult locomotory behavior [GO:0008344]; ephrin receptor signaling pathway [GO:0048013]; G protein-coupled receptor signaling pathway [GO:0007186]; habituation [GO:0046959]; intracellular signal transduction [GO:0035556]; lactation [GO:0007595]; maternal behavior [GO:0042711]; maternal process involved in parturition [GO:0060137]; memory [GO:0007613]; negative regulation of growth hormone secretion [GO:0060125]; nervous system development [GO:0007399]; neuromuscular junction development [GO:0007528]; positive regulation of apoptotic process [GO:0043065]; positive regulation of dendritic spine morphogenesis [GO:0061003]; positive regulation of GTPase activity [GO:0043547]; protein phosphorylation [GO:0006468]; regulation of Rho protein signal transduction [GO:0035023]; regulation of small GTPase mediated signal transduction [GO:0051056]; signal transduction [GO:0007165]; social behavior [GO:0035176]; vesicle-mediated transport [GO:0016192]\r\n",
      "P04052\tP04052\tRPB1_DROME\tDNA-directed RNA polymerase II subunit RPB1 (RNA polymerase II subunit B1) (EC 2.7.7.6) (DNA-directed RNA polymerase III largest subunit)\tGO:0003677; GO:0003899; GO:0005634; GO:0005665; GO:0005700; GO:0005703; GO:0006366; GO:0046872\tnucleus [GO:0005634]; polytene chromosome [GO:0005700]; polytene chromosome puff [GO:0005703]; RNA polymerase II, core complex [GO:0005665]; DNA binding [GO:0003677]; DNA-directed 5'-3' RNA polymerase activity [GO:0003899]; metal ion binding [GO:0046872]; transcription by RNA polymerase II [GO:0006366]\r\n"
     ]
    }
   ],
   "source": [
    "!head uniprot-20191031_GO.tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Protein names\r\n",
      "E3 ubiquitin-protein ligase rnf213-alpha (EC 2.3.2.27) (EC 3.6.4.-) (Mysterin-A) (Mysterin-alpha) (RING finger protein 213-A) (RING finger protein 213-alpha) (RING-type E3 ubiquitin transferase rnf213-alpha)\r\n",
      "Exocyst complex component 6B (Exocyst complex component Sec15B) (SEC15-like protein 2)\r\n",
      "Collectin-12 (Collectin placenta protein 1) (CL-P1)\r\n",
      "BRCA1-A complex subunit Abraxas 1 (Coiled-coil domain-containing protein 98) (Protein FAM175A)\r\n",
      "Centrosome-associated protein 350 (Cep350)\r\n",
      "Neurabin-2 (Neurabin-II) (Neural tissue-specific F-actin-binding protein II) (PP1bp134) (Protein phosphatase 1 regulatory subunit 9B) (Spinophilin) (p130)\r\n",
      "Ras GTPase-activating protein 4 (Calcium-promoted Ras inactivator) (Ras p21 protein activator 4) (RasGAP-activating-like protein 2)\r\n",
      "Kalirin (EC 2.7.11.1) (Huntingtin-associated protein-interacting protein) (Protein Duo) (Serine/threonine-protein kinase with Dbl- and pleckstrin homology domain)\r\n",
      "DNA-directed RNA polymerase II subunit RPB1 (RNA polymerase II subunit B1) (EC 2.7.7.6) (DNA-directed RNA polymerase III largest subunit)\r\n",
      "Beta-galactosidase (EC 3.2.1.23) (Acid beta-galactosidase) (Lactase)\r\n",
      "Thymidylate kinase (EC 2.7.4.9) (dTMP kinase)\r\n",
      "G-protein coupled receptor GRL101\r\n",
      "Glypican-5 [Cleaved into: Secreted glypican-5]\r\n",
      "Phthioceranic/hydroxyphthioceranic acid synthase (EC 2.3.1.-) (Polyketide synthase pks2)\r\n",
      "TSC22 domain family protein 1\r\n",
      "Rho GTPase-activating protein 19 (Rho-type GTPase-activating protein 19)\r\n",
      "Selenocysteine insertion sequence-binding protein 2-like (SECIS-binding protein 2-like)\r\n",
      "Solute carrier family 49 member 4 homolog (Disrupted in renal carcinoma protein 2 homolog)\r\n",
      "Eukaryotic translation initiation factor 3 subunit A (eIF3a) (Eukaryotic translation initiation factor 3 subunit 10) (eIF-3-theta)\r\n",
      "E3 ubiquitin-protein ligase TRIM56 (EC 2.3.2.27) (RING-type E3 ubiquitin transferase TRIM56) (Tripartite motif-containing protein 56)\r\n",
      "Peptidyl-prolyl cis-trans isomerase-like 4 (PPIase) (EC 5.2.1.8) (Cyclophilin-like protein PPIL4) (Rotamase PPIL4)\r\n",
      "Caspase-10 (CASP-10) (EC 3.4.22.63) (Apoptotic protease Mch-4) (FAS-associated death domain protein interleukin-1B-converting enzyme 2) (FLICE2) (ICE-like apoptotic protease 4) [Cleaved into: Caspase-10 subunit p23/17; Caspase-10 subunit p12]\r\n",
      "Sacsin (DnaJ homolog subfamily C member 29) (DNAJC29)\r\n",
      "Calcium-binding mitochondrial carrier protein Aralar2 (Citrin) (Mitochondrial aspartate glutamate carrier 2) (Solute carrier family 25 member 13)\r\n",
      "L-proline trans-4-hydroxylase (P4H) (EC 1.14.11.57)\r\n"
     ]
    }
   ],
   "source": [
    "!awk -F\"\\t\" '{print $4}' uniprot-20191031_GO.tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#now bind .geneANDspID.txt files with uniprot GO annotations\n",
    "#code explanation:\n",
    "#loop through each .geneANDspID.txt file (there is one file for each DMR comparison: e.g. Day10 all pHs)\n",
    "#tr opens the file and replaces all spaces in the .geneANDspID.txt file with tabs\n",
    "#awk -F\"\\t\" specifies input files are tab delimited \n",
    "#'NR==FNR' specifies comparing two files\n",
    "#{a[$1]=$2\"\\t\"$3\"\\t\"$4\"\\t\"$5;next} creates an array of the first input file where\n",
    "## the first field ($1) is the key and the 2nd, 3rd, 4th, and 5th field are the definition\n",
    "#{if($6 in a)print $0\"\\t\"a[$6];else print $0\"\\t\"\"NA\"\"\\t\"\"NA\"\"\\t\"\"NA\"\"\\t\"\"NA\"} \n",
    "## if the 6th field of the second file is in the array (which is the list of Entry IDs in the first file)\n",
    "#print the whole like of the second file followed by the definition of the Entry ID in the array\n",
    "# if the 6th field of the second file is NOT in the array (print NA in all fields)\n",
    "# uniprot-20191001_GO.tab - \\\n",
    "# the first file is uniprot-20191001_GO.tab and the second file is the .geneANDspID.txt file "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for f in *.geneANDspID.txt\n",
    "do\n",
    "tr < ${f} ' ' '\\t' |\\\n",
    "awk -F\"\\t\" 'NR==FNR{a[$2]=$3\"\\t\"$4\"\\t\"$5\"\\t\"$6;next}{if($7 in a)print $0\"\\t\"a[$7];else print $0\"\\t\"\"NA\"\"\\t\"\"NA\"\"\\t\"\"NA\"\"\\t\"\"NA\"}' \\\n",
    "uniprot-20191031_GO.tab - \\\n",
    "> $(basename ${f} .geneANDspID.txt).GO.txt\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PGA_scaffold10__49_contigs__length_53961475\t11112810\t11112940\t6\tgene\tPGEN_.00g208420\tNA\tNA\tNA\tNA\tNA\r\n",
      "PGA_scaffold10__49_contigs__length_53961475\t11112810\t11112940\t6\tmRNA\tPGEN_.00g208420.m01\tNA\tNA\tNA\tNA\tNA\r\n",
      "PGA_scaffold10__49_contigs__length_53961475\t11112810\t11112940\t6\texon\tPGEN_.00g208420.m01.exon01\tNA\tNA\tNA\tNA\tNA\r\n",
      "PGA_scaffold10__49_contigs__length_53961475\t11112810\t11112940\t6\tCDS\tPGEN_.00g208420.m01.CDS01\tNA\tNA\tNA\tNA\tNA\r\n",
      "PGA_scaffold11__79_contigs__length_51449921\t34242236\t34242266\t4\tgene\tPGEN_.00g236130\tP78333\tGPC5_HUMAN\tGlypican-5 [Cleaved into: Secreted glypican-5]\tGO:0001523; GO:0005576; GO:0005615; GO:0005796; GO:0005886; GO:0006024; GO:0006027; GO:0009986; GO:0016021; GO:0016477; GO:0043202; GO:0046658; GO:0062023; GO:0090263; GO:1905475\tanchored component of plasma membrane [GO:0046658]; cell surface [GO:0009986]; collagen-containing extracellular matrix [GO:0062023]; extracellular region [GO:0005576]; extracellular space [GO:0005615]; Golgi lumen [GO:0005796]; integral component of membrane [GO:0016021]; lysosomal lumen [GO:0043202]; plasma membrane [GO:0005886]; cell migration [GO:0016477]; glycosaminoglycan biosynthetic process [GO:0006024]; glycosaminoglycan catabolic process [GO:0006027]; positive regulation of canonical Wnt signaling pathway [GO:0090263]; regulation of protein localization to membrane [GO:1905475]; retinoid metabolic process [GO:0001523]\r\n",
      "PGA_scaffold11__79_contigs__length_51449921\t34242236\t34242266\t4\tmRNA\tPGEN_.00g236130.m01\tNA\tNA\tNA\tNA\tNA\r\n",
      "PGA_scaffold11__79_contigs__length_51449921\t46672468\t46672516\t4\tgene\tPGEN_.00g240600\tQ3MHL6\tT22D1_BOVIN\tTSC22 domain family protein 1\tGO:0003700; GO:0005634; GO:0005737\tcytoplasm [GO:0005737]; nucleus [GO:0005634]; DNA-binding transcription factor activity [GO:0003700]\r\n",
      "PGA_scaffold11__79_contigs__length_51449921\t46672468\t46672516\t4\tmRNA\tPGEN_.00g240600.m01\tNA\tNA\tNA\tNA\tNA\r\n",
      "PGA_scaffold12__71_contigs__length_50438331\t7000227\t7000231\t3\trepeat_region\t19535.GS22252516.PGEN_.repeat01127701\tNA\tNA\tNA\tNA\tNA\r\n",
      "PGA_scaffold12__71_contigs__length_50438331\t7000227\t7000253\t3\tgene\tPGEN_.00g245080\tQ6PCR7\tEIF3A_DANRE\tEukaryotic translation initiation factor 3 subunit A (eIF3a) (Eukaryotic translation initiation factor 3 subunit 10) (eIF-3-theta)\tGO:0001732; GO:0002188; GO:0003729; GO:0003743; GO:0005852; GO:0016282; GO:0033290; GO:0043614; GO:0071540; GO:0071541\teukaryotic 43S preinitiation complex [GO:0016282]; eukaryotic 48S preinitiation complex [GO:0033290]; eukaryotic translation initiation factor 3 complex [GO:0005852]; eukaryotic translation initiation factor 3 complex, eIF3e [GO:0071540]; eukaryotic translation initiation factor 3 complex, eIF3m [GO:0071541]; multi-eIF complex [GO:0043614]; mRNA binding [GO:0003729]; translation initiation factor activity [GO:0003743]; formation of cytoplasmic translation initiation complex [GO:0001732]; translation reinitiation [GO:0002188]\r\n"
     ]
    }
   ],
   "source": [
    "!head amb_AllTimes.GO.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#read GO.txt files into R and make some plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "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": "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}