step4_expanding-phylogeny.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Aim\n",
    "\n",
    "In step four, I'd like to address expanding the phylogeny with more R2R3 MYB subfamiles to address some open questions: \n",
    "* Can I add a subfamily label to the other Azfi sequences I found before with a custom hmm.\n",
    "* Can I infer the phylogeny to more confidently/explicitly draw the Salvinales VII branch as a part of the remaining VII branch by including other subfamiles and an outgroup\n",
    "\n",
    "In a later notebook, I'll address adding extra data like: \n",
    "Next steps: \n",
    "\n",
    "* Can I add intron/exon structures as a second line of evidence supporting the different subfamilies\n",
    "* Finally, I'd like to log2 expressen/fold change to the Azfi sequences; to illustrate how these subfamilies change expression upon transistion to sexual reproduction at least in _Azolla_\n",
    "\n",
    "Since retrieving sequences from the different databases is quite a pain, \n",
    "I'll first retrieve only NCBI listed sequences (which I can retrieve in batch).\n",
    "Subsequently, I may reitterate in this notebook, expanding on that analysis with sequences from other genome databases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Genes in new list files in the `data` directory were retrieved from NCBI. \n",
    "Acc. nr.s from the J&R study were mixed protein and gene id's, so batch retrieval wasn't possible.\n",
    "Fasta files with both acc nr's were created in the same directory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# prepare data\n",
    "\n",
    "Combining and linearising fasta files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in data/*fasta\n",
    "do  inseq=$(echo $i | cut -d '/' -f 2 | cut -d '.' -f 1)\n",
    "    if   [ ! -f data/\"$inseq\"_linear.fasta ] \n",
    "    then cat data/$inseq.fasta \\\n",
    "         | awk '/^>/ {printf(\"%s%s\\n\",(N>0?\"\\n\":\"\"),$0);N++;next;} {printf(\"%s\",$0);} END {printf(\"\\n\");}' \\\n",
    "         > data/\"$inseq\"_linear.fasta\n",
    "    fi\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data/ARP_sequences_linear.fasta\n",
      "data/Azfi-mybs-subfamVI-suspects_linear.fasta\n",
      "data/Azfi-v1-MYB-sequences_linear.fasta\n",
      "data/Azfi-v1-MYB-sequences_linear_linear.fasta\n",
      "data/CDC5-outgroup_sequences_linear.fasta\n",
      "data/combi_sequences_linear.fasta\n",
      "data/combi_sequences_linear_linear.fasta\n",
      "data/combi-VI-VII-Azfisuspects_linear.fasta\n",
      "data/III_sequences_linear.fasta\n",
      "data/II_sequences_linear.fasta\n",
      "data/I_sequences_linear.fasta\n",
      "data/IV_sequences_linear.fasta\n",
      "data/MYB33_ARATH_linear.fasta\n",
      "data/MYB33_ARATH_linear_linear.fasta\n",
      "data/R1R2R3_sequences_linear.fasta\n",
      "data/VIII_sequences_linear.fasta\n",
      "data/VII_sequences_linear.fasta\n",
      "data/VII_sequences_linear_linear.fasta\n",
      "data/VI_sequences_linear.fasta\n",
      "data/VI_sequences_linear_linear.fasta\n"
     ]
    }
   ],
   "source": [
    "ls data/*_linear.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat data/CDC5-outgroup_sequences_linear.fasta \\\n",
    "    data/I_sequences_linear.fasta    \\\n",
    "    data/II_sequences_linear.fasta   \\\n",
    "    data/III_sequences_linear.fasta  \\\n",
    "    data/IV_sequences_linear.fasta   \\\n",
    "    data/V_sequences_linear.fasta    \\\n",
    "    data/VI_sequences_linear.fasta   \\\n",
    "    data/VII_sequences_linear.fasta  \\\n",
    "    data/VIII_sequences_linear.fasta \\\n",
    "    data/ARP_sequences_linear.fasta  \\\n",
    "    data/Azfi-mybs-subfamVI-suspects_linear.fasta \\\n",
    "    data/R1R2R3_sequences_linear.fasta \\\n",
    "    data/MYB33_ARATH_linear.fasta \\\n",
    "    > data/combi-I-to-VIII-Azfi_sequences_linear.fasta\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "inseq='combi-I-to-VIII-Azfi_sequences_linear'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(phylogenetics) (phylogenetics) (phylogenetics) "
     ]
    }
   ],
   "source": [
    "conda activate phylogenetics\n",
    "if    [ ! -d ./data/alignments_raw/ ]\n",
    "then  mkdir  ./data/alignments_raw\n",
    "fi\n",
    "for   i in data/combi*sequences_linear.fasta\n",
    "do    if    [ ! -f \"./data/alignments_raw/$inseq\"_aligned-mafft-linsi.fasta ]\n",
    "      then  echo \"aligning $inseq\"\n",
    "            linsi --thread $(nproc) data/$inseq.fasta > ./data/alignments_raw/\"$inseq\"_aligned-mafft-linsi.fasta \\\n",
    "                                                      2> ./data/alignments_raw/\"$inseq\"_aligned-mafft-linsi.log\n",
    "      fi\n",
    "done\n",
    "conda deactivate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(phylogenetics) (phylogenetics) "
     ]
    }
   ],
   "source": [
    "conda activate phylogenetics\n",
    "einsi --thread $(nproc) data/$inseq.fasta > ./data/alignments_raw/\"$inseq\"_aligned-mafft-einsi.fasta \\\n",
    "                                                      2> ./data/alignments_raw/\"$inseq\"_aligned-mafft-einsi.log\n",
    "conda deactivate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(jalview) (jalview) "
     ]
    }
   ],
   "source": [
    "conda activate jalview\n",
    "for   i in data/alignments_raw/*.fasta\n",
    "do    prefix=$(echo $i | sed 's/\\.fasta//')\n",
    "      if    [ ! -f $prefix.png ]\n",
    "      then  jalview -nodisplay \\\n",
    "                    -open $prefix.fasta \\\n",
    "                    -colour CLUSTAL \\\n",
    "                    -png  $prefix.png > /dev/null 2> /dev/null\n",
    "      fi\n",
    "done\n",
    "conda deactivate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linsi alignment looks like this:\n",
    "![](./data/alignments_raw/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi.png)\n",
    "\n",
    "And the einsi alignment looks like this:\n",
    "![](./data/alignments_raw/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-einsi.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linsi all. is even more gappy than the alignments we've seen before, but that stands to reason since a lot more of sequence diversity is covered in this MSA. \n",
    "After trimming this will look just fine I'm sure.\n",
    "\n",
    "The einsi alignment seems to concentrate more sequence information into columns.\n",
    "From this overview perspective, that looks quite good to me."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(phylogenetics) (phylogenetics) (phylogenetics) (phylogenetics) (phylogenetics) (phylogenetics) trimming alignment data/alignments_raw/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-einsi.fasta\n",
      "[1] 172740\n",
      "trimming alignment data/alignments_raw/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi.fasta\n",
      "[2] 172747\n",
      "(phylogenetics) [1]-  Done                    trimal -in $a -out data/alignments_trimmed/\"$inseq\"_\"$appendix\"_\"$trimappendix\".fasta -gt .4 -htmlout data/alignments_trimmed/\"$inseq\"_\"$appendix\"_\"$trimappendix\".html\n",
      "[2]+  Done                    trimal -in $a -out data/alignments_trimmed/\"$inseq\"_\"$appendix\"_\"$trimappendix\".fasta -gt .4 -htmlout data/alignments_trimmed/\"$inseq\"_\"$appendix\"_\"$trimappendix\".html\n",
      "(phylogenetics) "
     ]
    }
   ],
   "source": [
    "conda activate phylogenetics\n",
    "if    [ ! -d data/alignments_trimmed ]\n",
    "then  mkdir  data/alignments_trimmed \n",
    "fi\n",
    "\n",
    "# define appendix only once here:\n",
    "trimappendix='trim-gt4'\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\"_\"$trimappendix\".fasta ]\n",
    "    then  echo \"trimming alignment $a\"\n",
    "          sed -i 's/ /_/g' $a\n",
    "          trimal -in $a   \\\n",
    "                 -out data/alignments_trimmed/\"$inseq\"_\"$appendix\"_\"$trimappendix\".fasta \\\n",
    "                 -gt .4 \\\n",
    "                 -htmlout data/alignments_trimmed/\"$inseq\"_\"$appendix\"_\"$trimappendix\".html &\n",
    "    fi\n",
    "done\n",
    "wait\n",
    "conda deactivate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(jalview) [1] 172804\n",
      "[2] 172809\n",
      "[3] 172823\n",
      "[4] 172835\n",
      "[5] 172842\n",
      "[6] 172857\n",
      "(jalview) [1]   Done                    jalview -nodisplay -open $prefix.fasta -colour CLUSTAL -png $prefix.png > /dev/null 2> /dev/null\n",
      "[5]-  Done                    jalview -nodisplay -open $prefix.fasta -colour CLUSTAL -png $prefix.png > /dev/null 2> /dev/null\n",
      "[2]   Done                    jalview -nodisplay -open $prefix.fasta -colour CLUSTAL -png $prefix.png > /dev/null 2> /dev/null\n",
      "[3]   Done                    jalview -nodisplay -open $prefix.fasta -colour CLUSTAL -png $prefix.png > /dev/null 2> /dev/null\n",
      "[4]-  Done                    jalview -nodisplay -open $prefix.fasta -colour CLUSTAL -png $prefix.png > /dev/null 2> /dev/null\n",
      "[6]+  Done                    jalview -nodisplay -open $prefix.fasta -colour CLUSTAL -png $prefix.png > /dev/null 2> /dev/null\n",
      "(jalview) "
     ]
    }
   ],
   "source": [
    "conda activate jalview\n",
    "for   i in data/alignments_trimmed/*.fasta\n",
    "do    prefix=$(echo $i | sed 's/\\.fasta//')\n",
    "      if    [ ! -f $prefix.png ]\n",
    "      then  jalview -nodisplay \\\n",
    "                    -open $prefix.fasta \\\n",
    "                    -colour CLUSTAL \\\n",
    "                    -png  $prefix.png > /dev/null 2> /dev/null &\n",
    "      fi\n",
    "done\n",
    "wait\n",
    "conda deactivate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 40%\n",
    "Starting at gap threshol 40%; the linsi alignment looks like this:\n",
    "![](./data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi_trim-gt4.png)\n",
    "and the einsi like so:\n",
    "![](./data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-einsi_trim-gt4.png)\n",
    "\n",
    "### 60%\n",
    "Moving on to 60%; the linsi alignment looks like this:\n",
    "![](./data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi_trim-gt6.png)\n",
    "and the einsi like so:\n",
    "![](./data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-einsi_trim-gt6.png)\n",
    "\n",
    "### 80%\n",
    "Finally at 80%; the linsi alignment looks like this:\n",
    "![](./data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi_trim-gt8.png)\n",
    "and the einsi like so:\n",
    "![](./data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-einsi_trim-gt8.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing these trimmings of both the einsi and the linsi allignment, I'm thinking that a raw or very leniently trimmed allignment retains the differentiating regions as described in J&R.\n",
    "However, I also think that the more strictly trimmed regions may be better for phylogeny purposes.\n",
    "Perhaps that the 80% is a bit too much, I'll loose some differentiating capability. \n",
    "So 60 is the middle ground then?\n",
    "I'll run trees for all since I can miss the time and see what method gets me the best signal.\n",
    "I'm putting my money on 60% einsi.\n",
    "\n",
    "~~But first, I'll remove those two spurious seqeunces. from the alignments and trimms.~~\n",
    "I have removed one spurious sequence and replaced the other with the correct one.\n",
    "It appeared I had extracted a wrong sequence from ncbi.\n",
    "\n",
    "# Tree inference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(phylogenetics) Making a directory aligned-mafft-einsi_trim-gt4 to store trees (name based on alignment filename)\n",
      "making a tree of file data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-einsi_trim-gt4.fasta\n",
      "Making a directory aligned-mafft-einsi_trim-gt6 to store trees (name based on alignment filename)\n",
      "making a tree of file data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-einsi_trim-gt6.fasta\n",
      "Making a directory aligned-mafft-einsi_trim-gt8 to store trees (name based on alignment filename)\n",
      "making a tree of file data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-einsi_trim-gt8.fasta\n",
      "making a tree of file data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi_trim-gt4.fasta\n",
      "Making a directory aligned-mafft-linsi_trim-gt6 to store trees (name based on alignment filename)\n",
      "making a tree of file data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi_trim-gt6.fasta\n",
      "Making a directory aligned-mafft-linsi_trim-gt8 to store trees (name based on alignment filename)\n",
      "making a tree of file data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi_trim-gt8.fasta\n",
      "(phylogenetics) "
     ]
    }
   ],
   "source": [
    "conda activate phylogenetics\n",
    "for a in data/alignments_trimmed/\"$inseq\"_aligned*gt*.fasta\n",
    "do  #iqpendix='iqtree-b100'\n",
    "    iqpendix='iqtree-bb2000-alrt2000'\n",
    "\n",
    "\n",
    "\n",
    "    file_appendix=$(echo $a | cut -d '/' -f 3- | sed \"s/$inseq\\_//\" | sed \"s/.fasta//\")\n",
    "    \n",
    "    if   [ ! -d    analyses/\"$inseq\"_trees/\"$file_appendix\" ]\n",
    "    then echo \"Making a directory $file_appendix to store trees (name based on alignment filename)\"\n",
    "         mkdir -p  analyses/\"$inseq\"_trees/\"$file_appendix\" \n",
    "    fi\n",
    "\n",
    "    iqprefix=analyses/\"$inseq\"_trees/\"$file_appendix\"/\"$inseq\"_\"$file_appendix\"_\"$iqpendix\"\n",
    "    if   [ ! -f \"$iqprefix\".treefile ]\n",
    "    then echo \"making a tree of file $a\"\n",
    "         nice iqtree -s $a \\\n",
    "                     -m MFP \\\n",
    "                     -bb 2000 -alrt 2000 \\\n",
    "                     -nt AUTO \\\n",
    "                     -ntmax $(nproc)  \\\n",
    "                     -pre  \"$iqprefix\" \\\n",
    "                     2>   \"$iqprefix\".stderr \\\n",
    "                      >    \"$iqprefix\".stdout\n",
    "    #cat \"$iqprefix\".log | mail -s \"IQtree_run $a\" laura\n",
    "    fi\n",
    "done\n",
    "conda deactivate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This [initial tree](https://itol.embl.de/tree/1312115964244391596537363#) based on the `combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi_trim-gt4` alignment, puts some sequences at odd places. \n",
    "I'm wondering if I have a too gappy alignment here, like I already suspected for earlier alignments.\n",
    "Alternativelly, I may try to align sequences with the einsi method in stead of the linsi method, since the big conserved domain is quite gappy and linsi may not align conserved sites outside of the main domain propperly.\n",
    "\n",
    "On second thought, I reaslise that some sequences are very oddly aligned.\n",
    "Perhaps I made a typo extracting these, or they could be incomplete."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# conclusions\n",
    "\n",
    "Comparing the different trimming percentages, it seems that 60% minumum sequence content per column gives me the better internal consistency and phylogenetic signal than the 40% alignments for both einsi and linsi.\n",
    "\n",
    "In the 80% tree, some branch lengths understandibly become versy short. \n",
    "Bootstraps are similar to the 60% tree.\n",
    "\n",
    "The ARP clade isn't consistently placed and since it is not of interest for me, I would opt to take it out.\n",
    "Similarly, one CDC5 sequence behaves odly: it is consistently placed with subfamily VIII sequences.\n",
    "I suspect it was miscategorised and I will omit it from this analysis.\n",
    "\n",
    "Comparing the einsi to the linsi alignment, I find that the linsi alignment places the ARP clade often on very early branches, while the einsi alignment puts it often with the VII category of sequences. \n",
    "Topology wise, thei einsi method reproduced the structure of older clades as shown in J&R, but the linsi method has better internal consistency in those regions. \n",
    "Alternativelly, the ARP clade may weaken internal consitency in whatever regions it pops up.\n",
    "In that case I'm overinterpreting a side effect of including this oddly behaving clade of sequences.\n",
    "I'll reinspect the MSAs and choose a method. \n",
    "Additionally, I'll run non-parametric trees on both data sets but remove the ARP sequences before.\n",
    "\n",
    "Biologically speaking, the Azfi sequences of interest end up in the same clades as they did before. \n",
    "The VII salviniales clade is placed nicely with the other VII sequences with OK support.\n",
    "The newly introduced R2R3 myb sequences from the Azfi genome (result from the hmm scan) are all placed in the VIII subfamily; answering that question. \n",
    "This meets the two aims I set for this notebook; but I'll double check still with two nonparametric tree runs.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[35m\u001b[Kaligned-mafft-einsi_trim-gt4/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-einsi_trim-gt4_iqtree-bb2000-alrt2000.log\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[KBest-fit model\u001b[m\u001b[K: VT+F+R5 chosen according to BIC\n",
      "\u001b[35m\u001b[Kaligned-mafft-einsi_trim-gt6/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-einsi_trim-gt6_iqtree-bb2000-alrt2000.log\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[KBest-fit model\u001b[m\u001b[K: VT+R5 chosen according to BIC\n",
      "\u001b[35m\u001b[Kaligned-mafft-einsi_trim-gt8/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-einsi_trim-gt8_iqtree-bb2000-alrt2000.log\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[KBest-fit model\u001b[m\u001b[K: LG+I+G4 chosen according to BIC\n",
      "\u001b[35m\u001b[Kaligned-mafft-linsi_trim-gt4/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi_trim-gt4_iqtree-bb2000-alrt2000.log\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[KBest-fit model\u001b[m\u001b[K: LG+F+R8 chosen according to BIC\n",
      "\u001b[35m\u001b[Kaligned-mafft-linsi_trim-gt6/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi_trim-gt6_iqtree-bb2000-alrt2000.log\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[KBest-fit model\u001b[m\u001b[K: LG+F+I+G4 chosen according to BIC\n",
      "\u001b[35m\u001b[Kaligned-mafft-linsi_trim-gt8/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi_trim-gt8_iqtree-bb2000-alrt2000.log\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[KBest-fit model\u001b[m\u001b[K: LG+G4 chosen according to BIC\n"
     ]
    }
   ],
   "source": [
    "cd ./analyses/combi-I-to-VIII-Azfi_sequences_linear_trees/\n",
    "grep 'Best-fit model' */*.log\n",
    "cd ../../"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(phylogenetics) Making a directory aligned-mafft-einsi_trim-gt6-seqrmmanual to store trees (name based on alignment filename)\n",
      "making a tree of file data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-einsi_trim-gt6-seqrmmanual.fasta\n",
      "Making a directory aligned-mafft-linsi_trim-gt6-seqrmmanual to store trees (name based on alignment filename)\n",
      "making a tree of file data/alignments_trimmed/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi_trim-gt6-seqrmmanual.fasta\n",
      "(phylogenetics) "
     ]
    }
   ],
   "source": [
    "conda activate phylogenetics\n",
    "for a in data/alignments_trimmed/\"$inseq\"_aligned*gt6-seqrmmanual.fasta\n",
    "do  iqpendix='iqtree-b200'\n",
    "    #iqpendix='iqtree-bb2000-alrt2000'\n",
    "\n",
    "    file_appendix=$(echo $a | cut -d '/' -f 3- | sed \"s/$inseq\\_//\" | sed \"s/.fasta//\")\n",
    "    \n",
    "    if   [ ! -d    analyses/\"$inseq\"_trees/\"$file_appendix\" ]\n",
    "    then echo \"Making a directory $file_appendix to store trees (name based on alignment filename)\"\n",
    "         mkdir -p  analyses/\"$inseq\"_trees/\"$file_appendix\" \n",
    "    fi\n",
    "\n",
    "    iqprefix=analyses/\"$inseq\"_trees/\"$file_appendix\"/\"$inseq\"_\"$file_appendix\"_\"$iqpendix\"\n",
    "    if   [ ! -f \"$iqprefix\".treefile ]\n",
    "    then echo \"making a tree of file $a\"\n",
    "         nice iqtree -s $a  \\\n",
    "                     -m MFP \\\n",
    "                     -b 200 \\\n",
    "                     -nt AUTO \\\n",
    "                     -ntmax $(nproc)  \\\n",
    "                     -pre  \"$iqprefix\" \\\n",
    "                     2>   \"$iqprefix\".stderr \\\n",
    "                      >    \"$iqprefix\".stdout\n",
    "    #cat \"$iqprefix\".log | mail -s \"IQtree_run $a\" laura\n",
    "    fi\n",
    "done\n",
    "conda deactivate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# conclusions\n",
    "Inspecting the two nonparametric trees, they are actually very similar and allow for some solid biological interpretation now.\n",
    "\n",
    "I am convinced that ferns do have the VI subfamily R2R3 Myb, in contradiction to J&R. \n",
    "The Azolla sequences of interest are confidently placed in the VII subfamily.\n",
    "It could be fun to see if there is more Azfi R2R3 mybs I could put in this tree; it's not often you can search a fern genome for such orthologous and examine their differential expressen in context of the fern life cycle.\n",
    "\n",
    "### technical stuff.\n",
    "The spurious CDC5 sequence should go out, it is all over the place and I doubt it is a CDC5 at all.\n",
    "\n",
    "I added one arabidopsis sequence, it may be nice to add a couple well known ones to the different clades to illustrate their possible function.\n",
    "\n",
    "Although the linsi alignment gives me trees with better supported nodes on average (regardless of their position in the tree) the einsi alignment gives me better support for the ancestral nodes.\n",
    "Looking at the alignments again (trimmed at 60% collumn content) I think to see that the einsi alignments more consistently aligns the part just after the main conserved domain.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's add some snapshots of the nonparametric trees in here, first the einsi based tree: [(link)](https://itol.embl.de/tree/1312115964248471596621955)\n",
    "![einsi](analyses/combi-I-to-VIII-Azfi_sequences_linear_trees/aligned-mafft-einsi_trim-gt6-seqrmmanual/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-einsi_trim-gt6-seqrmmanual_iqtree-b200.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then the linsi tree [link](https://itol.embl.de/tree/1312115964249831596621961)\n",
    "![einsi](analyses/combi-I-to-VIII-Azfi_sequences_linear_trees/aligned-mafft-linsi_trim-gt6-seqrmmanual/combi-I-to-VIII-Azfi_sequences_linear_aligned-mafft-linsi_trim-gt6-seqrmmanual_iqtree-b200.png)"
   ]
  }
 ],
 "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
}