tree_building_workflow_Guide-v4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Plan\n",
    "\n",
    "0. Make sure you have all tools installed. mafft, trimal, iqtree at the very least.\n",
    "1. Get a set of sequences to build a tree of.\n",
    "2. Align these to each other.\n",
    "4. Trim the alignment, removing gaps.\n",
    "3. Build a fast phylogentic tree.\n",
    "4. Build a thorough phylogenetic tree.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Composing your fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at what we have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ".\n",
      "├── analyses\n",
      "│   ├── 1kP_LAR_orthogroup_random700_guidev3_fasttrees\n",
      "│   │   └── aligned-mafft-linsi_trim-gt6-seq80\n",
      "│   │       ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_fasttree.log\n",
      "│   │       ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_fasttree.stderr\n",
      "│   │       └── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_fasttree.tree\n",
      "│   └── 1kP_LAR_orthogroup_random700_guidev3_trees\n",
      "│       └── aligned-mafft-linsi_trim-gt6-seq80\n",
      "│           ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000\n",
      "│           ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000.bionj\n",
      "│           ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000.ckp.gz\n",
      "│           ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000.contree\n",
      "│           ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000.iqtree\n",
      "│           ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000.log\n",
      "│           ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000.mldist\n",
      "│           ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000.model.gz\n",
      "│           ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000.splits.nex\n",
      "│           ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000.treefile\n",
      "│           └── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000.uniqueseq.phy\n",
      "├── data\n",
      "│   ├── 1kP_LAR_orthogroup.fasta\n",
      "│   ├── 1kP_LAR_orthogroup_random700.fasta\n",
      "│   ├── 1kP_LAR_orthogroup_random700_guidev3.fasta\n",
      "│   ├── 1kP_LAR_orthogroup_random700.txt\n",
      "│   ├── alignments_raw\n",
      "│   │   └── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi.fasta\n",
      "│   ├── alignments_trimmed\n",
      "│   │   ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6.html\n",
      "│   │   ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80.fasta\n",
      "│   │   ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80.html\n",
      "│   │   ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq90.html\n",
      "│   │   ├── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt7.html\n",
      "│   │   └── 1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt8.html\n",
      "│   └── Erbils_guide_v3.fasta\n",
      "├── tree_building_workflow_Guide-v4.ipynb\n",
      "└── tree_building_workflow.ipynb\n",
      "\n",
      "8 directories, 28 files\n"
     ]
    }
   ],
   "source": [
    "tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "inseq=1kP_LAR_orthogroup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">KTWL-2000328-Kaliphora_madagascariensis\n",
      "GGTGYIGKFLVEASAKAGHPTFAFAREITVSDPLKGKLVESFKNSGVTLVHGDLYDHSSLVKAIKQVDVVISTVGQGQLNDQGKIIDAIKEAGNVKKFYPSEFGVDVDRQHAVEPAKFTFTGKAQFRRAIEAAGIPHTFLVSNAFAGYVLPTLGQPGVSAPPHDKVVILGDGNAKAVVNEEHDIGTFTIKTVDDPRTLNKIVYIRPPKNTVSFNEIVAIWERTIGKTLEKEYVSEEQLLKQIKEAPPPLNYMLSIYHAYLVKGDSANFAIEPSGVEASELYPDVNYTTVEEYI\n",
      ">QSNJ-2064790-Taiwania_cryptomerioides\n",
      "KSKVLIFGATGYMGHFIAIASLDSDHPTFVFIRPSSLEDPARAQILSNLRSSGIHFITGSLEDHESLVRAIQGMDVVISVLGGGQIMEQLKIIDAIKQVGTVKRFLPSEFGHDVDRAEPVEPALGFYESKRRIRRATEEAKIGYTYICCNSIAGWPYHYHTHPSKMFPPTDKIHIYGNGNIKAYFMTGKDIGSYTMRAIEDPRTLNKKLHFRPADNFLTMNELAVMWQKQLGKKLPMAVISEEDLLKCARANVMPASIVAALTHDIFIKGCQYNFSIEGDEDEEACQLYPDIKYTKLEDYFSSF\n",
      ">URCP-2006931-Athyrium_filix_femina\n",
      "MASKSKVLVIGATGYIGQHIAQAGPAQGHPTFVLVRPSSLSSKATLIHSFKLAGITILEGSLDDYESLVAAIKQVDVVISALGSNTLDQLKIVEAIKEVGTIKRFFPSEFGNDVDRVKALEPASQVFGLKVQVRRAIEEAGIPYTYVVSNGFAGYFLRSLLQPGLQSPPRDKVSIYGSGDVKFITASEEDIGAYTIKAVDDPRTLNKTLHIRPAANIVTINSLVELWEKLIGKTLEKVTVTAEDIVKQIQETPFPNNIFPAILHDIVIQGDHCNFELGPKDVEATTLYPDHKYTTSEEYLAKFV\n",
      ">ROLB-2007609-Oenothera_elata-multiple_tissue_pool\n",
      "NKQRFLPSEFGNDVDRVNAVGPAKLAFGIKANIRRAIETEGIPYTYVVNNFFAGYFLPTLLQPGATAPPRDQVVIPGNGNPKAIFNKEEDIGTYTIKAVDDPRTLNKTYTIKAVDDPRTLNKILYLRPPKNIYSFNELVALWENKIGHTLKKIYVLEDQILKDIQEASDGRFRSLTMNAFQSINHSVFVKGDQTNFDIEPAFGVEASQ\n",
      ">RFSD-2005689-RFSD-Lactuca_graminifolia-2_samples_combined\n",
      "EEEIALTTIKTVDDRRTLNKALIFRPPGNTLSFNEIVSIWESKIGKTLVKTYVSAEQLLKNIQEAPFGLSIELSVLHSVLVNGSSTNFDIEPSFGVEASELYPEVKYTTIDEYLT\n"
     ]
    }
   ],
   "source": [
    "head data/$inseq.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7449\n"
     ]
    }
   ],
   "source": [
    "grep '>' data/$inseq.fasta -c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shuf: write error: Broken pipe\n",
      "shuf: write error\n"
     ]
    }
   ],
   "source": [
    "n=700\n",
    "grep '>' data/$inseq.fasta | tr -d '>' | shuf | head -n $n > data/\"$inseq\"_random\"$n\".txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "grep '>'  ../LAR_tree/alignments_raw/LAR_orthogroup_LARselectionExtraFernsAndLycophytes_guide2_aligned-mafft-linsi.fasta | tr -d '>' | tail -n +49 > data/\"$inseq\"_random\"$n\".txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "698 ./data/1kP_LAR_orthogroup_random700.txt\n"
     ]
    }
   ],
   "source": [
    "wc -l ./data/\"$inseq\"_random700.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's add all Azolla stuff in there, but first see what is there already"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LAR-like_Azolla_filiculoides\n",
      "CVEG-2005710-Azolla_cf_caroliniana\n",
      "CVEG-2016335-Azolla_cf_caroliniana\n",
      "CVEG-2016336-Azolla_cf_caroliniana\n",
      "CVEG-2005709-Azolla_cf_caroliniana\n",
      "AZOLLAFILICULOIDES\n"
     ]
    }
   ],
   "source": [
    "grep -i azolla data/\"$inseq\"_random\"$n\".txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">CVEG-2005710-Azolla_cf._caroliniana\n",
      ">CVEG-2016335-Azolla_cf._caroliniana\n",
      ">CVEG-2016336-Azolla_cf._caroliniana\n",
      ">CVEG-2005709-Azolla_cf._caroliniana\n",
      ">AZOLLAFILICULOIDES\n"
     ]
    }
   ],
   "source": [
    "grep -i azolla data/\"$inseq\".fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove double guide sequence AZOLLAFILICULOIDES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "sed -i '/AZOLLAFILICULOIDES/d' data/\"$inseq\".fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Klefl_v1_0-kfl00402_0160-Klebsormidium_flaccidum\n",
      "AKXB-2007298-Phaeomegaceros_coriaceus\n",
      "AKXB-2067781-Phaeomegaceros_coriaceus\n",
      "JKKI-2057061-Lobomonas_rostrata\n",
      "Klefl_v1_0-kfl00034_0150-Klebsormidium_flaccidum\n",
      "AKXB-2064674-Phaeomegaceros_coriaceus\n",
      "RCBT-2175720-Sphagnum_palustre\n",
      "RSOF-2008002-Glaucosphaera_vacuolata\n",
      "AKXB-2001848-Phaeomegaceros_coriaceus\n",
      "BAZF-2006540-Chaetopeltis_orbicularis\n"
     ]
    }
   ],
   "source": [
    "head data/\"$inseq\"_random\"$n\".txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "/opt/rnaseq/bin/extract_sequence_by_id.pl --ID=<(sort data/\"$inseq\"_random\"$n\".txt | uniq) \\\n",
    "                                          --fileformat=txt \\\n",
    "                                          --sequences=<(cat data/$inseq.fasta | tr -d '.')\\\n",
    "                                          --outfile=data/\"$inseq\"_random\"$n\".fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "666\n"
     ]
    }
   ],
   "source": [
    "grep '>' -c data/\"$inseq\"_random\"$n\".fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">CVEG-2005710-Azolla_cf_caroliniana\n",
      ">CVEG-2016335-Azolla_cf_caroliniana\n",
      ">CVEG-2016336-Azolla_cf_caroliniana\n",
      ">CVEG-2005709-Azolla_cf_caroliniana\n",
      ">AZOLLAFILICULOIDES\n"
     ]
    }
   ],
   "source": [
    "n=700\n",
    "grep -i azolla data/\"$inseq\"_random\"$n\".fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1kP_LAR_orthogroup.fasta                    alignments_raw\n",
      "1kP_LAR_orthogroup_random700.fasta          alignments_trimmed\n",
      "1kP_LAR_orthogroup_random700_guidev3.fasta  Erbils_guide_v3.fasta\n",
      "1kP_LAR_orthogroup_random700.txt\n"
     ]
    }
   ],
   "source": [
    "ls data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">YSRZ-2016241-Fouqueria_macdougalli\n",
      "MQKSKVLIVGGTGYLGRRLVKASLALEHETYVLRRPEIGVDIDKAQMLLSFKKEGAHLVPGSFDDHRSLVDAVKMVDVVIRSHQILLQLKLVDAIKEAGNVKRFLPSEFGTDPAKMENAMEPGRVTFDDKMVVRKAIEDAGIPHTYVSANCFAGYFIGGLCQPGQIFPSKDHVVLLGDGNQKAIYVDEDDIATYTIKTIDDPQTINKTLYLRPPENILSQREVIQIWEKLIGKELQKSSISSEQFLATMKEQDYAEQVGLTHYYHVCYEGCLANFEIGDE\n",
      ">HJMP-2011848-Astragalus_membranaceus\n",
      "TLNKTLLLRPSANTLTFNELVSLWENKIKSTIHKIYIPEEKLLKTIQETPFPGNLISALCHSILVKGDSANFDIDPSFCVEGTGLYPEVKYTTVDNYLNAFV\n",
      ">OSIP-2125743-Garcinia_livingstonei\n",
      "SPNILEACLGRTLIIGSTGFIGGFVAEASVETGRPTYLLIRPGTSKISTIRSLQEKGATILYGSINDQDIIERVLREHKIEYVISAVGGVSTVDQVNLVKAIKSVGSVKRFLPSEFGHDIDRTEPVEPGLTMYREKRMVRRAVEEAGIPYTYICCNSIAAWPYHDNIHPADVLPPLDRFLIYGDGSVKAYFVAGTDIGRFTMKVLDDDRTLNKTVHFQPPSNLLNVNELASLWEDKIGRKLPRFTVTEADLLTAAKEMPIPKSIVAALTHDIFINGCQINYTMDKEIDVEICSLYPDTPLRTIDQCFDDFLAKVVANRKPIKNPAAKN\n",
      ">TVCU-2006094-Ochna_mossambicensis\n",
      "PNVSEAGRTLVVGSTGFIGRFVTDACIDSGIPTYVLVRSGPVSASKDSAIKSFQDKGATIIHGSMDDQELMMKILKDYKIEFVISAVGGDSILDQIKLVNAVKEVGTVKRFLPSEFGHDIDRADPVEPGLTMYNEKRKVRRCIEEAGIPYTYICCNSIAAWPYHDNTHPADVPPSDQFHIYGDGSVKAYFVAGTDIGKFTVKALTDDRTLNKIVHFQPPSNMLNINELASLWEQKLGLKLPRVTITEDDLLAAAAEMKIPESIVASFTHDIFIKGCQLFSVANPNDVEVCSLYPDTPFRTMDECFTDFAAKIVAENPKIVEK\n",
      ">PEZP-2001925-Glycyrrhiza_glabra\n",
      "NVKRFFPSEFGLDVDRHDAVEPVREVFEDKARIRRVIEAEGVPYTYLCCHAFTGYFLRNLAQLDATVPPRDKVVYVTEADVGTYTIRAANDPRTLNKAVHIRLPANYLTANEVVSLWEKKIGK\n"
     ]
    }
   ],
   "source": [
    "head data/1kP_LAR_orthogroup_random700.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">Ocimum tenuiflorum EGS (AOC97444)\n",
      "MEENGMKSKILIFGGTGYIGNHMVKGSLKLGHPTYVFTRPNSTKTALLNEFQSMGAIIVKGEMGEHEKLV\n",
      "RLMKKVDVVISAVAIPQVLEQLKIVEAIKVAGNIKRFLPSDFGVEEDRISALPPFEACIEKKRIIRRAIE\n",
      "EANIPYTYVSANCFASYFINYLLHPSDPKDEITVYGTGEVKFAMNYEQDIGLYTIKVATDPRTLNRVVIY\n",
      "RPSTNITTQLELISKWEKKIGRNFKKIHIPEEEIVALTKELPEPENIPVAILHCVFIKGVTMSYDFKEND\n",
      "VEASTLYPELKFTTIDELLDIFLHDPPPPASAAF\n",
      "\n",
      ">Clarkia breweri IGS (ABR24112)\n",
      "MEKIIIYGGTGYIGKFMVRASLSFSHPTFIYARPLTPDSTPSSVQLREEFRSMGVTIIEGEMEEHEKMVS\n",
      "VLRQVDVVISALSVPMYPSQLLIIDAIKAAGNIKRFLPSEFGSEEDRIKPLPPFESVLEKKRIIRRAIEA\n"
     ]
    }
   ],
   "source": [
    "head data/Erbils_guide_v3.fasta "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">EGS Ocimum gratissimum (Koeduka et al. 2018)\n",
      ">EGS Ocimum kilimandscharicum (Koeduka et al. 2018)\n",
      ">EGS Ocimum basillicum (Koeduka et al. 2018)\n",
      ">EGS Gymnadenia conopsea (Koeduka et al. 2018)\n",
      ">EGS Gymnadenia odoratissima (Koeduka et al. 2018)\n",
      ">PLR Linum perenne (Koeduka et al. 2018)\n",
      ">P(L)R Arabidopsis thaliana (Koeduka et al. 2018)\n",
      ">LAR Vitis vinifera (Koeduka et al. 2018)\n",
      ">LAR Desmodium uncinatum\n",
      ">LAR Malus domestica (Koeduka et al. 2018)\n",
      ">LAR Fragaria x ananassa (Koeduka et al. 2018)\n",
      ">IFR-LIKE Oryza sativa (Kim et al. 2009)\n",
      ">IFR-LIKE Ginkgo biloba (Hua et al. 2013)\n",
      ">LAR Pinus taeda\n",
      ">LAR Picea abies\n",
      ">LAR Azolla caroliniana (CVEG-2016335)\n"
     ]
    }
   ],
   "source": [
    "grep '>' data/Erbils_verified_sequences_v2_nonredundantwithv3.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">Ocimum tenuiflorum EGS (AOC97444)\n",
      ">Clarkia breweri IGS (ABR24112)\n",
      ">Clarkia breweri EGS1 (ABR24113)\n",
      ">Piper regnellii PPS (AHA90809)\n",
      ">Larrea tridentata PPS (AHA90806)\n",
      ">Pimpinella anisum AIS (ACL13526)\n",
      ">Petunia hybrida IGS (ABD17322)\n",
      ">Gymnadenia densiflora EGS (AKB11749)\n",
      ">Fragaria ananassa EGS1 (AGV02007)\n",
      ">Larrea tridentata APS (AHA90804)\n",
      ">Lotus japonicus PTR (BAF34844)\n",
      ">Clarkia breweri EGS2 (ABR24114)\n",
      ">Populus trichorpa PCBER (CAA06707)\n",
      ">Forsythia intermedia PCBER (AAF64174)\n",
      ">Petunia hybrida EGS (ABR24115)\n",
      ">Nicotiana tabacum PCBER (BAG84267)\n",
      ">Fragaria ananassa EGS2 (AGV02008)\n",
      ">Piper regnellii APS (AHA90807)\n",
      ">Pinus taeda PCBER (AAC32591)\n",
      ">Tsuga heterophylla PCBER (AAF64179)\n",
      ">Cicer arietinum IFR (Q00016)\n",
      ">Pisum sativum IFR (P52576)\n",
      ">Medicago sativa IFR (CAA41106)\n",
      ">Tsuga heterophylla PLR (AAF64184)\n",
      ">Thuja plicata PLR (AAF63509)\n",
      ">Forsythia intermedia PLR (AAC49608)\n",
      ">Linum album PLR (CAH60857)\n",
      ">Desmodium uncinatum LAR (Q84V83)\n",
      ">Fragaria ananassa LAR (ABH07785)\n",
      ">Selaginella moellendorffii PLR-like1 (438560)\n",
      ">Selaginella moellendorffii PLR-like2 (230858)\n",
      ">Chlamydomonas reinhardtii PLR-like (379876)\n",
      ">Chlorella variabilis PLR-like (42998)\n",
      ">Coccomyxa subellipsoidea PLR-like (23909)\n",
      ">Klebsormidium nitens PLR-like (kfl00017_0610)\n",
      ">Azolla filiculoides LAR\n"
     ]
    }
   ],
   "source": [
    "grep '>' data/Erbils_guide_v3.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat data/Erbils_guide_v3.fasta > data/1kP_LAR_orthogroup_random700_guidev4.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat data/Erbils_verified_sequences_v2_nonredundantwithv3.fasta >> data/1kP_LAR_orthogroup_random700_guidev4.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "echo \"\\n\" >> data/1kP_LAR_orthogroup_random700_guidev4.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat data/1kP_LAR_orthogroup_random700.fasta >> data/1kP_LAR_orthogroup_random700_guidev4.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "718\n"
     ]
    }
   ],
   "source": [
    "grep '>' -c ./data/1kP_LAR_orthogroup_random700_guidev4.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1kP_LAR_orthogroup_random700_guidev4\n"
     ]
    }
   ],
   "source": [
    "inseq=\"$inseq\"_random\"$n\"_guidev4\n",
    "echo $inseq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">Ocimum tenuiflorum EGS (AOC97444)\n",
      "MEENGMKSKILIFGGTGYIGNHMVKGSLKLGHPTYVFTRPNSTKTALLNEFQSMGAIIVKGEMGEHEKLV\n",
      "RLMKKVDVVISAVAIPQVLEQLKIVEAIKVAGNIKRFLPSDFGVEEDRISALPPFEACIEKKRIIRRAIE\n",
      "EANIPYTYVSANCFASYFINYLLHPSDPKDEITVYGTGEVKFAMNYEQDIGLYTIKVATDPRTLNRVVIY\n",
      "RPSTNITTQLELISKWEKKIGRNFKKIHIPEEEIVALTKELPEPENIPVAILHCVFIKGVTMSYDFKEND\n",
      "VEASTLYPELKFTTIDELLDIFLHDPPPPASAAF\n",
      "\n",
      ">Clarkia breweri IGS (ABR24112)\n",
      "MEKIIIYGGTGYIGKFMVRASLSFSHPTFIYARPLTPDSTPSSVQLREEFRSMGVTIIEGEMEEHEKMVS\n",
      "VLRQVDVVISALSVPMYPSQLLIIDAIKAAGNIKRFLPSEFGSEEDRIKPLPPFESVLEKKRIIRRAIEA\n"
     ]
    }
   ],
   "source": [
    "head data/$inseq.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Aligning\n",
    "\n",
    "### Either run everything one by one..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First with MAFFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "------------------------------------------------------------------------------\n",
      "  MAFFT v7.429 (2019/Jul/1)\n",
      "  https://mafft.cbrc.jp/alignment/software/\n",
      "  MBE 30:772-780 (2013), NAR 30:3059-3066 (2002)\n",
      "------------------------------------------------------------------------------\n",
      "High speed:\n",
      "  % mafft in > out\n",
      "  % mafft --retree 1 in > out (fast)\n",
      "\n",
      "High accuracy (for <~200 sequences x <~2,000 aa/nt):\n",
      "  % mafft --maxiterate 1000 --localpair  in > out (% linsi in > out is also ok)\n",
      "  % mafft --maxiterate 1000 --genafpair  in > out (% einsi in > out)\n",
      "  % mafft --maxiterate 1000 --globalpair in > out (% ginsi in > out)\n",
      "\n",
      "If unsure which option to use:\n",
      "  % mafft --auto in > out\n",
      "\n",
      "--op # :         Gap opening penalty, default: 1.53\n",
      "--ep # :         Offset (works like gap extension penalty), default: 0.0\n",
      "--maxiterate # : Maximum number of iterative refinement, default: 0\n",
      "--clustalout :   Output: clustal format, default: fasta\n",
      "--reorder :      Outorder: aligned, default: input order\n",
      "--quiet :        Do not report progress\n",
      "--thread # :     Number of threads (if unsure, --thread -1)\n",
      "--dash :         Add structural information (Rozewicki et al, submitted)\n"
     ]
    },
    {
     "ename": "",
     "evalue": "1",
     "output_type": "error",
     "traceback": []
    }
   ],
   "source": [
    "mafft --help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is probably the most acurate mafft setting, which is turned off by default in normal or auto mafft for alignments bigger than 200 sequences. Since We're moving to a publication quality tree, let's do is this way and see. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(phylogenetics) "
     ]
    },
    {
     "ename": "",
     "evalue": "1",
     "output_type": "error",
     "traceback": []
    }
   ],
   "source": [
    "conda activate phylogenetics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(phylogenetics) (phylogenetics) outputhat23=16\n",
      "treein = 0\n",
      "compacttree = 0\n",
      "stacksize: 8192 kb\n",
      "rescale = 1\n",
      "All-to-all alignment.\n",
      "tbfast-pair (aa) Version 7.455\n",
      "alg=L, model=BLOSUM62, 2.00, -0.10, +0.10, noshift, amax=0.0\n",
      "6 thread(s)\n",
      "\n",
      "outputhat23=16\n",
      "Loading 'hat3.seed' ... \n",
      "done.\n",
      "Writing hat3 for iterative refinement\n",
      "rescale = 1\n",
      "Gap Penalty = -1.53, +0.00, +0.00\n",
      "tbutree = 1, compacttree = 0\n",
      "Constructing a UPGMA tree ... \n",
      "  130 / 137\n",
      "done.\n",
      "\n",
      "Progressive alignment ... \n",
      "STEP   135 /136 (thread    3) \n",
      "Reallocating (by thread 2) ..done. *alloclen = 1927\n",
      "STEP   136 /136 (thread    2) \n",
      "done.\n",
      "tbfast (aa) Version 7.455\n",
      "alg=A, model=BLOSUM62, 1.53, -0.00, -0.00, noshift, amax=0.0\n",
      "6 thread(s)\n",
      "\n",
      "minimumweight = 0.000010\n",
      "autosubalignment = 0.000000\n",
      "nthread = 6\n",
      "randomseed = 0\n",
      "blosum 62 / kimura 200\n",
      "poffset = 0\n",
      "niter = 16\n",
      "sueff_global = 0.100000\n",
      "nadd = 16\n",
      "Loading 'hat3' ... done.\n",
      "rescale = 1\n",
      "\n",
      "  130 / 137\n",
      "Segment   1/  1    1- 628\n",
      "004-0270-1 (thread    1) better        \n",
      "Converged2.\n",
      "done\n",
      "dvtditr (aa) Version 7.455\n",
      "alg=A, model=BLOSUM62, 1.53, -0.00, -0.00, noshift, amax=0.0\n",
      "6 thread(s)\n",
      "\n",
      "\n",
      "Strategy:\n",
      " L-INS-i (Probably most accurate, very slow)\n",
      " Iterative refinement method (<16) with LOCAL pairwise alignment information\n",
      "\n",
      "If unsure which option to use, try 'mafft --auto input > output'.\n",
      "For more information, see 'mafft --help', 'mafft --man' and the mafft page.\n",
      "\n",
      "The default gap scoring scheme has been changed in version 7.110 (2013 Oct).\n",
      "It tends to insert more gaps into gap-rich regions than previous versions.\n",
      "To disable this change, add the --leavegappyregion option.\n",
      "\n",
      "(phylogenetics) "
     ]
    },
    {
     "ename": "",
     "evalue": "1",
     "output_type": "error",
     "traceback": []
    }
   ],
   "source": [
    "#rm \"./data/alignments_raw/$inseq\"_aligned-mafft-linsi.fasta\n",
    "if    [ ! -d ./data/alignments_raw/ ]\n",
    "then  mkdir  ./data/alignments_raw\n",
    "fi\n",
    "if    [ ! -f \"./data/alignments_raw/$inseq\"_aligned-mafft.fasta ]\n",
    "then  linsi --thread 6 data/$inseq.fasta > \"./data/alignments_raw/$inseq\"_aligned-mafft-linsi.fasta\n",
    "fi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1kP_LAR_orthogroup_manual-selection-1_guidev4_aligned-mafft-linsi.fasta\n",
      "1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi.fasta\n",
      "1kP_LAR_orthogroup_random700_guidev4_aligned-mafft-linsi.fasta\n",
      "(phylogenetics) "
     ]
    },
    {
     "ename": "",
     "evalue": "1",
     "output_type": "error",
     "traceback": []
    }
   ],
   "source": [
    "ls ./data/alignments_raw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">Ocimum tenuiflorum EGS (AOC97444)\n",
      "M-----------------------------------------------------------\n",
      "-----------------------EENGMKSKILIFGGTGYIGNHMVKGSLK------LGH\n",
      "-----------------------------PTYVFTRP----------------NSTKTAL\n",
      "LNEFQSMGAIIVKGEMG-E----------------HEKLVRLMK----------------\n",
      "---------KVDVVISAVAIPQ----------------VLEQLKIVEAIKVA--------\n",
      "---G--NIKRFLPSDF-------GVEEDR--IS-ALP-PFEA----CIEKKRIIRRAIEE\n",
      "ANIPYTYVSANCFASYFINYLLHPS---------DPK----D-EITVYGTGEVKFAMNYE\n",
      "QDIGLYTIKVATDPRTLNRVVIYRPSTNIT---TQLELISKWEKK------IGRNFKKIH\n",
      "IPEEEIVALT--KELPEP-E-NIPVAILHCVF-IKGVTMSYDFKEN-D-VEASTLYPELK\n",
      "(phylogenetics) "
     ]
    },
    {
     "ename": "",
     "evalue": "1",
     "output_type": "error",
     "traceback": []
    }
   ],
   "source": [
    "head ./data/alignments_raw/\"$inseq\"_aligned-mafft-linsi.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alignment trimming\n",
    "\n",
    "odds are, your alignment is quite gappy which may confuse tree building algorithms. Often it is better to remove gappy columns in your alignment. Let's have a look at this with `trimAl`. Short for 'trim alignment' No Artificial intelegence stuff going on here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "trimAl v1.4.rev15 build[2013-12-17]. 2009-2013. Salvador Capella-Gutierrez and Toni Gabaldón.\n",
      "\n",
      "trimAl webpage: http://trimal.cgenomics.org\n",
      "\n",
      "This program is free software: you can redistribute it and/or modify \n",
      "it under the terms of the GNU General Public License as published by \n",
      "the Free Software Foundation, the last available version.\n",
      "\n",
      "Please cite:\n",
      "\t\ttrimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses.\n",
      "\t\tSalvador Capella-Gutierrez; Jose M. Silla-Martinez; Toni Gabaldon.\n",
      "\t\tBioinformatics 2009, 25:1972-1973.\n",
      "\n",
      "Basic usage\n",
      "\ttrimal -in <inputfile> -out <outputfile> -(other options).\n",
      "\n",
      "Common options (for a complete list please see the User Guide or visit http://trimal.cgenomics.org):\n",
      "\n",
      "    -h                       Print this information and show some examples.\n",
      "    --version                Print the trimAl version.\n",
      "\n",
      "    -in <inputfile>          Input file in several formats (clustal, fasta, NBRF/PIR, nexus, phylip3.2, phylip).\n",
      "\n",
      "    -compareset <inputfile>  Input list of paths for the files containing the alignments to compare.\n",
      "    -forceselect <inputfile> Force selection of the given input file in the files comparison method.\n",
      "\n",
      "    -backtrans <inputfile>   Use a Coding Sequences file to get a backtranslation for a given AA alignment\n",
      "    -ignorestopcodon         Ignore stop codons in the input coding sequences\n",
      "    -splitbystopcodon        Split input coding sequences up to first stop codon appearance\n",
      "\n",
      "    -matrix <inpufile>       Input file for user-defined similarity matrix (default is Blosum62).\n",
      "\n",
      "    -out <outputfile>        Output alignment in the same input format (default stdout). (default input format)\n",
      "    -htmlout <outputfile>    Get a summary of trimal's work in an HTML file.\n",
      "\n",
      "    -keepheader              Keep original sequence header including non-alphanumeric characters.\n",
      "                             Only available for input FASTA format files. (future versions will extend this feature)\n",
      "\n",
      "    -nbrf                    Output file in NBRF/PIR format\n",
      "    -mega                    Output file in MEGA format\n",
      "    -nexus                   Output file in NEXUS format\n",
      "    -clustal                 Output file in CLUSTAL format\n",
      "\n",
      "    -fasta                   Output file in FASTA format\n",
      "    -fasta_m10               Output file in FASTA format. Sequences name length up to 10 characters.\n",
      "\n",
      "    -phylip                  Output file in PHYLIP/PHYLIP4 format\n",
      "    -phylip_m10              Output file in PHYLIP/PHYLIP4 format. Sequences name length up to 10 characters.\n",
      "    -phylip_paml             Output file in PHYLIP format compatible with PAML\n",
      "    -phylip_paml_m10         Output file in PHYLIP format compatible with PAML. Sequences name length up to 10 characters.\n",
      "    -phylip3.2               Output file in PHYLIP3.2 format\n",
      "    -phylip3.2_m10           Output file in PHYLIP3.2 format. Sequences name length up to 10 characters.\n",
      "\n",
      "    -complementary           Get the complementary alignment.\n",
      "    -colnumbering            Get the relationship between the columns in the old and new alignment.\n",
      "\n",
      "    -selectcols { n,l,m-k }  Selection of columns to be removed from the alignment. Range: [0 - (Number of Columns - 1)]. (see User Guide).\n",
      "    -selectseqs { n,l,m-k }  Selection of sequences to be removed from the alignment. Range: [0 - (Number of Sequences - 1)]. (see User Guide).\n",
      "\n",
      "    -gt -gapthreshold <n>    1 - (fraction of sequences with a gap allowed). Range: [0 - 1]\n",
      "    -st -simthreshold <n>    Minimum average similarity allowed. Range: [0 - 1]\n",
      "    -ct -conthreshold <n>    Minimum consistency value allowed.Range: [0 - 1]\n",
      "    -cons <n>                Minimum percentage of the positions in the original alignment to conserve. Range: [0 - 100]\n",
      "\n",
      "    -nogaps                  Remove all positions with gaps in the alignment.\n",
      "    -noallgaps               Remove columns composed only by gaps.\n",
      "    -keepseqs                Keep sequences even if they are composed only by gaps.\n",
      "\n",
      "    -gappyout                Use automated selection on \"gappyout\" mode. This method only uses information based on gaps' distribution. (see User Guide).\n",
      "    -strict                  Use automated selection on \"strict\" mode. (see User Guide).\n",
      "    -strictplus              Use automated selection on \"strictplus\" mode. (see User Guide).\n",
      "                             (Optimized for Neighbour Joining phylogenetic tree reconstruction).\n",
      "\n",
      "    -automated1              Use a heuristic selection of the automatic method based on similarity statistics. (see User Guide). (Optimized for Maximum Likelihood phylogenetic tree reconstruction).\n",
      "\n",
      "    -terminalonly            Only columns out of internal boundaries (first and last column without gaps) are \n",
      "                             candidated to be trimmed depending on the applied method\n",
      "    -block <n>               Minimum column block size to be kept in the trimmed alignment. Available with manual and automatic (gappyout) methods\n",
      "\n",
      "    -resoverlap              Minimum overlap of a positions with other positions in the column to be considered a \"good position\". Range: [0 - 1]. (see User Guide).\n",
      "    -seqoverlap              Minimum percentage of \"good positions\" that a sequence must have in order to be conserved. Range: [0 - 100](see User Guide).\n",
      "\n",
      "    -clusters <n>            Get the most Nth representatives sequences from a given alignment. Range: [1 - (Number of sequences)]\n",
      "    -maxidentity <n>         Get the representatives sequences for a given identity threshold. Range: [0 - 1].\n",
      "\n",
      "    -w <n>                   (half) Window size, score of position i is the average of the window (i - n) to (i + n).\n",
      "    -gw <n>                  (half) Window size only applies to statistics/methods based on Gaps.\n",
      "    -sw <n>                  (half) Window size only applies to statistics/methods based on Similarity.\n",
      "    -cw <n>                  (half) Window size only applies to statistics/methods based on Consistency.\n",
      "\n",
      "    -sgc                     Print gap scores for each column in the input alignment.\n",
      "    -sgt                     Print accumulated gap scores for the input alignment.\n",
      "    -ssc                     Print similarity scores for each column in the input alignment.\n",
      "    -sst                     Print accumulated similarity scores for the input alignment.\n",
      "    -sfc                     Print sum-of-pairs scores for each column from the selected alignment\n",
      "    -sft                     Print accumulated sum-of-pairs scores for the selected alignment\n",
      "    -sident                  Print identity scores for all sequences in the input alignment. (see User Guide).\n",
      "\n",
      "Some Examples:\n",
      "\n",
      "1) Removes all positions in the alignment with gaps in 10% or more of\n",
      "   the sequences, unless this leaves less than 60% of original alignment. \n",
      "   In such case, print the 60% best (with less gaps) positions.\n",
      "\n",
      "   trimal -in <inputfile> -out <outputfile> -gt 0.9 -cons 60\n",
      "\n",
      "2) As above but, the gap score is averaged over a window starting\n",
      "   3 positions before and ending 3 positions after each column.\n",
      "\n",
      "   trimal -in <inputfile> -out <outputfile> -gt 0.9 -cons 60 -w 3\n",
      "\n",
      "3) Use an automatic method to decide optimal thresholds, based in the gap scores\n",
      "   from input alignment. (see User Guide for details).\n",
      "\n",
      "   trimal -in <inputfile> -out <outputfile> -gappyout\n",
      "\n",
      "4) Use automatic methods to decide optimal thresholds, based on the combination \n",
      "   of gap and similarity scores. (see User Guide for details).\n",
      "\n",
      "   trimal -in <inputfile> -out <outputfile> -strictplus\n",
      "\n",
      "5) Use an heuristic to decide the optimal method for trimming the alignment. \n",
      "   (see User Guide for details).\n",
      "\n",
      "   trimal -in <inputfile> -out <outputfile> -automated1\n",
      "\n",
      "6) Use residues and sequences overlap thresholds to delete some sequences from the \n",
      "   alignemnt. (see User Guide for details).\n",
      "\n",
      "   trimal -in <inputfile> -out <outputfile> -resoverlap 0.8 -seqoverlap 75\n",
      "\n",
      "7) Selection of columns to be deleted from the alignment. The selection can \n",
      "   be a column number or a column number interval. Start from 0\n",
      "\n",
      "   trimal -in <inputfile> -out <outputfile> -selectcols { 0,2,3,10,45-60,68,70-78 }\n",
      "\n",
      "8) Get the complementary alignment from the alignment previously trimmed.\n",
      "\n",
      "   trimal -in <inputfile> -out <outputfile> -selectcols { 0,2,3,10,45-60,68,70-78 } -complementary\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "9) Selection of sequences to be deleted from the alignment. Start in 0\n",
      "\n",
      "   trimal -in <inputfile> -out <outputfile> -selectseqs { 2,4,8-12 } \n",
      "\n",
      "10) Select the 5 most representative sequences from the alignment\n",
      "\n",
      "   trimal -in <inputfile> -out <outputfile> -clusters 5 \n",
      "\n",
      "(phylogenetics) "
     ]
    },
    {
     "ename": "",
     "evalue": "1",
     "output_type": "error",
     "traceback": []
    }
   ],
   "source": [
    "trimal -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'm trying some manual trimming as alternative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(phylogenetics) (phylogenetics) trimming alignment data/alignments_raw/1kP_LAR_orthogroup_manual-selection-1_guidev4_aligned-mafft-linsi.fasta\n",
      "\n",
      "WARNING: Removing column '50' composed only by gaps\n",
      "WARNING: Removing column '51' composed only by gaps\n",
      "WARNING: Removing column '52' composed only by gaps\n",
      "WARNING: Removing column '53' composed only by gaps\n",
      "WARNING: Removing column '54' composed only by gaps\n",
      "WARNING: Removing column '55' composed only by gaps\n",
      "WARNING: Removing column '56' composed only by gaps\n",
      "WARNING: Removing column '57' composed only by gaps\n",
      "WARNING: Removing column '58' composed only by gaps\n",
      "WARNING: Removing column '59' composed only by gaps\n",
      "WARNING: Removing column '60' composed only by gaps\n",
      "WARNING: Removing column '61' composed only by gaps\n",
      "WARNING: Removing column '62' composed only by gaps\n",
      "WARNING: Removing column '63' composed only by gaps\n",
      "WARNING: Removing column '64' composed only by gaps\n",
      "WARNING: Removing column '65' composed only by gaps\n",
      "WARNING: Removing column '66' composed only by gaps\n",
      "WARNING: Removing column '67' composed only by gaps\n",
      "WARNING: Removing column '68' composed only by gaps\n",
      "WARNING: Removing column '69' composed only by gaps\n",
      "WARNING: Removing column '70' composed only by gaps\n",
      "WARNING: Removing column '71' composed only by gaps\n",
      "WARNING: Removing column '72' composed only by gaps\n",
      "WARNING: Removing column '73' composed only by gaps\n",
      "WARNING: Removing column '74' composed only by gaps\n",
      "WARNING: Removing column '75' composed only by gaps\n",
      "WARNING: Removing column '76' composed only by gaps\n",
      "WARNING: Removing column '77' composed only by gaps\n",
      "WARNING: Removing column '78' composed only by gaps\n",
      "WARNING: Removing column '79' composed only by gaps\n",
      "WARNING: Removing column '111' composed only by gaps\n",
      "WARNING: Removing column '112' composed only by gaps\n",
      "WARNING: Removing column '113' composed only by gaps\n",
      "WARNING: Removing column '114' composed only by gaps\n",
      "WARNING: Removing column '115' composed only by gaps\n",
      "WARNING: Removing column '116' composed only by gaps\n",
      "WARNING: Removing column '120' composed only by gaps\n",
      "WARNING: Removing column '121' composed only by gaps\n",
      "WARNING: Removing column '122' composed only by gaps\n",
      "WARNING: Removing column '123' composed only by gaps\n",
      "WARNING: Removing column '124' composed only by gaps\n",
      "WARNING: Removing column '125' composed only by gaps\n",
      "WARNING: Removing column '126' composed only by gaps\n",
      "WARNING: Removing column '127' composed only by gaps\n",
      "WARNING: Removing column '128' composed only by gaps\n",
      "WARNING: Removing column '129' composed only by gaps\n",
      "WARNING: Removing column '130' composed only by gaps\n",
      "WARNING: Removing column '131' composed only by gaps\n",
      "WARNING: Removing column '132' composed only by gaps\n",
      "WARNING: Removing column '133' composed only by gaps\n",
      "WARNING: Removing column '134' composed only by gaps\n",
      "WARNING: Removing column '135' composed only by gaps\n",
      "WARNING: Removing column '136' composed only by gaps\n",
      "WARNING: Removing column '137' composed only by gaps\n",
      "WARNING: Removing column '138' composed only by gaps\n",
      "WARNING: Removing column '139' composed only by gaps\n",
      "WARNING: Removing column '140' composed only by gaps\n",
      "WARNING: Removing column '141' composed only by gaps\n",
      "WARNING: Removing column '142' composed only by gaps\n",
      "WARNING: Removing column '143' composed only by gaps\n",
      "WARNING: Removing column '144' composed only by gaps\n",
      "WARNING: Removing column '145' composed only by gaps\n",
      "WARNING: Removing column '146' composed only by gaps\n",
      "WARNING: Removing column '147' composed only by gaps\n",
      "WARNING: Removing column '199' composed only by gaps\n",
      "WARNING: Removing column '200' composed only by gaps\n",
      "WARNING: Removing column '201' composed only by gaps\n",
      "WARNING: Removing column '202' composed only by gaps\n",
      "WARNING: Removing column '203' composed only by gaps\n",
      "WARNING: Removing column '204' composed only by gaps\n",
      "WARNING: Removing column '205' composed only by gaps\n",
      "WARNING: Removing column '206' composed only by gaps\n",
      "WARNING: Removing column '207' composed only by gaps\n",
      "WARNING: Removing column '208' composed only by gaps\n",
      "WARNING: Removing column '209' composed only by gaps\n",
      "WARNING: Removing column '210' composed only by gaps\n",
      "WARNING: Removing column '211' composed only by gaps\n",
      "WARNING: Removing column '212' composed only by gaps\n",
      "WARNING: Removing column '213' composed only by gaps\n",
      "WARNING: Removing column '214' composed only by gaps\n",
      "WARNING: Removing column '225' composed only by gaps\n",
      "WARNING: Removing column '226' composed only by gaps\n",
      "WARNING: Removing column '227' composed only by gaps\n",
      "WARNING: Removing column '228' composed only by gaps\n",
      "WARNING: Removing column '229' composed only by gaps\n",
      "WARNING: Removing column '230' composed only by gaps\n",
      "WARNING: Removing column '231' composed only by gaps\n",
      "WARNING: Removing column '232' composed only by gaps\n",
      "WARNING: Removing column '233' composed only by gaps\n",
      "WARNING: Removing column '234' composed only by gaps\n",
      "WARNING: Removing column '235' composed only by gaps\n",
      "WARNING: Removing column '236' composed only by gaps\n",
      "WARNING: Removing column '237' composed only by gaps\n",
      "WARNING: Removing column '238' composed only by gaps\n",
      "WARNING: Removing column '239' composed only by gaps\n",
      "WARNING: Removing column '240' composed only by gaps\n",
      "WARNING: Removing column '241' composed only by gaps\n",
      "WARNING: Removing column '242' composed only by gaps\n",
      "WARNING: Removing column '243' composed only by gaps\n",
      "WARNING: Removing column '244' composed only by gaps\n",
      "WARNING: Removing column '245' composed only by gaps\n",
      "WARNING: Removing column '246' composed only by gaps\n",
      "WARNING: Removing column '247' composed only by gaps\n",
      "WARNING: Removing column '266' composed only by gaps\n",
      "WARNING: Removing column '267' composed only by gaps\n",
      "WARNING: Removing column '268' composed only by gaps\n",
      "WARNING: Removing column '269' composed only by gaps\n",
      "WARNING: Removing column '270' composed only by gaps\n",
      "WARNING: Removing column '271' composed only by gaps\n",
      "WARNING: Removing column '272' composed only by gaps\n",
      "WARNING: Removing column '273' composed only by gaps\n",
      "WARNING: Removing column '274' composed only by gaps\n",
      "WARNING: Removing column '275' composed only by gaps\n",
      "WARNING: Removing column '276' composed only by gaps\n",
      "WARNING: Removing column '297' composed only by gaps\n",
      "WARNING: Removing column '298' composed only by gaps\n",
      "WARNING: Removing column '299' composed only by gaps\n",
      "WARNING: Removing column '300' composed only by gaps\n",
      "WARNING: Removing column '301' composed only by gaps\n",
      "WARNING: Removing column '316' composed only by gaps\n",
      "WARNING: Removing column '317' composed only by gaps\n",
      "WARNING: Removing column '318' composed only by gaps\n",
      "WARNING: Removing column '319' composed only by gaps\n",
      "WARNING: Removing column '320' composed only by gaps\n",
      "WARNING: Removing column '321' composed only by gaps\n",
      "WARNING: Removing column '322' composed only by gaps\n",
      "WARNING: Removing column '450' composed only by gaps\n",
      "WARNING: Removing column '451' composed only by gaps\n",
      "WARNING: Removing column '452' composed only by gaps\n",
      "\n",
      "(phylogenetics) "
     ]
    },
    {
     "ename": "",
     "evalue": "1",
     "output_type": "error",
     "traceback": []
    }
   ],
   "source": [
    "mkdir data/alignments_trimmed 2> /dev/null\n",
    "\n",
    "for a in \"data/alignments_raw/$inseq\"_aligned*.fasta\n",
    "do  appendix=$(echo $a | cut -d '/' -f 3- | sed \"s/$inseq\\_//\" | sed \"s/.fasta//\")\n",
    "    if    [ ! -f data/alignments_trimmed/\"$inseq\"_\"$appendix\"_trim-gt6-seq80.fasta ]\n",
    "    then  echo \"trimming alignment $a\"\n",
    "          sed -i 's/ /_/g' $a\n",
    "          trimal -in $a   \\\n",
    "                 -out data/alignments_trimmed/\"$inseq\"_\"$appendix\"_trim-gt6-seq80.fasta \\\n",
    "                 -gt .6 \\\n",
    "                 -seqoverlap 80 \\\n",
    "                 -resoverlap .6 \\\n",
    "                 -htmlout data/alignments_trimmed/\"$inseq\"_\"$appendix\"_trim-gt6-seq80.html \n",
    "    fi\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6.html\n",
      "1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80.fasta\n",
      "1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80.html\n",
      "1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq90.html\n",
      "1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt7.html\n",
      "1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt8.html\n",
      "1kP_LAR_orthogroup_random700_guidev4_aligned-mafft-linsi_trim-gt6-seq80.fasta\n",
      "1kP_LAR_orthogroup_random700_guidev4_aligned-mafft-linsi_trim-gt6-seq80.html\n",
      "(phylogenetics) "
     ]
    },
    {
     "ename": "",
     "evalue": "1",
     "output_type": "error",
     "traceback": []
    }
   ],
   "source": [
    "ls data/alignments_trimmed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">Ocimum_tenuiflorum_EGS_(AOC97444)\n",
      "KILIFGGTGYIGNHMVKGSLKLGHPTYVFTRP----TKTALLNEFQSMGAIIVKGEMGEH\n",
      "EKLVRLMKKVDVVISAVAIPQVLEQLKIVEAIKVAGNIKRFLPSDFGVEEDRISALPPFE\n",
      "ACIEKKRIIRRAIEEANIPYTYVSANCFASYFINYLLHPSD---PKDEITVYGTGEVKFA\n",
      "MNYEQDIGLYTIKVATDPRTLNRVVIYRPSTNITTQLELISKWEKKIGRNFKKIHIPEEE\n",
      "IVALTKELPEPENIPVAILHCVFIKGVTMSYDFKENDVEASTLYPELKFTTIDELLDIF\n",
      ">Clarkia_breweri_IGS_(ABR24112)\n",
      "KIIIYGGTGYIGKFMVRASLSFSHPTFIYARPLTPSSSVQLREEFRSMGVTIIEGEMEEH\n",
      "EKMVSVLRQVDVVISALSVPMYPSQLLIIDAIKAAGNIKRFLPSEFGSEEDRIKPLPPFE\n",
      "SVLEKKRIIRRAIEAAELPYTYVSANCFGAYFVNYLLHPSN---RDDDIVIYGTGETKFV\n",
      "(phylogenetics) "
     ]
    },
    {
     "ename": "",
     "evalue": "1",
     "output_type": "error",
     "traceback": []
    }
   ],
   "source": [
    "head data/alignments_trimmed/\"$inseq\"_aligned*.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [],
   "source": [
    "conda deactivate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tree building\n",
    "We'll make fast trees (not so acurate, no bootstraps, but fast) (optional) and we'll make \"propper trees\" using the amazing `iqtree`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## fasttree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm -rf ./data/LAR_orthogroup_fasttrees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls analyses/1kP_LAR_orthogroup_random700_guidev3_fasttrees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for a in data/alignments_trimmed/\"$inseq\"_aligned*.fasta\n",
    "do echo \"making a fasttree of file $a\"\n",
    "   appendix=$(echo $a | cut -d '/' -f 3- | sed \"s/$inseq\\_//\" | sed \"s/.fasta//\")\n",
    "   echo $appendix\n",
    "   if   [ ! -d analyses/\"$inseq\"_fasttrees ]\n",
    "   then mkdir  analyses/\"$inseq\"_fasttrees\n",
    "   fi\n",
    "   if   [ ! -d analyses/\"$inseq\"_fasttrees/\"$appendix\" ]\n",
    "   then mkdir  analyses/\"$inseq\"_fasttrees/\"$appendix\" \n",
    "        fasttree -log analyses/\"$inseq\"_fasttrees/\"$appendix\"/\"$inseq\"_\"$appendix\"_fasttree.log \\\n",
    "                $a \\\n",
    "                > analyses/\"$inseq\"_fasttrees/\"$appendix\"/\"$inseq\"_\"$appendix\"_fasttree.tree \\\n",
    "                2> analyses/\"$inseq\"_fasttrees/\"$appendix\"/\"$inseq\"_\"$appendix\"_fasttree.stderr &\n",
    "    fi\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tail analyses/\"$inseq\"_fasttrees/\"$appendix\"/\"$inseq\"_\"$appendix\"_fasttree.log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ls analyses/\"$inseq\"_fasttrees/\"$appendix\"/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IQtree \n",
    "\n",
    "Now let's loop over all made aligments, making a tree. Choose your parameters wisely, this can take a long time\n",
    "\n",
    " * Are doing a first testrun? Have a look at the -fast option\n",
    " * No time to wait for propper bootstraps, use the ultrafast bootstrap option, -bb 1000 (Don't believe nodes with less than 95% bootstrap\n",
    " * Writing for your amazing manuscript? Do propper bootstraps! -b 1000\n",
    " \n",
    "Model finder, this is the best feature of iqtree. \n",
    "Evolution can happen in a lot of ways, and iqtree takes this into account. \n",
    "Use `-m TEST` to use modelfinder, or `-m MFP` for extended modelfinding (for your publication quality tree). \n",
    "Modelfinder, especially the extended version, can take a long time to calculate. \n",
    "So if you have done it once for a specific alignment, don't do it twice. ;)\n",
    "\n",
    "Now have a look at these examples and the manual:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### iqtree ultra fast bootstrap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data/alignments_trimmed/1kP_LAR_orthogroup_random700_guidev4_aligned-mafft-linsi_trim-gt6-seq80.fasta\n"
     ]
    }
   ],
   "source": [
    "ls data/alignments_trimmed/\"$inseq\"_aligned*gt*.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best-fit model: LG+R9 chosen according to BIC\n"
     ]
    }
   ],
   "source": [
    "grep 'Best-fit model'  ./analyses/1kP_LAR_orthogroup_random700_guidev3_trees/aligned-mafft-linsi_trim-gt6-seq80/1kP_LAR_orthogroup_random700_guidev3_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "making a tree of file data/alignments_trimmed/1kP_LAR_orthogroup_manual-selection-1_guidev4_aligned-mafft-linsi_trim-gt6-seq80.fasta\n",
      "aligned-mafft-linsi_trim-gt6-seq80\n",
      "[1] 13959\n"
     ]
    }
   ],
   "source": [
    "for a in data/alignments_trimmed/\"$inseq\"_aligned*gt*.fasta\n",
    "do echo \"making a tree of file $a\"\n",
    "   appendix=$(echo $a | cut -d '/' -f 3- | sed \"s/$inseq\\_//\" | sed \"s/.fasta//\")\n",
    "   echo $appendix\n",
    "   if   [ ! -d analyses/\"$inseq\"_trees ]\n",
    "   then mkdir  analyses/\"$inseq\"_trees\n",
    "   fi\n",
    "   if   [ ! -d analyses/\"$inseq\"_trees/\"$appendix\" ]\n",
    "   then mkdir  analyses/\"$inseq\"_trees/\"$appendix\" \n",
    "   fi\n",
    "   if   [ ! -f analyses/\"$inseq\"_trees/\"$appendix\"/\"$inseq\"_\"$appendix\"_iqtree-bb1000.tree ]\n",
    "   then     nice iqtree -s $a \\\n",
    "            -b 1000 \\\n",
    "            -nt AUTO \\\n",
    "            -ntmax 12  \\\n",
    "            -pre analyses/\"$inseq\"_trees/\"$appendix\"/\"$inseq\"_\"$appendix\"_iqtree-b1000 \\\n",
    "            2>  analyses/\"$inseq\"_trees/\"$appendix\"/\"$inseq\"_\"$appendix\"_iqtree-b1000 \\\n",
    "            >  analyses/\"$inseq\"_trees/\"$appendix\"/\"$inseq\"_\"$appendix\"_iqtree-b1000 \\\n",
    "            -m MFP && \\\n",
    "            cat analyses/\"$inseq\"_trees/\"$appendix\"/\"$inseq\"_\"$appendix\"_iqtree-b1000.log | mail -s LARguide4 laura.w.dijkhuizen@gmail.com &\n",
    "   fi\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1kP_LAR_orthogroup_manual-selection-1_guidev4_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000\n",
      "1kP_LAR_orthogroup_manual-selection-1_guidev4_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000.log\n",
      "1kP_LAR_orthogroup_manual-selection-1_guidev4_aligned-mafft-linsi_trim-gt6-seq80_iqtree-bb2000-shalrt2000.model.gz\n"
     ]
    }
   ],
   "source": [
    "ls analyses/\"$inseq\"_trees/aligned*/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Log-likelihood of BIONJ tree: -26054.931\n",
      "--------------------------------------------------------------------\n",
      "|             INITIALIZING CANDIDATE TREE SET                      |\n",
      "--------------------------------------------------------------------\n",
      "Generating 98 parsimony trees... 3.243 second\n",
      "Computing log-likelihood of 98 initial trees ... 3.834 seconds\n",
      "Current best score: -26018.011\n",
      "\n",
      "Do NNI search on 20 best initial trees\n",
      "Estimate model parameters (epsilon = 0.100)\n",
      "BETTER TREE FOUND at iteration 1: -25937.370\n",
      "Estimate model parameters (epsilon = 0.100)\n",
      "BETTER TREE FOUND at iteration 2: -25929.496\n",
      "Estimate model parameters (epsilon = 0.100)\n",
      "BETTER TREE FOUND at iteration 7: -25914.412\n",
      "Iteration 10 / LogL: -25979.020 / Time: 0h:0m:11s\n",
      "Iteration 20 / LogL: -25955.453 / Time: 0h:0m:14s\n",
      "Finish initializing candidate tree set (20)\n",
      "Current best tree score: -25914.412 / CPU time: 13.561\n",
      "Number of iterations: 20\n",
      "--------------------------------------------------------------------\n",
      "|               OPTIMIZING CANDIDATE TREE SET                      |\n",
      "--------------------------------------------------------------------\n",
      "Estimate model parameters (epsilon = 0.100)\n",
      "BETTER TREE FOUND at iteration 26: -25909.379\n",
      "Iteration 30 / LogL: -25937.171 / Time: 0h:0m:19s (0h:1m:3s left)\n",
      "Estimate model parameters (epsilon = 0.100)\n",
      "BETTER TREE FOUND at iteration 35: -25907.896\n",
      "Iteration 40 / LogL: -25916.054 / Time: 0h:0m:23s (0h:0m:57s left)\n",
      "Iteration 50 / LogL: -25909.251 / Time: 0h:0m:28s (0h:0m:49s left)\n",
      "Estimate model parameters (epsilon = 0.100)\n",
      "BETTER TREE FOUND at iteration 53: -25907.783\n",
      "Estimate model parameters (epsilon = 0.100)\n",
      "BETTER TREE FOUND at iteration 59: -25907.629\n",
      "Iteration 60 / LogL: -25932.891 / Time: 0h:0m:33s (0h:0m:55s left)\n",
      "Iteration 70 / LogL: -25926.043 / Time: 0h:0m:37s (0h:0m:48s left)\n",
      "UPDATE BEST LOG-LIKELIHOOD: -25907.623\n",
      "Iteration 80 / LogL: -25942.090 / Time: 0h:0m:42s (0h:0m:42s left)\n",
      "Iteration 90 / LogL: -25907.844 / Time: 0h:0m:46s (0h:0m:35s left)\n",
      "Iteration 100 / LogL: -25907.908 / Time: 0h:0m:50s (0h:0m:30s left)\n"
     ]
    }
   ],
   "source": [
    "tail -n 40 analyses/\"$inseq\"_trees/aligned-mafft-*_trim-gt*/\"$inseq\"_*_iqtree-b10*log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "207.710\n",
      "68.437\n",
      "160.606\n",
      "159.044\n",
      "59.495\n",
      "129.266\n",
      "104.288\n",
      "89.461\n",
      "75.402\n",
      "156.090\n",
      "127.875\n",
      "73.773\n",
      "164.789\n",
      "94.132\n",
      "89.039\n",
      "100.568\n",
      "91.968\n",
      "98.221\n",
      "222.949\n",
      "70.330\n",
      "102.509\n",
      "67.632\n",
      "204.925\n",
      "114.973\n",
      "162.648\n",
      "72.103\n",
      "128.652\n",
      "106.543\n",
      "70.537\n",
      "96.770\n",
      "73.556\n",
      "89.199\n",
      "106.254\n",
      "103.882\n",
      "170.749\n",
      "101.082\n",
      "130.859\n",
      "170.035\n",
      "91.127\n",
      "130.286\n",
      "160.742\n",
      "155.138\n",
      "154.546\n",
      "132.545\n",
      "75.157\n",
      "169.577\n",
      "102.092\n",
      "94.485\n",
      "120.175\n",
      "99.854\n",
      "91.182\n",
      "113.628\n",
      "74.572\n",
      "111.872\n",
      "79.561\n",
      "128.137\n",
      "165.133\n",
      "75.980\n",
      "94.105\n",
      "104.848\n",
      "168.801\n",
      "137.418\n",
      "117.718\n",
      "116.673\n",
      "100.831\n",
      "110.485\n",
      "73.302\n",
      "195.893\n",
      "94.550\n",
      "168.349\n",
      "93.181\n",
      "147.269\n",
      "78.141\n",
      "125.551\n",
      "101.201\n",
      "137.044\n",
      "255.232\n",
      "70.162\n",
      "157.198\n",
      "86.679\n",
      "67.720\n",
      "156.570\n",
      "127.188\n",
      "107.918\n",
      "97.744\n",
      "87.710\n",
      "114.682\n",
      "135.268\n",
      "101.482\n",
      "133.839\n",
      "121.002\n",
      "142.569\n",
      "86.912\n",
      "93.676\n",
      "92.130\n",
      "60.067\n",
      "70.612\n",
      "73.490\n",
      "102.705\n",
      "167.626\n",
      "104.834\n",
      "76.348\n",
      "94.770\n",
      "76.464\n",
      "106.696\n",
      "84.185\n",
      "92.886\n",
      "160.983\n",
      "108.668\n",
      "74.362\n",
      "156.502\n",
      "71.418\n",
      "129.707\n",
      "87.569\n",
      "107.085\n"
     ]
    }
   ],
   "source": [
    "grep START analyses/\"$inseq\"_trees/aligned-mafft-*_trim-gt*/\"$inseq\"_*_iqtree-b10*log -B 2 | grep 'Total' | cut -d ':' -f 2 | cut -d ' ' -f 2 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# tree storage\n",
    "\n",
    "find the trees in iTOL: https://itol.embl.de/shared/lauradijkhuizen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Refinement\n",
    "\n",
    "Based on the analysis above, I'd like to make a tree of a subset of this tree. I'll be using: \n",
    " * Erbils verified sequences\n",
    " * Representatives of the major clades in the LAR and PCBER section of the tree.\n",
    " * For those representatives I'll select several genera, and then take all sequences of that genus for my subset.\n",
    " \n",
    " These genera are\n",
    "  * Eucalyptus\n",
    "  * Pinus\n",
    "  * Linum\n",
    "  * Blechnum\n",
    "  * Davallia \n",
    "  * Azolla\n",
    "  * Dipteris\n",
    "  * Ceratopteris\n",
    "  * Osmunda\n",
    "  * Quercus\n",
    "  * Populus\n",
    "  * Freycinetia \n",
    "  * Lendenbergia\n",
    "  * Salix\n",
    "  * Galax\n",
    "  \n",
    "  Nevermind, that's way too much. I'm selecting on these spefific key words\n",
    "  \n",
    "  * Azolla\n",
    "  * Blechnum\n",
    "  * Ceratopteris\n",
    "  * Linum_tenuifolium\n",
    "  * Eucalyptus_leucoxylon\n",
    "  * Pinus_radiata\n",
    "  * Salix_acutifolia\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      4 Azolla_cf._caroliniana\n",
      "     16 Blechnum_spicant\n",
      "      9 BVOF\n",
      "     10 Ceratopteris_thalictroides\n",
      "     11 Davallia_fejeensis\n",
      "     18 Dipteris_conjugata\n",
      "     56 Eucalyptus_grandis\n",
      "     12 Eucalyptus_leucoxylon\n",
      "      9 Freycinetia_multiflora\n",
      "     17 Galax_urceolata\n",
      "      8 Linum_bienne\n",
      "      5 Linum_flavum\n",
      "      8 Linum_grandiflorum\n",
      "      7 Linum_hirsutum\n",
      "      6 Linum_leonii\n",
      "     10 Linum_lewisii\n",
      "      9 Linum_macraei\n",
      "      8 Linum_strictum\n",
      "      8 Linum_tenuifolium\n",
      "     13 Osmunda_javanica\n",
      "     12 Osmunda_sp.\n",
      "     17 Osmundastrum_cinnamomeum\n",
      "     21 Pinus_jeffreyi\n",
      "     13 Pinus_parviflora\n",
      "     22 Pinus_ponderosa\n",
      "     22 Pinus_radiata\n",
      "     62 Pinus_taeda\n",
      "     22 Populus_trichocarpa\n",
      "     11 POZS\n",
      "     13 Quercus_shumardii\n",
      "     13 Salix_acutifolia\n",
      "      9 Salix_dasyclados\n",
      "     10 Salix_eriocephala\n",
      "      9 Salix_fargesii\n",
      "      7 Salix_purpurea\n",
      "      7 Salix_sachalinensis\n",
      "     10 Salix_viminalis\n"
     ]
    }
   ],
   "source": [
    "grep '>' data/$inseq.fasta | tr -d '>' | grep -f <(echo -e \"Eucalyptus\\nPinus\\nLinum\\nBlechnum\\nDavallia\\nAzolla\\nDipteris\\nCeratopteris\\nOsmunda\\nQuercus\\nPopulus\\nFreycinetia\\nLendenbergia\\nSalix\\nGalax\" ) | cut -f 3 -d '-'| sort | uniq -c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "85\n"
     ]
    }
   ],
   "source": [
    "grep '>' data/$inseq.fasta | tr -d '>' | grep -f <(echo -e \"Azolla\\nBlechnum\\nCeratopteris\\nLinum_tenuifolium\\nEucalyptus_leucoxylon\\nPinus_radiata\\nSalix_acutifolia\" ) | sort | uniq | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So that's 85 sequences from 7 genera + all ERbils Guide sequences, this will be fun. First extract this into a new txt and then a new fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "grep '>' data/$inseq.fasta | tr -d '>' | grep -f <(echo -e \"Azolla\\nBlechnum\\nCeratopteris\\nLinum_tenuifolium\\nEucalyptus_leucoxylon\\nPinus_radiata\\nSalix_acutifolia\" ) | sort | uniq > data/1kP_LAR_orthogroup_manual-selection-1.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [],
   "source": [
    "/opt/rnaseq/bin/extract_sequence_by_id.pl --ID=data/1kP_LAR_orthogroup_manual-selection-1.txt \\\n",
    "                                          --fileformat=txt \\\n",
    "                                          --sequences=<(cat data/$inseq.fasta )\\\n",
    "                                          --outfile=data/\"$inseq\"_manual-selection-1.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "85\n"
     ]
    }
   ],
   "source": [
    "grep '>'  data/\"$inseq\"_manual-selection-1.fasta -c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat data/Erbils_guide_v3.fasta > data/\"$inseq\"_manual-selection-1_guidev4.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "echo \"\\n\" >> data/\"$inseq\"_manual-selection-1_guidev4.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat data/Erbils_verified_sequences_v2_nonredundantwithv3.fasta >> data/\"$inseq\"_manual-selection-1_guidev4.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat data/\"$inseq\"_manual-selection-1.fasta >> data/\"$inseq\"_manual-selection-1_guidev4.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "137\n"
     ]
    }
   ],
   "source": [
    "grep '>' -c data/\"$inseq\"_manual-selection-1_guidev4.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1kP_LAR_orthogroup_manual-selection-1_guidev4\n"
     ]
    }
   ],
   "source": [
    "inseq=\"$inseq\"_manual-selection-1_guidev4\n",
    "echo $inseq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I'm redoing the whole thing again from up-top.\n",
    "\n",
    "stats after linsi mapping and trimal trimming\n",
    "\n",
    "Selected Sequences:   100 /Selected Residues:     304\n",
    "Deleted Sequences:     37 /Deleted Residues:      324\n",
    "\n",
    "I'm quite happy with that, so let's start making a tree."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# other stuff\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Manually set parameters for the script"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copy pasting and adapting the command above to do \n",
    " * no modelfitting, I did that already\n",
    " * propper boostraps, just 200. (Minimum 100)\n",
    " * 12 threads\n",
    " * **new names**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iqtree -s ./alignments_trimmed/\"$inseq\"_\"$appendix\".fasta \\\n",
    "       -b 200 \\\n",
    "       -nt AUTO \\\n",
    "       -ntmax 12 \\\n",
    "       -pre \"$inseq\"_trees/\"$appendix\"/\"$inseq\"_\"$appendix\"_iqtree-b200 \\\n",
    "       -m LG+R9 \\\n",
    "       >      \"$inseq\"_trees/\"$appendix\"/\"$inseq\"_\"$appendix\"_iqtree-b200.stdout \\\n",
    "       && cat \"$inseq\"_trees/\"$appendix\"/\"$inseq\"_\"$appendix\"_iqtree-b200.stdout \\\n",
    "       | mail -s 'LAR tree lauraZolla Desktop' laura &\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "cat -n 80  \"$inseq\"_trees/\"$appendix\"/\"$inseq\"_\"$appendix\"_iqtree-b200.stdout"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### likelihood bootstrap thingy\n",
    "\n",
    "http://www.iqtree.org/doc/Tutorial\n",
    "\n",
    "https://www.ncbi.nlm.nih.gov/pubmed/16785212\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nice iqtree -s ./alignments_trimmed/\"$inseq\"_\"$appendix\".fasta \\\n",
    "       -bb   1000 \\\n",
    "       -alrt 1000 \\\n",
    "       -nt AUTO \\\n",
    "       -ntmax $(nproc) \\\n",
    "       -pre \"$inseq\"_trees/\"$appendix\"/\"$inseq\"_\"$appendix\"_iqtree-bb1000-alrt1000 \\\n",
    "       -m LG+R9 \\\n",
    "       >      \"$inseq\"_trees/\"$appendix\"/\"$inseq\"_\"$appendix\"_iqtree-bb1000-alrt1000.stdout \\\n",
    "       && cat \"$inseq\"_trees/\"$appendix\"/\"$inseq\"_\"$appendix\"_iqtree-bb1000-alrt1000.stdout \\\n",
    "       | mail -s 'LAR tree alrt1000 lauraZolla Desktop' laura &\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tail -n 40 \"$inseq\"_trees/\"$appendix\"/\"$inseq\"_\"$appendix\"_iqtree-bb1000-alrt1000.stdout"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}