Basic_HybSeq_pipeline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A very basic and quick pipeline to derive conservatively called consensus fasta seqeunce from Hyb-Seq data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Requirements: \n",
    "\n",
    "    Trimmomatic-0.36  \n",
    "    bwa mem \n",
    "    samtools  \n",
    "    bcftools  \n",
    "    mafft \n",
    "    trimal\n",
    "    amas  \n",
    "    FastTree  \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Raw reads are trimmed using standard trimmomatic settings, aligned to the baits using BWA mem and a vcf file derived from the bam.\n",
    "The vcf is filtered for quality of at least 30 and to obly include homozygous alternate allele SNPS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Raw reads should in the folder:\n",
    "\n",
    "    ~/Process/Hairdrier/\n",
    "    \n",
    "The expected file name format is:\n",
    "\n",
    "    SAMPLE1_1.fastq.gz, SAMPLE1_2.fastq.gz\n",
    "    \n",
    "The BWA index of the baits should be prepared with the following command:\n",
    "    \n",
    "    bwa index Baits.fna\n",
    "\n",
    "Consensus fastas for each bait for each sample can be produced by making a list of the samples to use and running the following while loop:  \n",
    "    \n",
    "    while read f ; do ./raw_to_fasta_BWA.sh \"$f\" ; done < sample_list\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This produces a multi fasta file for each sample with the form:\n",
    "    >Bait_name1\n",
    "    TCTCACCTACGAATTGATAATGCTGTTCGATGCAATAGCGTTCTCCACTTTTTCAGCATA\n",
    "    ACTCAACGGTCCACAATGTTTATAAAAATAAAAAATAAAAATAACAAACAGTGACTTGAC\n",
    "    >Bait_name2\n",
    "    TAGGAACTTGAGTCGAATCAAGAAAAATCTCGCTTTATCATACGAAGAAGCACCTATCCA\n",
    "    TACCAAAAATAATAAAGTTGGGAAAAAGTTCATTACAAATGGTCCGGGAAACATATCTCG\n",
    "    \n",
    "    \n",
    "This needs to be converted to a set of multifastas for each bait using the script Switch_multi_fasta.py.  Make a folder called By_locus and two test files - one with the list of samples, one with the list of bait names.\n",
    "\n",
    "    mkdir By_locus/\n",
    "    ls *_consensus.fna > fasta_files\n",
    "    grep \">\" SAMPLE_consensus.fna | sed 's/>//g' > locus_list\n",
    "    \n",
    "    python3 Switch_multi_fasta.py\n",
    "\n",
    "This generates many empty files in your working folder which need to be removed with:\n",
    "\n",
    "    find . -size 0 -delete\n",
    "    \n",
    "The by-locus fasta files are in By_locus/ with the form:\n",
    "\n",
    "BAIT_NAME.fasta\n",
    "\n",
    "\\>Sample1.fast_BAIT_NAME\n",
    "TTATCCTCCTTCCCACAGCACACTCATTGACTTCTGTCTCGTGATGCCGTGCTACTCGTTGAAACCTGCTACTGTAGGTACATGAACATGTCCTGAAGCGGCGACATCTCATTTCCACTCTGAAAATGGCTTTCCATCGTCGCATCCATGTAGTAGTTATAAGGAAATTGCAGCGATTTCTCCAATTCCTTCCATTTGGGCTCGCTCTGCACTT\n",
    "\\>Sample2.fast_BAIT_NAME\n",
    "TTATCCTCCTTCCCACAGCACACTCATTGACTTCTGTCTCGTGATGCCGTGCTACTCGTTGAAACCTGCTACTGTAGGTACATGAACATGTCCTGAAGCGGCGACATCTCATTTCCACTCTGAAAATGGCTTTCCATCGTCGCATCCATGTAGTAGTTATAAGGAAATTGCAGCGATTTCTCCAATTCCTTCCATTTGGGCT\n",
    "\n",
    "The fasta headers of these files have the locus name in.  If you want to concatenate these you will want to remove this, leaving just the accession name:\n",
    "\n",
    "    while read f ; do sed ’s/‘\\.fast_\"$f\"’//g’ \"$f\".fasta > \"$f\".fna ; done < ../locus_list\n",
    "    \n",
    "\\>Sample1\n",
    "TTATCCTCCTTCCCACAGCACACTCATTGACTTCTGTCTCGTGATGCCGTGCTACTCGTTGAAACCTGCTACTGTAGGTACATGAACATGTCCTGAAGCGGCGACATCTCATTTCCACTCTGAAAATGGCTTTCCATCGTCGCATCCATGTAGTAGTTATAAGGAAATTGCAGCGATTTCTCCAATTCCTTCCA\n",
    "\\>Sample2\n",
    "TTATCCTCCTTCCCACAGCACACTCATTGACTTCTGTCTCGTGATGCCGGGCTACTCGTAGAAACCTGCTACTGTAGGTACATGAACATGTCCTGGAGCGGCGACATCTCATTTCCACTTTGAAAATGGCTTTCCATCGTCGTATCCATGTAGTAGTTATAAGGAAATTGCAGTGATTTCTCCCATTCCTTCCATTTGGGCTCGCTCTGCACTTC\n",
    "\n",
    "The multifastas can be individually aligned using mafft:\n",
    "\n",
    "    while read f ; do mafft --thread 4 \"$f\".fna > \"$f\"_mafft.fna ; done < locus_list\n",
    "  \n",
    "Strict-trimming of the alignments:  \n",
    "        \n",
    "    while read f ; do trimal -in \"$f\"_mafft.fna -out \"$f\"_mafft_trim.fna -fasta -strict ; done < locus_list \n",
    "        \n",
    "Summarising each alignment:\n",
    "\n",
    "    while read f ; do python3 ~/amas-0.93/amas/AMAS.py -f fasta -d dna -i \"$f\"_mafft_trim.fna -s -o \"$f\".summary.txt ; done < ../locus_list\n",
    "\n",
    "The alignments can be concatenated and a quick fast-tree alignment doen to check the assembly before further analysis\n",
    "The lengths of the alignments are put into the header by trimal.  These need removing to allow concatenation:\n",
    "\n",
    "    while read f ; do sed 's/ [0-9]* bp//g' \"$f\"_mafft_trim.fna > \"$f\".fasta ; done < ../locus_list\n",
    "\n",
    "Concatenating using AMAS.  Put all the alignemed an trimmed fastas into a new folder.  Move into that folder and call AMAS:\n",
    "\n",
    "    python3 ~/amas-0.93/amas/AMAS.py -f fasta -d dna -i *.fasta -c\n",
    "\n",
    "\n",
    "Quick FastTree for checking:\n",
    "\n",
    "    fasttree -gtr -nt  < concat.fasta > concat.tre  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}