Concatenate_alignments.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/work/kelsey/Ribosomal Alignments Plastid Project\n"
     ]
    }
   ],
   "source": [
    "from Bio import AlignIO, SeqIO, Align, Alphabet\n",
    "import pandas as pd\n",
    "import os, re, sys\n",
    "from copy import deepcopy\n",
    "\n",
    "%cd /work/kelsey/Ribosomal\\ Alignments\\ Plastid\\ Project/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "present_genes = []\n",
    "for aln in os.listdir('.'):\n",
    "    if not aln.endswith('.aln') or 'Mito_' in aln:\n",
    "        continue\n",
    "    present_genes.append(aln.replace('.aln', ''))\n",
    "\n",
    "for aln in os.listdir('/work/kelsey/Cyano_Clock_alignments/'):\n",
    "    if not aln.startswith('Analysis14') or not aln.endswith('.fas.fas'):\n",
    "        continue\n",
    "    \n",
    "    gene_name = re.search('Analysis14_(.*).fas.fas', aln).group(1)\n",
    "    if gene_name in present_genes:\n",
    "        continue\n",
    "    \n",
    "    os.system('cp %s/%s ./%s.aln' % ('/work/kelsey/Cyano_Clock_alignments',\n",
    "                                 aln,\n",
    "                                 gene_name)\n",
    "             )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "aln_alphabet = Alphabet.Gapped(Alphabet.IUPAC.ambiguous_dna)\n",
    "\n",
    "genomes    = {}\n",
    "for aln in os.listdir('.'):\n",
    "    if not aln.endswith('.aln') or 'Mito_' in aln:\n",
    "        continue\n",
    "    alignment    = AlignIO.read(aln, 'fasta')\n",
    "    genomes[aln] = set()\n",
    "    for entry in alignment:\n",
    "        if '[' in entry.description and ']' in entry.description:\n",
    "            genome = re.search('\\[(.*)\\]$', entry.description, re.M).group(1).replace(' ', '_')\n",
    "            genome += '_mitochondria'\n",
    "        else:\n",
    "            genome = entry.name\n",
    "\n",
    "        if genome in genomes[aln]:\n",
    "            #print ('\\t**Error, duplicated genome in %s: %s' %(aln, genome))\n",
    "            sys.exit('\\t**Error, duplicated genome in %s: %s' %(aln, genome))\n",
    "\n",
    "        genomes[aln].add(genome)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "genome_union  = set.union(*genomes.values())\n",
    "missing_genes = {} # just to keep track of the number of missing marker genes in each genome\n",
    "concatenation = {}\n",
    "for genome in genome_union:\n",
    "    missing_genes[genome]             = 0\n",
    "    concatenation[genome]             = Align.SeqRecord( Align.Seq('', aln_alphabet) )\n",
    "    concatenation[genome].name        = genome\n",
    "    concatenation[genome].id          = genome\n",
    "    concatenation[genome].description = genome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L22\n",
      "L23\n",
      "S19\n",
      "L18\n",
      "L24\n",
      "S9\n",
      "S8\n",
      "L6\n",
      "L4\n",
      "L5\n",
      "L1\n",
      "L2\n",
      "L3\n",
      "S13\n",
      "S12\n",
      "S7\n",
      "L14\n",
      "S5\n",
      "S10\n",
      "S11\n",
      "S4\n",
      "L29\n",
      "L15\n",
      "S15\n",
      "S14\n",
      "L10\n",
      "S3\n",
      "S17\n",
      "S2\n",
      "L13\n"
     ]
    }
   ],
   "source": [
    "#\n",
    "# fill the handles with the marker sequences from each genome\n",
    "total_genes      = 0.0 # keep track of the number of genes added to the concatenation\n",
    "current_position = 1\n",
    "partitions       = open('concatenated_partitions', 'w')\n",
    "for aln in os.listdir('.'):\n",
    "    if not aln.endswith('.aln') or 'Mito_' in aln:\n",
    "        continue\n",
    "    \n",
    "    print(aln.replace('.aln', ''))\n",
    "    tmp_aln      = AlignIO.read(aln, 'fasta' )\n",
    "    aln_length   = tmp_aln.get_alignment_length() # get the expected size of the alignment so you can compare if all have the same size\n",
    "    total_genes += aln_length\n",
    "\n",
    "    for entry in tmp_aln:\n",
    "        # if this alignment has a different size from the rest, something is reaaaaaly wrong!\n",
    "        if len(entry) != aln_length:\n",
    "            sys.exit('\\t**Error, block \"%s\" has a different length than the rest of the MSA: %s' %(entry.name, aln))\n",
    "\n",
    "        if '[' in entry.description and ']' in entry.description:\n",
    "            genome = re.search('\\[(.*)\\]$', entry.description, re.M).group(1).replace(' ', '_')\n",
    "            genome += '_mitochondria'\n",
    "        else:\n",
    "            genome = entry.name\n",
    "        concatenation[genome] += deepcopy(entry.seq)\n",
    "\n",
    "    partitions.write('LG, %s = %i-%i\\n' %(aln.replace('.aln', ''), current_position, current_position+aln_length-1) )\n",
    "    current_position += aln_length\n",
    "\n",
    "    #\n",
    "    # add gaps for those genomes missing this gene (same size as the expected alignment)\n",
    "    for genome in genome_union.difference(genomes[aln]):\n",
    "        concatenation[genome] += Align.Seq( '-' * aln_length, aln_alphabet )\n",
    "        missing_genes[genome] += aln_length\n",
    "partitions.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t\t**Marchantia_paleacea_subsp_diptera: missing 58.10% from concatenated alignment!\n",
      "\t\t**Lilium_tsingtauense: missing 58.10% from concatenated alignment!\n",
      "\t\t**Mesostigma_viride: missing 51.70% from concatenated alignment!\n",
      "\t\t**Danio_rerio_mitochondria: missing 68.35% from concatenated alignment!\n",
      "\t\t**Cattleya_liliputana: missing 60.65% from concatenated alignment!\n",
      "\t\t**Volvox_carteri_f._nagariensis_mitochondria: missing 68.35% from concatenated alignment!\n",
      "\t\t**Chlamydomonas_reinhardtii_mitochondria: missing 68.35% from concatenated alignment!\n",
      "\t\t**Zygnema_circumcarinatum: missing 56.99% from concatenated alignment!\n",
      "\t\t**Saccharomyces_cerevisiae_S288C_mitochondria: missing 68.35% from concatenated alignment!\n",
      "\t\t**Cyanidioschyzon_merolae_strain_10D_mitochondria: missing 68.35% from concatenated alignment!\n",
      "\t\t**Koliella_corcontica: missing 84.42% from concatenated alignment!\n",
      "\t\t**Physcomitrella_patens: missing 58.10% from concatenated alignment!\n",
      "\t\t**Chara_vulgaris: missing 55.61% from concatenated alignment!\n",
      "\t\t**Oryza_rufipogon: missing 58.10% from concatenated alignment!\n",
      "\t\t**Adiantum_capillus-veneris: missing 58.10% from concatenated alignment!\n",
      "\t\t**Volvox_carteri_f_nagariensis: missing 58.36% from concatenated alignment!\n",
      "\t\t**Amborella_trichopoda: missing 58.10% from concatenated alignment!\n",
      "\t\t**Huperzia_lucidula: missing 59.70% from concatenated alignment!\n",
      "\t\t**Acer_davidii: missing 63.18% from concatenated alignment!\n",
      "19 genomes missing more than 10%\n"
     ]
    }
   ],
   "source": [
    "#\n",
    "# remove genomes missing more than 20% of the marker genes\n",
    "counter = 0\n",
    "for genome, num_missing_genes in missing_genes.items():\n",
    "    if num_missing_genes/total_genes > 0.5:\n",
    "        print('\\t\\t**%s: missing %.2f%% from concatenated alignment!' %(genome, (num_missing_genes/total_genes)*100))\n",
    "        counter += 1\n",
    "\n",
    "print('%i genomes missing more than 10%%' %counter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AlignIO.write(Align.MultipleSeqAlignment(concatenation.values() ), 'concatenated_alignment.aln', 'fasta')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}