notebooks/2018-05-22-Gonad-Methylation-Full-Samples.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bismark\n",
    "\n",
    "In this notebook I'll use Bismark on a subset of my data to ensure the pipeline is reproducible. The steps from the [user guide](https://rawgit.com/FelixKrueger/Bismark/master/Docs/Bismark_User_Guide.html#i-bismark-genome-preparation) are as follows:\n",
    "\n",
    "1. Genome Preparation (already done in [this Jupyter notebook](https://github.com/RobertsLab/project-virginica-oa/blob/master/notebooks/2018-04-27-Gonad-Methylation-Bismark.ipynb))\n",
    "2. Alignment\n",
    "3. Deduplication\n",
    "4. Methlyation Extractor\n",
    "5. HTML Processing Report\n",
    "6. Summary Report"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. Set working directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/Users/yaamini/Documents/project-virginica-oa/notebooks'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/yaamini/Documents/project-virginica-oa/analyses\n"
     ]
    }
   ],
   "source": [
    "cd ../analyses/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mkdir 2018-05-22-Bismark-Full-Samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0516_bismark.err                     \u001b[34m2018-05-01-MethylKit\u001b[m\u001b[m/\r\n",
      "0516_bme.err                         \u001b[34m2018-05-04-Bismark-Full-Samples\u001b[m\u001b[m/\r\n",
      "0516_dedup.err                       \u001b[34m2018-05-22-Bismark-Full-Samples\u001b[m\u001b[m/\r\n",
      "\u001b[34m2018-01-23-MBDSeq-Labwork\u001b[m\u001b[m/           \u001b[34m2018-05-22-Bismark-Subset\u001b[m\u001b[m/\r\n",
      "\u001b[34m2018-04-26-Gonad-Methylation-FastQC\u001b[m\u001b[m/ README.md\r\n",
      "\u001b[34m2018-04-27-Bismark\u001b[m\u001b[m/                  \u001b[34mSamples\u001b[m\u001b[m/\r\n"
     ]
    }
   ],
   "source": [
    "ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples\n"
     ]
    }
   ],
   "source": [
    "cd 2018-05-22-Bismark-Full-Samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Genome Preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This step was already completed in [this Jupyter notebook](https://github.com/RobertsLab/project-virginica-oa/blob/master/notebooks/2018-04-27-Gonad-Methylation-Bismark.ipynb). The genome only needs to be prepared once. I will move on to the second step, alignment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Alignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\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, either version 3 of the License, or\n",
      "     (at your option) any later version.\n",
      "\n",
      "     This program is distributed in the hope that it will be useful,\n",
      "     but WITHOUT ANY WARRANTY; without even the implied warranty of\n",
      "     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n",
      "     GNU General Public License for more details.\n",
      "     You should have received a copy of the GNU General Public License\n",
      "     along with this program.  If not, see <http://www.gnu.org/licenses/>.\n",
      "\n",
      "\n",
      "\n",
      "DESCRIPTION\n",
      "\n",
      "\n",
      "The following is a brief description of command line options and arguments to control the Bismark\n",
      "bisulfite mapper and methylation caller. Bismark takes in FastA or FastQ files and aligns the\n",
      "reads to a specified bisulfite genome. Sequence reads are transformed into a bisulfite converted forward strand\n",
      "version (C->T conversion) or into a bisulfite treated reverse strand (G->A conversion of the forward strand).\n",
      "Each of these reads are then aligned to bisulfite treated forward strand index of a reference genome\n",
      "(C->T converted) and a bisulfite treated reverse strand index of the genome (G->A conversion of the\n",
      "forward strand, by doing this alignments will produce the same positions). These 4 instances of Bowtie (1 or 2)\n",
      "are run in parallel. The sequence file(s) are then read in again sequence by sequence to pull out the original\n",
      "sequence from the genome and determine if there were any protected C's present or not.\n",
      "\n",
      "As of version 0.7.0 Bismark will only run 2 alignment threads for OT and OB in parallel, the 4 strand mode can be\n",
      "re-enabled by using --non_directional.\n",
      "\n",
      "The final output of Bismark is in SAM format by default. For Bowtie 1 one can alos choose to report the old\n",
      "'vanilla' output format, which is a single tab delimited file with all sequences that have a unique best\n",
      "alignment to any of the 4 possible strands of a bisulfite PCR product. Both formats are described in more detail below.\n",
      "\n",
      "\n",
      "USAGE: bismark [options] <genome_folder> {-1 <mates1> -2 <mates2> | <singles>}\n",
      "\n",
      "\n",
      "ARGUMENTS:\n",
      "\n",
      "<genome_folder>          The path to the folder containing the unmodified reference genome\n",
      "                         as well as the subfolders created by the Bismark_Genome_Preparation\n",
      "                         script (/Bisulfite_Genome/CT_conversion/ and /Bisulfite_Genome/GA_conversion/).\n",
      "                         Bismark expects one or more fastA files in this folder (file extension: .fa\n",
      "                         or .fasta). The path can be relative or absolute. The path may also be set as\n",
      "                         '--genome_folder /path/to/genome/folder/'.\n",
      "\n",
      "-1 <mates1>              Comma-separated list of files containing the #1 mates (filename usually includes\n",
      "                         \"_1\"), e.g. flyA_1.fq,flyB_1.fq). Sequences specified with this option must\n",
      "                         correspond file-for-file and read-for-read with those specified in <mates2>.\n",
      "                         Reads may be a mix of different lengths. Bismark will produce one mapping result\n",
      "                         and one report file per paired-end input file pair.\n",
      "\n",
      "-2 <mates2>              Comma-separated list of files containing the #2 mates (filename usually includes\n",
      "                         \"_2\"), e.g. flyA_1.fq,flyB_1.fq). Sequences specified with this option must\n",
      "                         correspond file-for-file and read-for-read with those specified in <mates1>.\n",
      "                         Reads may be a mix of different lengths.\n",
      "\n",
      "<singles>                A comma- or space-separated list of files containing the reads to be aligned (e.g.\n",
      "                         lane1.fq,lane2.fq lane3.fq). Reads may be a mix of different lengths. Bismark will\n",
      "                         produce one mapping result and one report file per input file.\n",
      "\n",
      "\n",
      "OPTIONS:\n",
      "\n",
      "\n",
      "Input:\n",
      "\n",
      "--se/--single_end <list> Sets single-end mapping mode explicitly giving a list of file names as <list>.\n",
      "                         The filenames may be provided as a comma [,] or colon [:] separated list.\n",
      "\n",
      "-q/--fastq               The query input files (specified as <mate1>,<mate2> or <singles> are FASTQ\n",
      "                         files (usually having extension .fg or .fastq). This is the default. See also\n",
      "                         --solexa-quals.\n",
      "\n",
      "-f/--fasta               The query input files (specified as <mate1>,<mate2> or <singles> are FASTA\n",
      "                         files (usually having extensions .fa, .mfa, .fna or similar). All quality values\n",
      "                         are assumed to be 40 on the Phred scale. FASTA files are expected to contain both\n",
      "                         the read name and the sequence on a single line (and not spread over several lines).\n",
      "\n",
      "-s/--skip <int>          Skip (i.e. do not align) the first <int> reads or read pairs from the input.\n",
      "\n",
      "-u/--upto <int>          Only aligns the first <int> reads or read pairs from the input. Default: no limit.\n",
      "\n",
      "--phred33-quals          FASTQ qualities are ASCII chars equal to the Phred quality plus 33. Default: on.\n",
      "\n",
      "--phred64-quals          FASTQ qualities are ASCII chars equal to the Phred quality plus 64. Default: off.\n",
      "\n",
      "--solexa-quals           Convert FASTQ qualities from solexa-scaled (which can be negative) to phred-scaled\n",
      "                         (which can't). The formula for conversion is:\n",
      "                         phred-qual = 10 * log(1 + 10 ** (solexa-qual/10.0)) / log(10). Used with -q. This\n",
      "                         is usually the right option for use with (unconverted) reads emitted by the GA\n",
      "                         Pipeline versions prior to 1.3. Works only for Bowtie 1. Default: off.\n",
      "\n",
      "--solexa1.3-quals        Same as --phred64-quals. This is usually the right option for use with (unconverted)\n",
      "                         reads emitted by GA Pipeline version 1.3 or later. Default: off.\n",
      "\n",
      "--path_to_bowtie         The full path </../../> to the Bowtie (1 or 2) installation on your system. If not\n",
      "                         specified it is assumed that Bowtie (1 or 2) is in the PATH.\n",
      "\n",
      "\n",
      "Alignment:\n",
      "\n",
      "-n/--seedmms <int>       The maximum number of mismatches permitted in the \"seed\", i.e. the first L base pairs\n",
      "                         of the read (where L is set with -l/--seedlen). This may be 0, 1, 2 or 3 and the\n",
      "                         default is 1. This option is only available for Bowtie 1 (for Bowtie 2 see -N).\n",
      "\n",
      "-l/--seedlen             The \"seed length\"; i.e., the number of bases of the high quality end of the read to\n",
      "                         which the -n ceiling applies. The default is 28. Bowtie (and thus Bismark) is faster for\n",
      "                         larger values of -l. This option is only available for Bowtie 1 (for Bowtie 2 see -L).\n",
      "\n",
      "-e/--maqerr <int>        Maximum permitted total of quality values at all mismatched read positions throughout\n",
      "                         the entire alignment, not just in the \"seed\". The default is 70. Like Maq, bowtie rounds\n",
      "                         quality values to the nearest 10 and saturates at 30. This value is not relevant for\n",
      "                         Bowtie 2.\n",
      "\n",
      "--chunkmbs <int>         The number of megabytes of memory a given thread is given to store path descriptors in\n",
      "                         --best mode. Best-first search must keep track of many paths at once to ensure it is\n",
      "                         always extending the path with the lowest cumulative cost. Bowtie tries to minimize the\n",
      "                         memory impact of the descriptors, but they can still grow very large in some cases. If\n",
      "                         you receive an error message saying that chunk memory has been exhausted in --best mode,\n",
      "                         try adjusting this parameter up to dedicate more memory to the descriptors. This value\n",
      "                         is not relevant for Bowtie 2. Default: 512.\n",
      "\n",
      "-I/--minins <int>        The minimum insert size for valid paired-end alignments. E.g. if -I 60 is specified and\n",
      "                         a paired-end alignment consists of two 20-bp alignments in the appropriate orientation\n",
      "                         with a 20-bp gap between them, that alignment is considered valid (as long as -X is also\n",
      "                         satisfied). A 19-bp gap would not be valid in that case. Default: 0.\n",
      "\n",
      "-X/--maxins <int>        The maximum insert size for valid paired-end alignments. E.g. if -X 100 is specified and\n",
      "                         a paired-end alignment consists of two 20-bp alignments in the proper orientation with a\n",
      "                         60-bp gap between them, that alignment is considered valid (as long as -I is also satisfied).\n",
      "                         A 61-bp gap would not be valid in that case. Default: 500.\n",
      "\n",
      "--parallel <int>         (May also be --multicore <int>) Sets the number of parallel instances of Bismark to be run concurrently.\n",
      "                         This forks the Bismark alignment step very early on so that each individual Spawn of Bismark processes\n",
      "                         only every n-th sequence (n being set by --parallel). Once all processes have completed,\n",
      "                         the individual BAM files, mapping reports, unmapped or ambiguous FastQ files are merged\n",
      "                         into single files in very much the same way as they would have been generated running Bismark\n",
      "                         conventionally with only a single instance.\n",
      "\n",
      "                         If system resources are plentiful this is a viable option to speed up the alignment process\n",
      "                         (we observed a near linear speed increase for up to --parallel 8 tested). However, please note\n",
      "                         that a typical Bismark run will use several cores already (Bismark itself, 2 or 4 threads of\n",
      "                         Bowtie/Bowtie2, Samtools, gzip etc...) and ~10-16GB of memory depending on the choice of aligner\n",
      "                         and genome. WARNING: Bismark Parallel (BP?) is resource hungry! Each value of --parallel specified\n",
      "                         will effectively lead to a linear increase in compute and memory requirements, so --parallel 4 for\n",
      "                         e.g. the GRCm38 mouse genome will probably use ~20 cores and eat ~40GB or RAM, but at the same time\n",
      "                         reduce the alignment time to ~25-30%. You have been warned.\n",
      "\n",
      "\n",
      "\n",
      "Bowtie 1 Reporting:\n",
      "\n",
      "-k <2>                   Due to the way Bismark works Bowtie will report up to 2 valid alignments. This option\n",
      "                         will be used by default.\n",
      "\n",
      "--best                   Make Bowtie guarantee that reported singleton alignments are \"best\" in terms of stratum\n",
      "                         (i.e. number of mismatches, or mismatches in the seed in the case if -n mode) and in\n",
      "                         terms of the quality; e.g. a 1-mismatch alignment where the mismatch position has Phred\n",
      "                         quality 40 is preferred over a 2-mismatch alignment where the mismatched positions both\n",
      "                         have Phred quality 10. When --best is not specified, Bowtie may report alignments that\n",
      "                         are sub-optimal in terms of stratum and/or quality (though an effort is made to report\n",
      "                         the best alignment). --best mode also removes all strand bias. Note that --best does not\n",
      "                         affect which alignments are considered \"valid\" by Bowtie, only which valid alignments\n",
      "                         are reported by Bowtie. Bowtie is about 1-2.5 times slower when --best is specified.\n",
      "                         Default: on.\n",
      "\n",
      "--no_best                Disables the --best option which is on by default. This can speed up the alignment process,\n",
      "                         e.g. for testing purposes, but for credible results it is not recommended to disable --best.\n",
      "\n",
      "\n",
      "Output:\n",
      "\n",
      "--non_directional        The sequencing library was constructed in a non strand-specific manner, alignments to all four\n",
      "                         bisulfite strands will be reported. Default: OFF.\n",
      "\n",
      "                         (The current Illumina protocol for BS-Seq is directional, in which case the strands complementary\n",
      "                         to the original strands are merely theoretical and should not exist in reality. Specifying directional\n",
      "                         alignments (which is the default) will only run 2 alignment threads to the original top (OT)\n",
      "                         or bottom (OB) strands in parallel and report these alignments. This is the recommended option\n",
      "                         for sprand-specific libraries).\n",
      "\n",
      "--pbat                   This options may be used for PBAT-Seq libraries (Post-Bisulfite Adapter Tagging; Kobayashi et al.,\n",
      "                         PLoS Genetics, 2012). This is essentially the exact opposite of alignments in 'directional' mode,\n",
      "                         as it will only launch two alignment threads to the CTOT and CTOB strands instead of the normal OT\n",
      "                         and OB ones. Use this option only if you are certain that your libraries were constructed following\n",
      "                         a PBAT protocol (if you don't know what PBAT-Seq is you should not specify this option). The option\n",
      "                         --pbat works only for FastQ files (in both Bowtie and Bowtie 2 mode) and using uncompressed\n",
      "                         temporary files only).\n",
      "\n",
      "--sam-no-hd              Suppress SAM header lines (starting with @). This might be useful when very large input files are\n",
      "                         split up into several smaller files to run concurrently and the output files are to be merged.\n",
      "\n",
      "--rg_tag                 Write out a Read Group tag to the resulting SAM/BAM file. This will write the following line to the\n",
      "                         SAM header: @RG PL: ILLUMINA ID:SAMPLE SM:SAMPLE ; to set ID and SM see --rg_id and --rg_sample.\n",
      "                         In addition each read receives an RG:Z:RG-ID tag. Default: OFF.\n",
      "\n",
      "--rg_id <string>         Sets the ID field in the @RG header line. The default is 'SAMPLE'.\n",
      "\n",
      "--rg_sample <string>     Sets the SM field in the @RG header line; can't be set without setting --rg_id as well. The default is\n",
      "                         'SAMPLE'.\n",
      "\n",
      "--quiet                  Print nothing besides alignments.\n",
      "\n",
      "--vanilla                Performs bisulfite mapping with Bowtie 1 and prints the 'old' output (as in Bismark 0.5.X) instead\n",
      "                         of SAM format output.\n",
      "\n",
      "-un/--unmapped           Write all reads that could not be aligned to a file in the output directory. Written reads will\n",
      "                         appear as they did in the input, without any translation of quality values that may have\n",
      "                         taken place within Bowtie or Bismark. Paired-end reads will be written to two parallel files with _1\n",
      "                         and _2 inserted in their filenames, i.e. _unmapped_reads_1.txt and unmapped_reads_2.txt. Reads\n",
      "                         with more than one valid alignment with the same number of lowest mismatches (ambiguous mapping)\n",
      "                         are also written to _unmapped_reads.txt unless the option --ambiguous is specified as well.\n",
      "\n",
      "--ambiguous              Write all reads which produce more than one valid alignment with the same number of lowest\n",
      "                         mismatches or other reads that fail to align uniquely to a file in the output directory.\n",
      "                         Written reads will appear as. they did in the input, without any of the translation of quality\n",
      "                         values that may have taken place within Bowtie or Bismark. Paired-end reads will be written to two\n",
      "                         parallel files with _1 and _2 inserted in their filenames, i.e. _ambiguous_reads_1.txt and\n",
      "                         _ambiguous_reads_2.txt. These reads are not written to the file specified with --un.\n",
      "\n",
      "-o/--output_dir <dir>    Write all output files into this directory. By default the output files will be written into\n",
      "                         the same folder as the input file(s). If the specified folder does not exist, Bismark will attempt\n",
      "                         to create it first. The path to the output folder can be either relative or absolute.\n",
      "\n",
      "--temp_dir <dir>         Write temporary files to this directory instead of into the same directory as the input files. If\n",
      "                         the specified folder does not exist, Bismark will attempt to create it first. The path to the\n",
      "                         temporary folder can be either relative or absolute.\n",
      "\n",
      "--non_bs_mm              Optionally outputs an extra column specifying the number of non-bisulfite mismatches a read during the\n",
      "                         alignment step. This option is only available for SAM format. In Bowtie 2 context, this value is\n",
      "                         just the number of actual non-bisulfite mismatches and ignores potential insertions or deletions.\n",
      "                         The format for single-end reads and read 1 of paired-end reads is 'XA:Z:number of mismatches'\n",
      "                         and 'XB:Z:number of mismatches' for read 2 of paired-end reads.\n",
      "\n",
      "--gzip                   Temporary bisulfite conversion files will be written out in a GZIP compressed form to save disk\n",
      "                         space. This option is available for most alignment modes but is not available for paired-end FastA\n",
      "                         files. This option might be somewhat slower than writing out uncompressed files, but this awaits\n",
      "                         further testing.\n",
      "\n",
      "--sam                    The output will be written out in SAM format instead of the default BAM format. Bismark will\n",
      "                         attempt to use the path to Samtools that was specified with '--samtools_path', or, if it hasn't\n",
      "                         been specified, attempt to find Samtools in the PATH. If no installation of Samtools can be found,\n",
      "                         the SAM output will be compressed with GZIP instead (yielding a .sam.gz output file).\n",
      "\n",
      "--cram                   Writes the output to a CRAM file instead of BAM. This requires the use of Samtools 1.2 or higher.\n",
      "\n",
      "--cram_ref <ref_file>    CRAM output requires you to specify a reference genome as a single FastA file. If this single-FastA\n",
      "                         reference file is not supplied explicitly it will be regenerated from the genome .fa sequence(s)\n",
      "                         used for the Bismark run and written to a file called 'Bismark_genome_CRAM_reference.mfa' into the\n",
      "                         oputput directory.\n",
      "\n",
      "--samtools_path          The path to your Samtools installation, e.g. /home/user/samtools/. Does not need to be specified\n",
      "                         explicitly if Samtools is in the PATH already.\n",
      "\n",
      "--prefix <prefix>        Prefixes <prefix> to the output filenames. Trailing dots will be replaced by a single one. For\n",
      "                         example, '--prefix test' with 'file.fq' would result in the output file 'test.file.fq_bismark.sam' etc.\n",
      "\n",
      "-B/--basename <basename> Write all output to files starting with this base file name. For example, '--basename foo'\n",
      "                         would result in the files 'foo.bam' and 'foo_SE_report.txt' (or its paired-end equivalent). Takes\n",
      "                         precedence over --prefix.\n",
      "\n",
      "--old_flag               Only in paired-end SAM mode, uses the FLAG values used by Bismark v0.8.2 and before. In addition,\n",
      "                         this options appends /1 and /2 to the read IDs for reads 1 and 2 relative to the input file. Since\n",
      "                         both the appended read IDs and custom FLAG values may cause problems with some downstream tools\n",
      "                         such as Picard, new defaults were implemented as of version 0.8.3.\n",
      "\n",
      "\n",
      "                                             default                         old_flag\n",
      "                                       ===================              ===================\n",
      "                                       Read 1       Read 2              Read 1       Read 2\n",
      "\n",
      "                              OT:         99          147                  67          131\n",
      "\n",
      "                              OB:         83          163                 115          179\n",
      "\n",
      "                              CTOT:      147           99                  67          131\n",
      "\n",
      "                              CTOB:      163           83                 115          179\n",
      "\n",
      "--ambig_bam              For reads that have multiple alignments a random alignment is written out to a special file ending in\n",
      "                         '.ambiguous.bam'. The alignments are in Bowtie2 format and do not any contain Bismark specific\n",
      "                         entries such as the methylation call etc. These ambiguous BAM files are intended to be used as\n",
      "                         coverage estimators for variant callers.\n",
      "\n",
      "--nucleotide_coverage    Calculates the mono- and di-nucleotide sequence composition of covered positions in the analysed BAM\n",
      "                         file and compares it to the genomic average composition once alignments are complete by calling 'bam2nuc'.\n",
      "                         Since this calculation may take a while, bam2nuc attempts to write the genomic sequence composition\n",
      "                         into a file called 'genomic_nucleotide_frequencies.txt' indside the reference genome folder so it can\n",
      "                         be re-used the next time round instead of calculating it once again. If a file 'nucleotide_stats.txt' is\n",
      "                         found with the Bismark reports it will be automatically detected and used for the Bismark HTML report.\n",
      "                         This option works only for BAM or CRAM files.\n",
      "\n",
      "\n",
      "Other:\n",
      "\n",
      "-h/--help                Displays this help file.\n",
      "\n",
      "-v/--version             Displays version information.\n",
      "\n",
      "\n",
      "BOWTIE 2 SPECIFIC OPTIONS\n",
      "\n",
      "--bowtie1                Uses Bowtie 1 instead of Bowtie 2, which might be a good choice for faster and very short\n",
      "                         alignments. Bismark assumes that raw sequence data is adapter and/or quality trimmed where\n",
      "                         appropriate. Default: off.\n",
      "\n",
      "--bowtie2                Default: ON. Uses Bowtie 2 instead of Bowtie 1. Bismark limits Bowtie 2 to only perform end-to-end\n",
      "                         alignments, i.e. searches for alignments involving all read characters (also called\n",
      "                         untrimmed or unclipped alignments). Bismark assumes that raw sequence data is adapter\n",
      "                         and/or quality trimmed where appropriate. Both small (.bt2) and large (.bt2l) Bowtie 2\n",
      "                         indexes are supported.\n",
      "\n",
      "Bowtie 2 alignment options:\n",
      "\n",
      "-N <int>                 Sets the number of mismatches to allowed in a seed alignment during multiseed alignment.\n",
      "                         Can be set to 0 or 1. Setting this higher makes alignment slower (often much slower)\n",
      "                         but increases sensitivity. Default: 0. This option is only available for Bowtie 2 (for\n",
      "                         Bowtie 1 see -n).\n",
      "\n",
      "-L <int>                 Sets the length of the seed substrings to align during multiseed alignment. Smaller values\n",
      "                         make alignment slower but more senstive. Default: the --sensitive preset of Bowtie 2 is\n",
      "                         used by default, which sets -L to 20. maximum of L can be set to 32. The length of the seed\n",
      "                         would effect the alignment speed dramatically while the larger L, the faster the aligment.\n",
      "                         This option is only available for Bowtie 2 (for Bowtie 1 see -l).\n",
      "\n",
      "--ignore-quals           When calculating a mismatch penalty, always consider the quality value at the mismatched\n",
      "                         position to be the highest possible, regardless of the actual value. I.e. input is treated\n",
      "                         as though all quality values are high. This is also the default behavior when the input\n",
      "                         doesn't specify quality values (e.g. in -f mode). This option is invariable and on by default.\n",
      "\n",
      "\n",
      "Bowtie 2 paired-end options:\n",
      "\n",
      "--no-mixed               This option disables Bowtie 2's behavior to try to find alignments for the individual mates if\n",
      "                         it cannot find a concordant or discordant alignment for a pair. This option is invariable and\n",
      "                         and on by default.\n",
      "\n",
      "--no-discordant          Normally, Bowtie 2 looks for discordant alignments if it cannot find any concordant alignments.\n",
      "                         A discordant alignment is an alignment where both mates align uniquely, but that does not\n",
      "                         satisfy the paired-end constraints (--fr/--rf/--ff, -I, -X). This option disables that behavior\n",
      "                         and it is on by default.\n",
      "\n",
      "--no_dovetail            It is possible, though unusual, for the mates to \"dovetail\", with the mates seemingly extending\n",
      "                         \"past\" each other as in this example:\n",
      "\n",
      "                         Mate 1:                 GTCAGCTACGATATTGTTTGGGGTGACACATTACGC\n",
      "                         Mate 2:            TATGAGTCAGCTACGATATTGTTTGGGGTGACACAT\n",
      "                         Reference: GCAGATTATATGAGTCAGCTACGATATTGTTTGGGGTGACACATTACGCGTCTTTGAC\n",
      "\n",
      "                         By default, dovetailing is considered inconsistent with concordant alignment, but by default\n",
      "                         Bismark calls Bowtie 2 with --dovetail, causing it to consider dovetailing alignments as\n",
      "                         concordant. This becomes relevant whenever reads are clipped from their 5' end prior to mapping,\n",
      "                         e.g. because of quality or bias issues.\n",
      "\n",
      "                         Specify --no_dovetail to turn off this behaviour for paired-end libraries. Default: OFF.\n",
      "\n",
      "\n",
      "Bowtie 2 effort options:\n",
      "\n",
      "-D <int>                 Up to <int> consecutive seed extension attempts can \"fail\" before Bowtie 2 moves on, using\n",
      "                         the alignments found so far. A seed extension \"fails\" if it does not yield a new best or a\n",
      "                         new second-best alignment. Default: 15.\n",
      "\n",
      "-R <int>                 <int> is the maximum number of times Bowtie 2 will \"re-seed\" reads with repetitive seeds.\n",
      "                         When \"re-seeding,\" Bowtie 2 simply chooses a new set of reads (same length, same number of\n",
      "                         mismatches allowed) at different offsets and searches for more alignments. A read is considered\n",
      "                         to have repetitive seeds if the total number of seed hits divided by the number of seeds\n",
      "                         that aligned at least once is greater than 300. Default: 2.\n",
      "\n",
      "Bowtie 2 parallelization options:\n",
      "\n",
      "\n",
      "-p NTHREADS              Launch NTHREADS parallel search threads (default: 1). Threads will run on separate processors/cores\n",
      "                         and synchronize when parsing reads and outputting alignments. Searching for alignments is highly\n",
      "                         parallel, and speedup is close to linear. Increasing -p increases Bowtie 2's memory footprint.\n",
      "                         E.g. when aligning to a human genome index, increasing -p from 1 to 8 increases the memory footprint\n",
      "                         by a few hundred megabytes. This option is only available if bowtie is linked with the pthreads\n",
      "                         library (i.e. if BOWTIE_PTHREADS=0 is not specified at build time). In addition, this option will\n",
      "                         automatically use the option '--reorder', which guarantees that output SAM records are printed in\n",
      "                         an order corresponding to the order of the reads in the original input file, even when -p is set\n",
      "                         greater than 1 (Bismark requires the Bowtie 2 output to be this way). Specifying --reorder and\n",
      "                         setting -p greater than 1 causes Bowtie 2 to run somewhat slower and use somewhat more memory then\n",
      "                         if --reorder were not specified. Has no effect if -p is set to 1, since output order will naturally\n",
      "                         correspond to input order in that case.\n",
      "\n",
      "Bowtie 2 Scoring options:\n",
      "\n",
      "--score_min <func>       Sets a function governing the minimum alignment score needed for an alignment to be considered\n",
      "                         \"valid\" (i.e. good enough to report). This is a function of read length. For instance, specifying\n",
      "                         L,0,-0.2 sets the minimum-score function f to f(x) = 0 + -0.2 * x, where x is the read length.\n",
      "                         See also: setting function options at http://bowtie-bio.sourceforge.net/bowtie2. The default is\n",
      "                         L,0,-0.2.\n",
      "\n",
      "--rdg <int1>,<int2>      Sets the read gap open (<int1>) and extend (<int2>) penalties. A read gap of length N gets a penalty\n",
      "                         of <int1> + N * <int2>. Default: 5, 3.\n",
      "\n",
      "--rfg <int1>,<int2>      Sets the reference gap open (<int1>) and extend (<int2>) penalties. A reference gap of length N gets\n",
      "                         a penalty of <int1> + N * <int2>. Default: 5, 3.\n",
      "\n",
      "\n",
      "Bowtie 2 Reporting options:\n",
      "\n",
      "-most_valid_alignments <int> This used to be the Bowtie 2 parameter -M. As of Bowtie 2 version 2.0.0 beta7 the option -M is\n",
      "                         deprecated. It will be removed in subsequent versions. What used to be called -M mode is still the\n",
      "                         default mode, but adjusting the -M setting is deprecated.  Use the -D and -R options to adjust the\n",
      "                         effort expended to find valid alignments.\n",
      "\n",
      "                         For reference, this used to be the old (now deprecated) description of -M:\n",
      "                         Bowtie 2 searches for at most <int>+1 distinct, valid alignments for each read. The search terminates when it\n",
      "                         can't find more distinct valid alignments, or when it finds <int>+1 distinct alignments, whichever\n",
      "                         happens first. Only the best alignment is reported. Information from the other alignments is used to\n",
      "                         estimate mapping quality and to set SAM optional fields, such as AS:i and XS:i. Increasing -M makes\n",
      "                         Bowtie 2 slower, but increases the likelihood that it will pick the correct alignment for a read that\n",
      "                         aligns many places. For reads that have more than <int>+1 distinct, valid alignments, Bowtie 2 does not\n",
      "                         guarantee that the alignment reported is the best possible in terms of alignment score. -M is\n",
      "                         always used and its default value is set to 10.\n",
      "\n",
      "\n",
      "'VANILLA' Bismark  OUTPUT:\n",
      "\n",
      "Single-end output format (tab-separated):\n",
      "\n",
      " (1) <seq-ID>\n",
      " (2) <read alignment strand>\n",
      " (3) <chromosome>\n",
      " (4) <start position>\n",
      " (5) <end position>\n",
      " (6) <observed bisulfite sequence>\n",
      " (7) <equivalent genomic sequence>\n",
      " (8) <methylation call>\n",
      " (9) <read conversion\n",
      "(10) <genome conversion>\n",
      "(11) <read quality score (Phred33)>\n",
      "\n",
      "\n",
      "Paired-end output format (tab-separated):\n",
      " (1) <seq-ID>\n",
      " (2) <read 1 alignment strand>\n",
      " (3) <chromosome>\n",
      " (4) <start position>\n",
      " (5) <end position>\n",
      " (6) <observed bisulfite sequence 1>\n",
      " (7) <equivalent genomic sequence 1>\n",
      " (8) <methylation call 1>\n",
      " (9) <observed bisulfite sequence 2>\n",
      "(10) <equivalent genomic sequence 2>\n",
      "(11) <methylation call 2>\n",
      "(12) <read 1 conversion\n",
      "(13) <genome conversion>\n",
      "(14) <read 1 quality score (Phred33)>\n",
      "(15) <read 2 quality score (Phred33)>\n",
      "\n",
      "\n",
      "Bismark SAM OUTPUT (default):\n",
      "\n",
      " (1) QNAME  (seq-ID)\n",
      " (2) FLAG   (this flag tries to take the strand a bisulfite read originated from into account (this is different from ordinary DNA alignment flags!))\n",
      " (3) RNAME  (chromosome)\n",
      " (4) POS    (start position)\n",
      " (5) MAPQ   (always 255 for use with Bowtie)\n",
      " (6) CIGAR\n",
      " (7) RNEXT\n",
      " (8) PNEXT\n",
      " (9) TLEN\n",
      "(10) SEQ\n",
      "(11) QUAL   (Phred33 scale)\n",
      "(12) NM-tag (edit distance to the reference)\n",
      "(13) MD-tag (base-by-base mismatches to the reference (handles indels)\n",
      "(14) XM-tag (methylation call string)\n",
      "(15) XR-tag (read conversion state for the alignment)\n",
      "(16) XG-tag (genome conversion state for the alignment)\n",
      "(17) XA/XB-tag (non-bisulfite mismatches) (optional!)\n",
      "\n",
      "Each read of paired-end alignments is written out in a separate line in the above format.\n",
      "\n",
      "\n",
      "Last modified on 10 October 2017\n"
     ]
    }
   ],
   "source": [
    "! ../../../../../Shared/Apps/Bismark_v0.19.0/bismark -help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " The trimmed *C. virginica* sequences I need to use are in [this OWL folder](http://owl.fish.washington.edu/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/). To run the `find` + `xargs` commands I want, I need to specify this path. I'm going to run this piece of code first just to confirm that I'm selecting the right files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_10_s1_R1_val_1.fq.gz\n",
      "/Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_1_s1_R1_val_1.fq.gz\n",
      "/Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_2_s1_R1_val_1.fq.gz\n",
      "/Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_3_s1_R1_val_1.fq.gz\n",
      "/Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_4_s1_R1_val_1.fq.gz\n",
      "/Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_5_s1_R1_val_1.fq.gz\n",
      "/Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_6_s1_R1_val_1.fq.gz\n",
      "/Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_7_s1_R1_val_1.fq.gz\n",
      "/Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_8_s1_R1_val_1.fq.gz\n",
      "/Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_9_s1_R1_val_1.fq.gz\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "find /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_*R1*.fq.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zr2096_10\n",
      "zr2096_1\n",
      "zr2096_2\n",
      "zr2096_3\n",
      "zr2096_4\n",
      "zr2096_5\n",
      "zr2096_6\n",
      "zr2096_7\n",
      "zr2096_8\n",
      "zr2096_9\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "find /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_*R1*.fq.gz \\\n",
    "| xargs basename -s _s1_R1_val_1.fq.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What I need for this command:\n",
    "\n",
    "1. Find all the *C. virginica* sequence files using `find`\n",
    "2. Isolate the base name from each of these files using `xargs` that I will then pipe into...\n",
    "\n",
    "*drumroll*\n",
    "\n",
    "1. ...The path to `bismark`\n",
    "2. --non_directional: See [this issue](https://github.com/RobertsLab/resources/issues/216)\n",
    "3. -p 4: Number of simultaneous threads to run\n",
    "4. --genome + path to the folder with the .fa genome, which also has all of the bisulfite genome directories\n",
    "5. -1 + Path to first paired file in the Athaliana folder\n",
    "6. -2 + Path to second paired file in the Athaliana folder\n",
    "7. Path to redirect standard error output.\n",
    "\n",
    "Note that I'm not including a `-score_min` option, since I want to just run this alignment with the default setting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FastQ format assumed (by default)\n",
      "Each Bowtie 2 instance is going to be run with 4 threads. Please monitor performance closely and tune down if needed!\n",
      "chr NC_035780.1 (65668440 bp)\n",
      "chr NC_035781.1 (61752955 bp)\n",
      "chr NC_035782.1 (77061148 bp)\n",
      "chr NC_035783.1 (59691872 bp)\n",
      "chr NC_035784.1 (98698416 bp)\n",
      "chr NC_035785.1 (51258098 bp)\n",
      "chr NC_035786.1 (57830854 bp)\n",
      "chr NC_035787.1 (75944018 bp)\n",
      "chr NC_035788.1 (104168038 bp)\n",
      "chr NC_035789.1 (32650045 bp)\n",
      "chr NC_007175.2 (17244 bp)\n",
      "\n",
      "Number of paired-end alignments with a unique best hit:\t3585012\n",
      "Mapping efficiency:\t20.5%\n",
      "\n",
      "Sequence pairs with no alignments under any condition:\t11875219\n",
      "Sequence pairs did not map uniquely:\t1988652\n",
      "Sequence pairs which were discarded because genomic sequence could not be extracted:\t0\n",
      "\n",
      "Number of sequence pairs with unique best (first) alignment came from the bowtie output:\n",
      "CT/GA/CT:\t910350\t((converted) top strand)\n",
      "GA/CT/CT:\t880499\t(complementary to (converted) top strand)\n",
      "GA/CT/GA:\t881086\t(complementary to (converted) bottom strand)\n",
      "CT/GA/GA:\t913077\t((converted) bottom strand)\n",
      "\n",
      "FastQ format assumed (by default)\n",
      "Each Bowtie 2 instance is going to be run with 4 threads. Please monitor performance closely and tune down if needed!\n",
      "chr NC_035780.1 (65668440 bp)\n",
      "chr NC_035781.1 (61752955 bp)\n",
      "chr NC_035782.1 (77061148 bp)\n",
      "chr NC_035783.1 (59691872 bp)\n",
      "chr NC_035784.1 (98698416 bp)\n",
      "chr NC_035785.1 (51258098 bp)\n",
      "chr NC_035786.1 (57830854 bp)\n",
      "chr NC_035787.1 (75944018 bp)\n",
      "chr NC_035788.1 (104168038 bp)\n",
      "chr NC_035789.1 (32650045 bp)\n",
      "chr NC_007175.2 (17244 bp)\n",
      "\n",
      "Number of paired-end alignments with a unique best hit:\t2299993\n",
      "Mapping efficiency:\t8.0%\n",
      "\n",
      "Sequence pairs with no alignments under any condition:\t24953333\n",
      "Sequence pairs did not map uniquely:\t1350020\n",
      "Sequence pairs which were discarded because genomic sequence could not be extracted:\t0\n",
      "\n",
      "Number of sequence pairs with unique best (first) alignment came from the bowtie output:\n",
      "CT/GA/CT:\t572719\t((converted) top strand)\n",
      "GA/CT/CT:\t575468\t(complementary to (converted) top strand)\n",
      "GA/CT/GA:\t578504\t(complementary to (converted) bottom strand)\n",
      "CT/GA/GA:\t573302\t((converted) bottom strand)\n",
      "\n",
      "FastQ format assumed (by default)\n",
      "Each Bowtie 2 instance is going to be run with 4 threads. Please monitor performance closely and tune down if needed!\n",
      "chr NC_035780.1 (65668440 bp)\n",
      "chr NC_035781.1 (61752955 bp)\n",
      "chr NC_035782.1 (77061148 bp)\n",
      "chr NC_035783.1 (59691872 bp)\n",
      "chr NC_035784.1 (98698416 bp)\n",
      "chr NC_035785.1 (51258098 bp)\n",
      "chr NC_035786.1 (57830854 bp)\n",
      "chr NC_035787.1 (75944018 bp)\n",
      "chr NC_035788.1 (104168038 bp)\n",
      "chr NC_035789.1 (32650045 bp)\n",
      "chr NC_007175.2 (17244 bp)\n",
      "\n",
      "Number of paired-end alignments with a unique best hit:\t5898474\n",
      "Mapping efficiency:\t19.5%\n",
      "\n",
      "Sequence pairs with no alignments under any condition:\t20848304\n",
      "Sequence pairs did not map uniquely:\t3578828\n",
      "Sequence pairs which were discarded because genomic sequence could not be extracted:\t1\n",
      "\n",
      "Number of sequence pairs with unique best (first) alignment came from the bowtie output:\n",
      "CT/GA/CT:\t1484169\t((converted) top strand)\n",
      "GA/CT/CT:\t1462149\t(complementary to (converted) top strand)\n",
      "GA/CT/GA:\t1465295\t(complementary to (converted) bottom strand)\n",
      "CT/GA/GA:\t1486860\t((converted) bottom strand)\n",
      "\n",
      "FastQ format assumed (by default)\n",
      "Each Bowtie 2 instance is going to be run with 4 threads. Please monitor performance closely and tune down if needed!\n",
      "chr NC_035780.1 (65668440 bp)\n",
      "chr NC_035781.1 (61752955 bp)\n",
      "chr NC_035782.1 (77061148 bp)\n",
      "chr NC_035783.1 (59691872 bp)\n",
      "chr NC_035784.1 (98698416 bp)\n",
      "chr NC_035785.1 (51258098 bp)\n",
      "chr NC_035786.1 (57830854 bp)\n",
      "chr NC_035787.1 (75944018 bp)\n",
      "chr NC_035788.1 (104168038 bp)\n",
      "chr NC_035789.1 (32650045 bp)\n",
      "chr NC_007175.2 (17244 bp)\n",
      "\n",
      "Number of paired-end alignments with a unique best hit:\t6143399\n",
      "Mapping efficiency:\t20.8%\n",
      "\n",
      "Sequence pairs with no alignments under any condition:\t19650427\n",
      "Sequence pairs did not map uniquely:\t3754927\n",
      "Sequence pairs which were discarded because genomic sequence could not be extracted:\t0\n",
      "\n",
      "Number of sequence pairs with unique best (first) alignment came from the bowtie output:\n",
      "CT/GA/CT:\t1528672\t((converted) top strand)\n",
      "GA/CT/CT:\t1538077\t(complementary to (converted) top strand)\n",
      "GA/CT/GA:\t1543167\t(complementary to (converted) bottom strand)\n",
      "CT/GA/GA:\t1533483\t((converted) bottom strand)\n",
      "\n",
      "FastQ format assumed (by default)\n",
      "Each Bowtie 2 instance is going to be run with 4 threads. Please monitor performance closely and tune down if needed!\n",
      "chr NC_035780.1 (65668440 bp)\n",
      "chr NC_035781.1 (61752955 bp)\n",
      "chr NC_035782.1 (77061148 bp)\n",
      "chr NC_035783.1 (59691872 bp)\n",
      "chr NC_035784.1 (98698416 bp)\n",
      "chr NC_035785.1 (51258098 bp)\n",
      "chr NC_035786.1 (57830854 bp)\n",
      "chr NC_035787.1 (75944018 bp)\n",
      "chr NC_035788.1 (104168038 bp)\n",
      "chr NC_035789.1 (32650045 bp)\n",
      "chr NC_007175.2 (17244 bp)\n",
      "\n",
      "Number of paired-end alignments with a unique best hit:\t5194561\n",
      "Mapping efficiency:\t21.7%\n",
      "\n",
      "Sequence pairs with no alignments under any condition:\t15690218\n",
      "Sequence pairs did not map uniquely:\t3085737\n",
      "Sequence pairs which were discarded because genomic sequence could not be extracted:\t0\n",
      "\n",
      "Number of sequence pairs with unique best (first) alignment came from the bowtie output:\n",
      "CT/GA/CT:\t1312800\t((converted) top strand)\n",
      "GA/CT/CT:\t1280122\t(complementary to (converted) top strand)\n",
      "GA/CT/GA:\t1282809\t(complementary to (converted) bottom strand)\n",
      "CT/GA/GA:\t1318830\t((converted) bottom strand)\n",
      "\n",
      "FastQ format assumed (by default)\n",
      "Each Bowtie 2 instance is going to be run with 4 threads. Please monitor performance closely and tune down if needed!\n",
      "chr NC_035780.1 (65668440 bp)\n",
      "chr NC_035781.1 (61752955 bp)\n",
      "chr NC_035782.1 (77061148 bp)\n",
      "chr NC_035783.1 (59691872 bp)\n",
      "chr NC_035784.1 (98698416 bp)\n",
      "chr NC_035785.1 (51258098 bp)\n",
      "chr NC_035786.1 (57830854 bp)\n",
      "chr NC_035787.1 (75944018 bp)\n",
      "chr NC_035788.1 (104168038 bp)\n",
      "chr NC_035789.1 (32650045 bp)\n",
      "chr NC_007175.2 (17244 bp)\n",
      "\n",
      "Number of paired-end alignments with a unique best hit:\t6213355\n",
      "Mapping efficiency:\t19.7%\n",
      "\n",
      "Sequence pairs with no alignments under any condition:\t21809843\n",
      "Sequence pairs did not map uniquely:\t3480083\n",
      "Sequence pairs which were discarded because genomic sequence could not be extracted:\t0\n",
      "\n",
      "Number of sequence pairs with unique best (first) alignment came from the bowtie output:\n",
      "CT/GA/CT:\t1583415\t((converted) top strand)\n",
      "GA/CT/CT:\t1528936\t(complementary to (converted) top strand)\n",
      "GA/CT/GA:\t1522571\t(complementary to (converted) bottom strand)\n",
      "CT/GA/GA:\t1578433\t((converted) bottom strand)\n",
      "\n",
      "FastQ format assumed (by default)\n",
      "Each Bowtie 2 instance is going to be run with 4 threads. Please monitor performance closely and tune down if needed!\n",
      "chr NC_035780.1 (65668440 bp)\n",
      "chr NC_035781.1 (61752955 bp)\n",
      "chr NC_035782.1 (77061148 bp)\n",
      "chr NC_035783.1 (59691872 bp)\n",
      "chr NC_035784.1 (98698416 bp)\n",
      "chr NC_035785.1 (51258098 bp)\n",
      "chr NC_035786.1 (57830854 bp)\n",
      "chr NC_035787.1 (75944018 bp)\n",
      "chr NC_035788.1 (104168038 bp)\n",
      "chr NC_035789.1 (32650045 bp)\n",
      "chr NC_007175.2 (17244 bp)\n",
      "\n",
      "Number of paired-end alignments with a unique best hit:\t5127457\n",
      "Mapping efficiency:\t21.4%\n",
      "\n",
      "Sequence pairs with no alignments under any condition:\t15685916\n",
      "Sequence pairs did not map uniquely:\t3096120\n",
      "Sequence pairs which were discarded because genomic sequence could not be extracted:\t0\n",
      "\n",
      "Number of sequence pairs with unique best (first) alignment came from the bowtie output:\n",
      "CT/GA/CT:\t1329957\t((converted) top strand)\n",
      "GA/CT/CT:\t1233621\t(complementary to (converted) top strand)\n",
      "GA/CT/GA:\t1233251\t(complementary to (converted) bottom strand)\n",
      "CT/GA/GA:\t1330628\t((converted) bottom strand)\n",
      "\n",
      "FastQ format assumed (by default)\n",
      "Each Bowtie 2 instance is going to be run with 4 threads. Please monitor performance closely and tune down if needed!\n",
      "chr NC_035780.1 (65668440 bp)\n",
      "chr NC_035781.1 (61752955 bp)\n",
      "chr NC_035782.1 (77061148 bp)\n",
      "chr NC_035783.1 (59691872 bp)\n",
      "chr NC_035784.1 (98698416 bp)\n",
      "chr NC_035785.1 (51258098 bp)\n",
      "chr NC_035786.1 (57830854 bp)\n",
      "chr NC_035787.1 (75944018 bp)\n",
      "chr NC_035788.1 (104168038 bp)\n",
      "chr NC_035789.1 (32650045 bp)\n",
      "chr NC_007175.2 (17244 bp)\n",
      "\n",
      "Number of paired-end alignments with a unique best hit:\t5560803\n",
      "Mapping efficiency:\t19.0%\n",
      "\n",
      "Sequence pairs with no alignments under any condition:\t20689615\n",
      "Sequence pairs did not map uniquely:\t3023217\n",
      "Sequence pairs which were discarded because genomic sequence could not be extracted:\t0\n",
      "\n",
      "Number of sequence pairs with unique best (first) alignment came from the bowtie output:\n",
      "CT/GA/CT:\t1412253\t((converted) top strand)\n",
      "GA/CT/CT:\t1368227\t(complementary to (converted) top strand)\n",
      "GA/CT/GA:\t1366318\t(complementary to (converted) bottom strand)\n",
      "CT/GA/GA:\t1414005\t((converted) bottom strand)\n",
      "\n",
      "FastQ format assumed (by default)\n",
      "Each Bowtie 2 instance is going to be run with 4 threads. Please monitor performance closely and tune down if needed!\n",
      "chr NC_035780.1 (65668440 bp)\n",
      "chr NC_035781.1 (61752955 bp)\n",
      "chr NC_035782.1 (77061148 bp)\n",
      "chr NC_035783.1 (59691872 bp)\n",
      "chr NC_035784.1 (98698416 bp)\n",
      "chr NC_035785.1 (51258098 bp)\n",
      "chr NC_035786.1 (57830854 bp)\n",
      "chr NC_035787.1 (75944018 bp)\n",
      "chr NC_035788.1 (104168038 bp)\n",
      "chr NC_035789.1 (32650045 bp)\n",
      "chr NC_007175.2 (17244 bp)\n",
      "\n",
      "Number of paired-end alignments with a unique best hit:\t5162019\n",
      "Mapping efficiency:\t17.5%\n",
      "\n",
      "Sequence pairs with no alignments under any condition:\t21488196\n",
      "Sequence pairs did not map uniquely:\t2833003\n",
      "Sequence pairs which were discarded because genomic sequence could not be extracted:\t0\n",
      "\n",
      "Number of sequence pairs with unique best (first) alignment came from the bowtie output:\n",
      "CT/GA/CT:\t1320745\t((converted) top strand)\n",
      "GA/CT/CT:\t1271713\t(complementary to (converted) top strand)\n",
      "GA/CT/GA:\t1258935\t(complementary to (converted) bottom strand)\n",
      "CT/GA/GA:\t1310626\t((converted) bottom strand)\n",
      "\n",
      "FastQ format assumed (by default)\n",
      "Each Bowtie 2 instance is going to be run with 4 threads. Please monitor performance closely and tune down if needed!\n",
      "chr NC_035780.1 (65668440 bp)\n",
      "chr NC_035781.1 (61752955 bp)\n",
      "chr NC_035782.1 (77061148 bp)\n",
      "chr NC_035783.1 (59691872 bp)\n",
      "chr NC_035784.1 (98698416 bp)\n",
      "chr NC_035785.1 (51258098 bp)\n",
      "chr NC_035786.1 (57830854 bp)\n",
      "chr NC_035787.1 (75944018 bp)\n",
      "chr NC_035788.1 (104168038 bp)\n",
      "chr NC_035789.1 (32650045 bp)\n",
      "chr NC_007175.2 (17244 bp)\n",
      "\n",
      "Number of paired-end alignments with a unique best hit:\t6030681\n",
      "Mapping efficiency:\t18.9%\n",
      "\n",
      "Sequence pairs with no alignments under any condition:\t22269696\n",
      "Sequence pairs did not map uniquely:\t3547164\n",
      "Sequence pairs which were discarded because genomic sequence could not be extracted:\t0\n",
      "\n",
      "Number of sequence pairs with unique best (first) alignment came from the bowtie output:\n",
      "CT/GA/CT:\t1609858\t((converted) top strand)\n",
      "GA/CT/CT:\t1407149\t(complementary to (converted) top strand)\n",
      "GA/CT/GA:\t1404003\t(complementary to (converted) bottom strand)\n",
      "CT/GA/GA:\t1609671\t((converted) bottom strand)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "\n",
    "find /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_*R1*.fq.gz \\\n",
    "| xargs basename -s _s1_R1_val_1.fq.gz | xargs -I{} /Users/Shared/Apps/Bismark_v0.19.0/bismark \\\n",
    "--non_directional \\\n",
    "-p 4 \\\n",
    "--genome /Users/yaamini/Documents/project-virginica-oa/analyses/2018-04-27-Bismark/2018-04-27-Bismark-Inputs/ \\\n",
    "-1 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/{}_s1_R1_val_1.fq.gz \\\n",
    "-2 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/{}_s1_R2_val_2.fq.gz \\\n",
    "2> bismark-full-0522.err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bismark-full-0522.err\r\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.bam\r\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.bam\r\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.bam\r\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.bam\r\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.bam\r\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.bam\r\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.bam\r\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.bam\r\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.bam\r\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.bam\r\n"
     ]
    }
   ],
   "source": [
    "!ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The alignment is done! Now to deduplication."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Deduplication, Sorting, and Indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3a. Deduplication"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This step wasn't explicitly laid out in the User Guide. The deduplication step removes any PCR artifacts from the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n",
      "===============================================================================================================\r\n",
      "\r\n",
      "This script is supposed to remove alignments to the same position in the genome from the Bismark mapping output\r\n",
      "(both single and paired-end SAM files), which can arise by e.g. excessive PCR amplification. If sequences align\r\n",
      "to the same genomic position but on different strands they will be scored individually.\r\n",
      "\r\n",
      "Note that deduplication is not recommended for RRBS-type experiments!\r\n",
      "\r\n",
      "In the default mode, the first alignment to a given position will be used irrespective of its methylation call\r\n",
      "(this is the fastest option, and as the alignments are not ordered in any way this is also near enough random).\r\n",
      "\r\n",
      "For single-end alignments only use the start coordinate of a read will be used for deduplication.\r\n",
      "\r\n",
      "For paired-end alignments the start-coordinate of the first read and the end coordinate of the second\r\n",
      "read will be used for deduplication. This script expects the Bismark output to be in SAM format\r\n",
      "(Bismark v0.6.x or higher). To deduplicate the old custom Bismark output please specify '--vanilla'.\r\n",
      "\r\n",
      "*** Please note that for paired-end BAM files the deduplication script expects Read1 and Read2 to\r\n",
      "follow each other in consecutive lines! If the file has been sorted by position make sure that you resort it\r\n",
      "by read name first (e.g. using samtools sort -n)  ***\r\n",
      "\r\n",
      "===============================================================================================================\r\n",
      "\r\n",
      ">>> USAGE: ./deduplicate_bismark_alignment_output.pl [options] filename(s) <<<\r\n",
      "\r\n",
      "-s/--single\t\tdeduplicate single-end Bismark files (default format: SAM)\r\n",
      "-p/--paired\t\tdeduplicate paired-end Bismark files (default format: SAM)\r\n",
      "\r\n",
      "--output_dir [path]\tOutput directory, either relative or absolute. Output is written to the current directory if not\r\n",
      "\t\t\tspecified explicitly.\r\n",
      "\r\n",
      "--vanilla\t\tThe input file is in the old custom Bismark format and not in SAM format\r\n",
      "\r\n",
      "--barcode\t\tIn addition to chromosome, start position and orientation this will also take a potential barcode into\r\n",
      "                        consideration while deduplicating. The barcode needs to be the last element of the read ID and separated\r\n",
      "                        by a ':', e.g.: MISEQ:14:000000000-A55D0:1:1101:18024:2858_1:N:0:CTCCT\r\n",
      "\r\n",
      "--bam\t\t\tThe output will be written out in BAM format instead of the default SAM format. This script will\r\n",
      "\t\t\tattempt to use the path to Samtools that was specified with '--samtools_path', or, if it hasn't\r\n",
      "\t\t\tbeen specified, attempt to find Samtools in the PATH. If no installation of Samtools can be found,\r\n",
      "\t\t\tthe SAM output will be compressed with GZIP instead (yielding a .sam.gz output file)\r\n",
      "\r\n",
      "--multiple\t\tAll specified input files are treated as one sample and concatenated together for deduplication.\r\n",
      "\t\t\tThis uses Unix 'cat' for SAM files and 'samtools cat' for BAM files. Additional notes for BAM files:\r\n",
      "\t\t\tAlthough this works on either BAM or CRAM, all input files must be the same format as each other.\r\n",
      "\t\t\tThe sequence dictionary of each input file must be identical, although this command does not check this.\r\n",
      "\t\t\tBy default the header is taken from the first file to be concatenated.\r\n",
      "\r\n",
      "--samtools_path\t\tThe path to your Samtools installation, e.g. /home/user/samtools/. Does not need to be specified\r\n",
      "\t\t\texplicitly if Samtools is in the PATH already\r\n",
      "\r\n",
      "--version\t\tPrint version information and exit\r\n",
      "===============================================================================================================\r\n",
      "\r\n",
      "This script was last modified on 07 September 2017\r\n",
      "\r\n"
     ]
    }
   ],
   "source": [
    "! ../../../../../Shared/Apps/Bismark_v0.19.0/deduplicate_bismark -help"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.bam\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.bam\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.bam\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.bam\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.bam\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.bam\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.bam\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.bam\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.bam\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.bam\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "\n",
    "find zr2096_*_s1_R1_val_1_bismark_bt2_pe.bam \\"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what I need to run this command:\n",
    "\n",
    "1. Path to `deduplicate_bismark`\n",
    "2. -p: Input data is paired\n",
    "3. --bam: Write the output as a .bam file\n",
    "4. Path to input files. It's important to specify only the files from the alignment! If I just used .bam, the process would become recursive, and the command will start to deduplicate the .deduplicated.bam files that are already deduplicated!\n",
    "5. Path to standard error output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Now testing Bismark result file zr2096_10_s1_R1_val_1_bismark_bt2_pe.bam for positional sorting (which would be bad...)\tskipping header line:\t@HD\tVN:1.0\tSO:unsorted\n",
      "skipping header line:\t@SQ\tSN:NC_035780.1\tLN:65668440\n",
      "skipping header line:\t@SQ\tSN:NC_035781.1\tLN:61752955\n",
      "skipping header line:\t@SQ\tSN:NC_035782.1\tLN:77061148\n",
      "skipping header line:\t@SQ\tSN:NC_035783.1\tLN:59691872\n",
      "skipping header line:\t@SQ\tSN:NC_035784.1\tLN:98698416\n",
      "skipping header line:\t@SQ\tSN:NC_035785.1\tLN:51258098\n",
      "skipping header line:\t@SQ\tSN:NC_035786.1\tLN:57830854\n",
      "skipping header line:\t@SQ\tSN:NC_035787.1\tLN:75944018\n",
      "skipping header line:\t@SQ\tSN:NC_035788.1\tLN:104168038\n",
      "skipping header line:\t@SQ\tSN:NC_035789.1\tLN:32650045\n",
      "skipping header line:\t@SQ\tSN:NC_007175.2\tLN:17244\n",
      "skipping header line:\t@PG\tID:Bismark\tVN:v0.19.0\tCL:\"bismark --non_directional -p 4 --genome /Users/yaamini/Documents/project-virginica-oa/analyses/2018-04-27-Bismark/2018-04-27-Bismark-Inputs/ -1 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_10_s1_R1_val_1.fq.gz -2 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_10_s1_R2_val_2.fq.gz\"\n",
      "\n",
      "Now testing Bismark result file zr2096_1_s1_R1_val_1_bismark_bt2_pe.bam for positional sorting (which would be bad...)\tskipping header line:\t@HD\tVN:1.0\tSO:unsorted\n",
      "skipping header line:\t@SQ\tSN:NC_035780.1\tLN:65668440\n",
      "skipping header line:\t@SQ\tSN:NC_035781.1\tLN:61752955\n",
      "skipping header line:\t@SQ\tSN:NC_035782.1\tLN:77061148\n",
      "skipping header line:\t@SQ\tSN:NC_035783.1\tLN:59691872\n",
      "skipping header line:\t@SQ\tSN:NC_035784.1\tLN:98698416\n",
      "skipping header line:\t@SQ\tSN:NC_035785.1\tLN:51258098\n",
      "skipping header line:\t@SQ\tSN:NC_035786.1\tLN:57830854\n",
      "skipping header line:\t@SQ\tSN:NC_035787.1\tLN:75944018\n",
      "skipping header line:\t@SQ\tSN:NC_035788.1\tLN:104168038\n",
      "skipping header line:\t@SQ\tSN:NC_035789.1\tLN:32650045\n",
      "skipping header line:\t@SQ\tSN:NC_007175.2\tLN:17244\n",
      "skipping header line:\t@PG\tID:Bismark\tVN:v0.19.0\tCL:\"bismark --non_directional -p 4 --genome /Users/yaamini/Documents/project-virginica-oa/analyses/2018-04-27-Bismark/2018-04-27-Bismark-Inputs/ -1 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_1_s1_R1_val_1.fq.gz -2 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_1_s1_R2_val_2.fq.gz\"\n",
      "\n",
      "Now testing Bismark result file zr2096_2_s1_R1_val_1_bismark_bt2_pe.bam for positional sorting (which would be bad...)\tskipping header line:\t@HD\tVN:1.0\tSO:unsorted\n",
      "skipping header line:\t@SQ\tSN:NC_035780.1\tLN:65668440\n",
      "skipping header line:\t@SQ\tSN:NC_035781.1\tLN:61752955\n",
      "skipping header line:\t@SQ\tSN:NC_035782.1\tLN:77061148\n",
      "skipping header line:\t@SQ\tSN:NC_035783.1\tLN:59691872\n",
      "skipping header line:\t@SQ\tSN:NC_035784.1\tLN:98698416\n",
      "skipping header line:\t@SQ\tSN:NC_035785.1\tLN:51258098\n",
      "skipping header line:\t@SQ\tSN:NC_035786.1\tLN:57830854\n",
      "skipping header line:\t@SQ\tSN:NC_035787.1\tLN:75944018\n",
      "skipping header line:\t@SQ\tSN:NC_035788.1\tLN:104168038\n",
      "skipping header line:\t@SQ\tSN:NC_035789.1\tLN:32650045\n",
      "skipping header line:\t@SQ\tSN:NC_007175.2\tLN:17244\n",
      "skipping header line:\t@PG\tID:Bismark\tVN:v0.19.0\tCL:\"bismark --non_directional -p 4 --genome /Users/yaamini/Documents/project-virginica-oa/analyses/2018-04-27-Bismark/2018-04-27-Bismark-Inputs/ -1 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_2_s1_R1_val_1.fq.gz -2 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_2_s1_R2_val_2.fq.gz\"\n",
      "\n",
      "Now testing Bismark result file zr2096_3_s1_R1_val_1_bismark_bt2_pe.bam for positional sorting (which would be bad...)\tskipping header line:\t@HD\tVN:1.0\tSO:unsorted\n",
      "skipping header line:\t@SQ\tSN:NC_035780.1\tLN:65668440\n",
      "skipping header line:\t@SQ\tSN:NC_035781.1\tLN:61752955\n",
      "skipping header line:\t@SQ\tSN:NC_035782.1\tLN:77061148\n",
      "skipping header line:\t@SQ\tSN:NC_035783.1\tLN:59691872\n",
      "skipping header line:\t@SQ\tSN:NC_035784.1\tLN:98698416\n",
      "skipping header line:\t@SQ\tSN:NC_035785.1\tLN:51258098\n",
      "skipping header line:\t@SQ\tSN:NC_035786.1\tLN:57830854\n",
      "skipping header line:\t@SQ\tSN:NC_035787.1\tLN:75944018\n",
      "skipping header line:\t@SQ\tSN:NC_035788.1\tLN:104168038\n",
      "skipping header line:\t@SQ\tSN:NC_035789.1\tLN:32650045\n",
      "skipping header line:\t@SQ\tSN:NC_007175.2\tLN:17244\n",
      "skipping header line:\t@PG\tID:Bismark\tVN:v0.19.0\tCL:\"bismark --non_directional -p 4 --genome /Users/yaamini/Documents/project-virginica-oa/analyses/2018-04-27-Bismark/2018-04-27-Bismark-Inputs/ -1 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_3_s1_R1_val_1.fq.gz -2 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_3_s1_R2_val_2.fq.gz\"\n",
      "\n",
      "Now testing Bismark result file zr2096_4_s1_R1_val_1_bismark_bt2_pe.bam for positional sorting (which would be bad...)\tskipping header line:\t@HD\tVN:1.0\tSO:unsorted\n",
      "skipping header line:\t@SQ\tSN:NC_035780.1\tLN:65668440\n",
      "skipping header line:\t@SQ\tSN:NC_035781.1\tLN:61752955\n",
      "skipping header line:\t@SQ\tSN:NC_035782.1\tLN:77061148\n",
      "skipping header line:\t@SQ\tSN:NC_035783.1\tLN:59691872\n",
      "skipping header line:\t@SQ\tSN:NC_035784.1\tLN:98698416\n",
      "skipping header line:\t@SQ\tSN:NC_035785.1\tLN:51258098\n",
      "skipping header line:\t@SQ\tSN:NC_035786.1\tLN:57830854\n",
      "skipping header line:\t@SQ\tSN:NC_035787.1\tLN:75944018\n",
      "skipping header line:\t@SQ\tSN:NC_035788.1\tLN:104168038\n",
      "skipping header line:\t@SQ\tSN:NC_035789.1\tLN:32650045\n",
      "skipping header line:\t@SQ\tSN:NC_007175.2\tLN:17244\n",
      "skipping header line:\t@PG\tID:Bismark\tVN:v0.19.0\tCL:\"bismark --non_directional -p 4 --genome /Users/yaamini/Documents/project-virginica-oa/analyses/2018-04-27-Bismark/2018-04-27-Bismark-Inputs/ -1 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_4_s1_R1_val_1.fq.gz -2 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_4_s1_R2_val_2.fq.gz\"\n",
      "\n",
      "Now testing Bismark result file zr2096_5_s1_R1_val_1_bismark_bt2_pe.bam for positional sorting (which would be bad...)\tskipping header line:\t@HD\tVN:1.0\tSO:unsorted\n",
      "skipping header line:\t@SQ\tSN:NC_035780.1\tLN:65668440\n",
      "skipping header line:\t@SQ\tSN:NC_035781.1\tLN:61752955\n",
      "skipping header line:\t@SQ\tSN:NC_035782.1\tLN:77061148\n",
      "skipping header line:\t@SQ\tSN:NC_035783.1\tLN:59691872\n",
      "skipping header line:\t@SQ\tSN:NC_035784.1\tLN:98698416\n",
      "skipping header line:\t@SQ\tSN:NC_035785.1\tLN:51258098\n",
      "skipping header line:\t@SQ\tSN:NC_035786.1\tLN:57830854\n",
      "skipping header line:\t@SQ\tSN:NC_035787.1\tLN:75944018\n",
      "skipping header line:\t@SQ\tSN:NC_035788.1\tLN:104168038\n",
      "skipping header line:\t@SQ\tSN:NC_035789.1\tLN:32650045\n",
      "skipping header line:\t@SQ\tSN:NC_007175.2\tLN:17244\n",
      "skipping header line:\t@PG\tID:Bismark\tVN:v0.19.0\tCL:\"bismark --non_directional -p 4 --genome /Users/yaamini/Documents/project-virginica-oa/analyses/2018-04-27-Bismark/2018-04-27-Bismark-Inputs/ -1 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_5_s1_R1_val_1.fq.gz -2 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_5_s1_R2_val_2.fq.gz\"\n",
      "\n",
      "Now testing Bismark result file zr2096_6_s1_R1_val_1_bismark_bt2_pe.bam for positional sorting (which would be bad...)\tskipping header line:\t@HD\tVN:1.0\tSO:unsorted\n",
      "skipping header line:\t@SQ\tSN:NC_035780.1\tLN:65668440\n",
      "skipping header line:\t@SQ\tSN:NC_035781.1\tLN:61752955\n",
      "skipping header line:\t@SQ\tSN:NC_035782.1\tLN:77061148\n",
      "skipping header line:\t@SQ\tSN:NC_035783.1\tLN:59691872\n",
      "skipping header line:\t@SQ\tSN:NC_035784.1\tLN:98698416\n",
      "skipping header line:\t@SQ\tSN:NC_035785.1\tLN:51258098\n",
      "skipping header line:\t@SQ\tSN:NC_035786.1\tLN:57830854\n",
      "skipping header line:\t@SQ\tSN:NC_035787.1\tLN:75944018\n",
      "skipping header line:\t@SQ\tSN:NC_035788.1\tLN:104168038\n",
      "skipping header line:\t@SQ\tSN:NC_035789.1\tLN:32650045\n",
      "skipping header line:\t@SQ\tSN:NC_007175.2\tLN:17244\n",
      "skipping header line:\t@PG\tID:Bismark\tVN:v0.19.0\tCL:\"bismark --non_directional -p 4 --genome /Users/yaamini/Documents/project-virginica-oa/analyses/2018-04-27-Bismark/2018-04-27-Bismark-Inputs/ -1 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_6_s1_R1_val_1.fq.gz -2 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_6_s1_R2_val_2.fq.gz\"\n",
      "\n",
      "Now testing Bismark result file zr2096_7_s1_R1_val_1_bismark_bt2_pe.bam for positional sorting (which would be bad...)\tskipping header line:\t@HD\tVN:1.0\tSO:unsorted\n",
      "skipping header line:\t@SQ\tSN:NC_035780.1\tLN:65668440\n",
      "skipping header line:\t@SQ\tSN:NC_035781.1\tLN:61752955\n",
      "skipping header line:\t@SQ\tSN:NC_035782.1\tLN:77061148\n",
      "skipping header line:\t@SQ\tSN:NC_035783.1\tLN:59691872\n",
      "skipping header line:\t@SQ\tSN:NC_035784.1\tLN:98698416\n",
      "skipping header line:\t@SQ\tSN:NC_035785.1\tLN:51258098\n",
      "skipping header line:\t@SQ\tSN:NC_035786.1\tLN:57830854\n",
      "skipping header line:\t@SQ\tSN:NC_035787.1\tLN:75944018\n",
      "skipping header line:\t@SQ\tSN:NC_035788.1\tLN:104168038\n",
      "skipping header line:\t@SQ\tSN:NC_035789.1\tLN:32650045\n",
      "skipping header line:\t@SQ\tSN:NC_007175.2\tLN:17244\n",
      "skipping header line:\t@PG\tID:Bismark\tVN:v0.19.0\tCL:\"bismark --non_directional -p 4 --genome /Users/yaamini/Documents/project-virginica-oa/analyses/2018-04-27-Bismark/2018-04-27-Bismark-Inputs/ -1 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_7_s1_R1_val_1.fq.gz -2 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_7_s1_R2_val_2.fq.gz\"\n",
      "\n",
      "Now testing Bismark result file zr2096_8_s1_R1_val_1_bismark_bt2_pe.bam for positional sorting (which would be bad...)\tskipping header line:\t@HD\tVN:1.0\tSO:unsorted\n",
      "skipping header line:\t@SQ\tSN:NC_035780.1\tLN:65668440\n",
      "skipping header line:\t@SQ\tSN:NC_035781.1\tLN:61752955\n",
      "skipping header line:\t@SQ\tSN:NC_035782.1\tLN:77061148\n",
      "skipping header line:\t@SQ\tSN:NC_035783.1\tLN:59691872\n",
      "skipping header line:\t@SQ\tSN:NC_035784.1\tLN:98698416\n",
      "skipping header line:\t@SQ\tSN:NC_035785.1\tLN:51258098\n",
      "skipping header line:\t@SQ\tSN:NC_035786.1\tLN:57830854\n",
      "skipping header line:\t@SQ\tSN:NC_035787.1\tLN:75944018\n",
      "skipping header line:\t@SQ\tSN:NC_035788.1\tLN:104168038\n",
      "skipping header line:\t@SQ\tSN:NC_035789.1\tLN:32650045\n",
      "skipping header line:\t@SQ\tSN:NC_007175.2\tLN:17244\n",
      "skipping header line:\t@PG\tID:Bismark\tVN:v0.19.0\tCL:\"bismark --non_directional -p 4 --genome /Users/yaamini/Documents/project-virginica-oa/analyses/2018-04-27-Bismark/2018-04-27-Bismark-Inputs/ -1 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_8_s1_R1_val_1.fq.gz -2 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_8_s1_R2_val_2.fq.gz\"\n",
      "\n",
      "Now testing Bismark result file zr2096_9_s1_R1_val_1_bismark_bt2_pe.bam for positional sorting (which would be bad...)\tskipping header line:\t@HD\tVN:1.0\tSO:unsorted\n",
      "skipping header line:\t@SQ\tSN:NC_035780.1\tLN:65668440\n",
      "skipping header line:\t@SQ\tSN:NC_035781.1\tLN:61752955\n",
      "skipping header line:\t@SQ\tSN:NC_035782.1\tLN:77061148\n",
      "skipping header line:\t@SQ\tSN:NC_035783.1\tLN:59691872\n",
      "skipping header line:\t@SQ\tSN:NC_035784.1\tLN:98698416\n",
      "skipping header line:\t@SQ\tSN:NC_035785.1\tLN:51258098\n",
      "skipping header line:\t@SQ\tSN:NC_035786.1\tLN:57830854\n",
      "skipping header line:\t@SQ\tSN:NC_035787.1\tLN:75944018\n",
      "skipping header line:\t@SQ\tSN:NC_035788.1\tLN:104168038\n",
      "skipping header line:\t@SQ\tSN:NC_035789.1\tLN:32650045\n",
      "skipping header line:\t@SQ\tSN:NC_007175.2\tLN:17244\n",
      "skipping header line:\t@PG\tID:Bismark\tVN:v0.19.0\tCL:\"bismark --non_directional -p 4 --genome /Users/yaamini/Documents/project-virginica-oa/analyses/2018-04-27-Bismark/2018-04-27-Bismark-Inputs/ -1 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_9_s1_R1_val_1.fq.gz -2 /Volumes/web/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/zr2096_9_s1_R2_val_2.fq.gz\"\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "\n",
    "/Users/Shared/Apps/Bismark_v0.19.0/deduplicate_bismark \\\n",
    "-p \\\n",
    "--bam \\\n",
    "zr2096_*_s1_R1_val_1_bismark_bt2_pe.bam \\\n",
    "2> dedup-full-0524.err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dedup-full-0524.err\r\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt\r\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt\r\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt\r\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt\r\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt\r\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt\r\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt\r\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt\r\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt\r\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt\r\n"
     ]
    }
   ],
   "source": [
    "!find *dedup*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam\r\n"
     ]
    }
   ],
   "source": [
    "!find *deduplicated.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 3b. Sorting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next two steps will allow me to use my .bam files in IGV. According to the [`igvtools`](https://software.broadinstitute.org/software/igv/igvtools_commandline) website, I need to `sort` then `index` all of my .bam files. Instead of `igvtools`, I will use `samtools`, since `igvtools` may not have the capacity to deal with paired-end read data (ie. my data).\n",
    "\n",
    "I will use a similar `find` and `xargs` approach to loop through all of my files.\n",
    "\n",
    "1. Identify all deduplicated .bam files using `find`\n",
    "2. Isolate the basename with `xargs` to be fed into `igvtools`\n",
    "\n",
    "THEN\n",
    "\n",
    "1. Path to `samtools`\n",
    "2. `sort`\n",
    "3. Path to input .bam file\n",
    "4. -o + Path to output file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[bam_sort_core] merging from 3 files and 1 in-memory blocks...\n",
      "[bam_sort_core] merging from 1 files and 1 in-memory blocks...\n",
      "[bam_sort_core] merging from 4 files and 1 in-memory blocks...\n",
      "[bam_sort_core] merging from 5 files and 1 in-memory blocks...\n",
      "[bam_sort_core] merging from 4 files and 1 in-memory blocks...\n",
      "[bam_sort_core] merging from 5 files and 1 in-memory blocks...\n",
      "[bam_sort_core] merging from 4 files and 1 in-memory blocks...\n",
      "[bam_sort_core] merging from 4 files and 1 in-memory blocks...\n",
      "[bam_sort_core] merging from 4 files and 1 in-memory blocks...\n",
      "[bam_sort_core] merging from 4 files and 1 in-memory blocks...\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "find *deduplicated.bam \\\n",
    "| xargs basename -s _s1_R1_val_1_bismark_bt2_pe.deduplicated.bam | xargs -I{} /Users/Shared/Apps/samtools-1.8/samtools \\\n",
    "sort {}_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam \\\n",
    "-o {}_dedup.sorted.bam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zr2096_10_dedup.sorted.bam\r\n",
      "zr2096_1_dedup.sorted.bam\r\n",
      "zr2096_2_dedup.sorted.bam\r\n",
      "zr2096_3_dedup.sorted.bam\r\n",
      "zr2096_4_dedup.sorted.bam\r\n",
      "zr2096_5_dedup.sorted.bam\r\n",
      "zr2096_6_dedup.sorted.bam\r\n",
      "zr2096_7_dedup.sorted.bam\r\n",
      "zr2096_8_dedup.sorted.bam\r\n",
      "zr2096_9_dedup.sorted.bam\r\n"
     ]
    }
   ],
   "source": [
    "!find *dedup.sorted.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3c. Indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`sort` worked! Now to use `index`.\n",
    "\n",
    "1. Identify all deduplicated .bam files using `find`\n",
    "2. Isolate the basename with `xargs` to be fed into `igvtools`\n",
    "\n",
    "THEN\n",
    "\n",
    "1. Path to `samtools`\n",
    "2. `index`\n",
    "3. Path to input .bam file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "find *dedup.sorted.bam \\\n",
    "| xargs basename -s _dedup.sorted.bam | xargs -I{} ../../../../../Shared/Apps/samtools-1.8/samtools \\\n",
    "index {}_dedup.sorted.bam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zr2096_10_dedup.sorted.bam.bai\r\n",
      "zr2096_1_dedup.sorted.bam.bai\r\n",
      "zr2096_2_dedup.sorted.bam.bai\r\n",
      "zr2096_3_dedup.sorted.bam.bai\r\n",
      "zr2096_4_dedup.sorted.bam.bai\r\n",
      "zr2096_5_dedup.sorted.bam.bai\r\n",
      "zr2096_6_dedup.sorted.bam.bai\r\n",
      "zr2096_7_dedup.sorted.bam.bai\r\n",
      "zr2096_8_dedup.sorted.bam.bai\r\n",
      "zr2096_9_dedup.sorted.bam.bai\r\n"
     ]
    }
   ],
   "source": [
    "!find *.bai"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sorted files need to be used in `methylKit`, while the sorted and indexed files are for IGV. Before I delve into these analayses, I will finish using `bismark`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Methylation Extractor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "DESCRIPTION\n",
      "\n",
      "The following is a brief description of all options to control the Bismark\n",
      "methylation extractor. The script reads in a bisulfite read alignment results file \n",
      "produced by the Bismark bisulfite mapper (in BAM/CRAM/SAM format) and extracts the\n",
      "methylation information for individual cytosines. This information is found in the\n",
      "methylation call field which can contain the following characters:\n",
      "\n",
      "       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
      "       ~~~   X   for methylated C in CHG context                      ~~~\n",
      "       ~~~   x   for not methylated C CHG                             ~~~\n",
      "       ~~~   H   for methylated C in CHH context                      ~~~\n",
      "       ~~~   h   for not methylated C in CHH context                  ~~~\n",
      "       ~~~   Z   for methylated C in CpG context                      ~~~\n",
      "       ~~~   z   for not methylated C in CpG context                  ~~~\n",
      "       ~~~   U   for methylated C in Unknown context (CN or CHN       ~~~\n",
      "       ~~~   u   for not methylated C in Unknown context (CN or CHN)  ~~~\n",
      "       ~~~   .   for any bases not involving cytosines                ~~~\n",
      "       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
      "\n",
      "The methylation extractor outputs result files for cytosines in CpG, CHG and CHH\n",
      "context (this distinction is actually already made in Bismark itself). As the methylation\n",
      "information for every C analysed can produce files which easily have tens or even hundreds of\n",
      "millions of lines, file sizes can become very large and more difficult to handle. The C\n",
      "methylation info additionally splits cytosine methylation calls up into one of the four possible\n",
      "strands a given bisulfite read aligned against:\n",
      "\n",
      "             OT      original top strand\n",
      "             CTOT    complementary to original top strand\n",
      "\n",
      "             OB      original bottom strand\n",
      "             CTOB    complementary to original bottom strand\n",
      "\n",
      "Thus, by default twelve individual output files are being generated per input file (unless\n",
      "--comprehensive is specified, see below). The output files can be imported into a genome\n",
      "viewer, such as SeqMonk, and re-combined into a single data group if desired (in fact\n",
      "unless the bisulfite reads were generated preserving directionality it doesn't make any\n",
      "sense to look at the data in a strand-specific manner). Strand-specific output files can\n",
      "optionally be skipped, in which case only three output files for CpG, CHG or CHH context\n",
      "will be generated. For both the strand-specific and comprehensive outputs there is also\n",
      "the option to merge both non-CpG contexts (CHG and CHH) into one single non-CpG context.\n",
      "\n",
      "\n",
      "The output files are in the following format (tab delimited):\n",
      "\n",
      "<sequence_id>     <strand>      <chromosome>     <position>     <methylation call>\n",
      "\n",
      "\n",
      "USAGE: bismark_methylation_extractor [options] <filenames>\n",
      "\n",
      "\n",
      "ARGUMENTS:\n",
      "==========\n",
      "\n",
      "<filenames>              A space-separated list of Bismark result files in SAM format from\n",
      "                         which methylation information is extracted for every cytosine in\n",
      "                         the reads. For alignment files in the older custom Bismark output\n",
      "                         see option '--vanilla'.\n",
      "\n",
      "OPTIONS:\n",
      "\n",
      "-s/--single-end          Input file(s) are Bismark result file(s) generated from single-end\n",
      "                         read data. If neither -s nor -p is set the type of experiment will\n",
      "                         be determined automatically.\n",
      "\n",
      "-p/--paired-end          Input file(s) are Bismark result file(s) generated from paired-end\n",
      "                         read data. If neither -s nor -p is set the type of experiment will\n",
      "                         be determined automatically.\n",
      "\n",
      "--vanilla                The Bismark result input file(s) are in the old custom Bismark format\n",
      "                         (up to version 0.5.x) and not in SAM format which is the default as\n",
      "                         of Bismark version 0.6.x or higher. Default: OFF.\n",
      "\n",
      "--no_overlap             For paired-end reads it is theoretically possible that read_1 and\n",
      "                         read_2 overlap. This option avoids scoring overlapping methylation\n",
      "                         calls twice (only methylation calls of read 1 are used for in the process\n",
      "                         since read 1 has historically higher quality basecalls than read 2).\n",
      "                         Whilst this option removes a bias towards more methylation calls\n",
      "                         in the center of sequenced fragments it may de facto remove a sizable\n",
      "                         proportion of the data. This option is on by default for paired-end data\n",
      "                         but can be disabled using --include_overlap. Default: ON.\n",
      "\n",
      "--include_overlap        For paired-end data all methylation calls will be extracted irrespective of\n",
      "                         of whether they overlap or not. Default: OFF.\n",
      "\n",
      "--ignore <int>           Ignore the first <int> bp from the 5' end of Read 1 (or single-end alignment\n",
      "                         files) when processing the methylation call string. This can remove e.g. a\n",
      "                         restriction enzyme site at the start of each read or any other source of\n",
      "                         bias (such as PBAT-Seq data).\n",
      "\n",
      "--ignore_r2 <int>        Ignore the first <int> bp from the 5' end of Read 2 of paired-end sequencing\n",
      "                         results only. Since the first couple of bases in Read 2 of BS-Seq experiments\n",
      "                         show a severe bias towards non-methylation as a result of end-repairing\n",
      "                         sonicated fragments with unmethylated cytosines (see M-bias plot), it is\n",
      "                         recommended that the first couple of bp of Read 2 are removed before\n",
      "                         starting downstream analysis. Please see the section on M-bias plots in the\n",
      "                         Bismark User Guide for more details.\n",
      "\n",
      "--ignore_3prime <int>    Ignore the last <int> bp from the 3' end of Read 1 (or single-end alignment\n",
      "                         files) when processing the methylation call string. This can remove unwanted\n",
      "                         biases from the end of reads.\n",
      "\n",
      "--ignore_3prime_r2 <int> Ignore the last <int> bp from the 3' end of Read 2 of paired-end sequencing\n",
      "                         results only. This can remove unwanted biases from the end of reads.\n",
      "\n",
      "--comprehensive          Specifying this option will merge all four possible strand-specific\n",
      "                         methylation info into context-dependent output files. The default\n",
      "\n",
      "                         contexts are:\n",
      "                          - CpG context\n",
      "                          - CHG context\n",
      "                          - CHH context\n",
      "\n",
      "--merge_non_CpG          This will produce two output files (in --comprehensive mode) or eight\n",
      "                         strand-specific output files (default) for Cs in\n",
      "                          - CpG context\n",
      "                          - non-CpG context\n",
      "\n",
      "--report                 Prints out a short methylation summary as well as the paramaters used to run\n",
      "                         this script. Default: ON.\n",
      "\n",
      "--no_header              Suppresses the Bismark version header line in all output files for more convenient\n",
      "                         batch processing.\n",
      "\n",
      "-o/--output DIR          Allows specification of a different output directory (absolute or relative\n",
      "                         path). If not specified explicitly, the output will be written to the current directory.\n",
      "\n",
      "--samtools_path          The path to your Samtools installation, e.g. /home/user/samtools/. Does not need to be specified\n",
      "                         explicitly if Samtools is in the PATH already.\n",
      "\n",
      "--gzip                   The methylation extractor files (CpG_OT_..., CpG_OB_... etc) will be written out in\n",
      "                         a GZIP compressed form to save disk space. This option is also passed on to the genome-wide\n",
      "                         cytosine report. BedGraph and coverage files will be written out as .gz by default.\n",
      "\n",
      "--version                Displays version information.\n",
      "\n",
      "-h/--help                Displays this help file and exits.\n",
      "\n",
      "--mbias_only             The methylation extractor will read the entire file but only output the M-bias table and plots as \n",
      "                         well as a report (optional) and then quit. Default: OFF.\n",
      "\n",
      "--mbias_off              The methylation extractor will process the entire file as usual but doesn't write out any M-bias report.\n",
      "                         Only recommended for users who deliberately want to keep an earlier version of the M-bias report. \n",
      "                         Default: OFF.\n",
      "\n",
      "--parallel <int>         May also be --multicore <int>. Sets the number of cores to be used for the methylation extraction process.\n",
      "                         If system resources are plentiful this is a viable option to speed up the extraction process (we observed a\n",
      "                         near linear speed increase for up to 10 cores used). Please note that a typical process of extracting a BAM file\n",
      "                         and writing out '.gz' output streams will in fact use ~3 cores per value of --parallel <int>\n",
      "                         specified (1 for the methylation extractor itself, 1 for a Samtools stream, 1 for GZIP stream), so\n",
      "                         --parallel 10 is likely to use around 30 cores of system resources. This option has no bearing\n",
      "                         on the bismark2bedGraph or genome-wide cytosine report processes.\n",
      "\n",
      "--yacht                  This option (Yet Another Context Hunting Tool) writes out additional information about the read a methylation\n",
      "                         call belongs to, and its output is meant to be fed into the NOMe_filtering script. This option writes out\n",
      "                         a single 'any_C_context' file that contains all methylation calls for a read consecutively. Its intended use\n",
      "                         is single-cell NOMe-Seq data, and thus this option works only in single-end mode (paired-end reads often suffer\n",
      "                         from chimaera problems...)\n",
      "\n",
      "                         --yacht will add three additional columns to the standard methylation call files:\n",
      "                                    <read start>  <read end>  <read orientation>\n",
      "                         For forward reads (+ orientation) the start position is the left-most position wheras for reverse reads\n",
      "                         (- orientation) it is the rightmost position.\n",
      "\n",
      "\n",
      "bedGraph specific options:\n",
      "==========================\n",
      "\n",
      "--bedGraph               After finishing the methylation extraction, the methylation output is written into a\n",
      "                         sorted bedGraph file that reports the position of a given cytosine and its methylation \n",
      "                         state (in %, see details below). The methylation extractor output is temporarily split up into\n",
      "                         temporary files, one per chromosome (written into the current directory or folder\n",
      "                         specified with -o/--output); these temp files are then used for sorting and deleted\n",
      "                         afterwards. By default, only cytosines in CpG context will be sorted. The option\n",
      "                         '--CX_context' may be used to report all cytosines irrespective of sequence context\n",
      "                         (this will take MUCH longer!). The default folder for temporary files during the sorting\n",
      "                         process is the output directory. The bedGraph conversion step is performed by the external\n",
      "                         module 'bismark2bedGraph'; this script needs to reside in the same folder as the \n",
      "                         bismark_methylation_extractor itself.\n",
      "\n",
      "--zero_based             Write out an additional coverage file (ending in .zero.cov) that uses 0-based genomic start\n",
      "                         and 1-based genomic end coordinates (zero-based, half-open), like used in the bedGraph file,\n",
      "                         instead of using 1-based coordinates throughout. Default: OFF.\n",
      "\n",
      "\n",
      "--cutoff [threshold]     The minimum number of times any methylation state (methylated or unmethylated) has to be seen\n",
      "                         for a nucleotide before its methylation percentage is reported. Default: 1.\n",
      "\n",
      "--remove_spaces          Replaces whitespaces in the sequence ID field with underscores to allow sorting.\n",
      "\n",
      "--CX/--CX_context        The sorted bedGraph output file contains information on every single cytosine that was covered\n",
      "                         in the experiment irrespective of its sequence context. This applies to both forward and\n",
      "                         reverse strands. Please be aware that this option may generate large temporary and output files\n",
      "                         and may take a long time to sort (up to many hours). Default: OFF.\n",
      "                         (i.e. Default = CpG context only).\n",
      "\n",
      "--buffer_size <string>   This allows you to specify the main memory sort buffer when sorting the methylation information.\n",
      "                         Either specify a percentage of physical memory by appending % (e.g. --buffer_size 50%) or\n",
      "\t\t\t a multiple of 1024 bytes, e.g. 'K' multiplies by 1024, 'M' by 1048576 and so on for 'T' etc. \n",
      "                         (e.g. --buffer_size 20G). For more information on sort type 'info sort' on a command line.\n",
      "                         Defaults to 2G.\n",
      "\n",
      "--scaffolds/--gazillion  Users working with unfinished genomes sporting tens or even hundreds of thousands of\n",
      "                         scaffolds/contigs/chromosomes frequently encountered errors with pre-sorting reads to\n",
      "                         individual chromosome files. These errors were caused by the operating system's limit\n",
      "                         of the number of filehandle that can be written to at any one time (typically 1024; to\n",
      "                         find out this limit on Linux, type: ulimit -a).\n",
      "                         To bypass the limitation of open filehandles, the option --scaffolds does not pre-sort\n",
      "                         methylation calls into individual chromosome files. Instead, all input files are\n",
      "                         temporarily merged into a single file (unless there is only a single file), and this\n",
      "                         file will then be sorted by both chromosome AND position using the Unix sort command.\n",
      "                         Please be aware that this option might take a looooong time to complete, depending on\n",
      "                         the size of the input files, and the memory you allocate to this process (see --buffer_size).\n",
      "                         Nevertheless, it seems to be working.\n",
      "\n",
      "--ample_memory           Using this option will not sort chromosomal positions using the UNIX 'sort' command, but will\n",
      "                         instead use two arrays to sort methylated and unmethylated calls. This may result in a faster\n",
      "                         sorting process of very large files, but this comes at the cost of a larger memory footprint\n",
      "                         (two arrays of the length of the largest human chromosome 1 (~250M bp) consume around 16GB\n",
      "                         of RAM). Due to overheads in creating and looping through these arrays it seems that it will\n",
      "                         actually be *slower* for small files (few million alignments), and we are currently testing at\n",
      "                         which point it is advisable to use this option. Note that --ample_memory is not compatible\n",
      "                         with options '--scaffolds/--gazillion' (as it requires pre-sorted files to begin with).\n",
      "\n",
      "\n",
      "\n",
      "Genome-wide cytosine methylation report specific options:\n",
      "=========================================================\n",
      "\n",
      "--cytosine_report        After the conversion to bedGraph has completed, the option '--cytosine_report' produces a\n",
      "                         genome-wide methylation report for all cytosines in the genome. By default, the output uses 1-based\n",
      "                         chromosome coordinates (zero-based start coords are optional) and reports CpG context only (all\n",
      "                         cytosine context is optional). The output considers all Cs on both forward and reverse strands and\n",
      "                         reports their position, strand, trinucleotide content and methylation state (counts are 0 if not\n",
      "                         covered). The cytosine report conversion step is performed by the external module\n",
      "                         'coverage2cytosine'; this script needs to reside in the same folder as the bismark_methylation_extractor\n",
      "                         itself.\n",
      "\n",
      "--CX/--CX_context        The output file contains information on every single cytosine in the genome irrespective of\n",
      "                         its context. This applies to both forward and reverse strands. Please be aware that this will\n",
      "                         generate output files with > 1.1 billion lines for a mammalian genome such as human or mouse.\n",
      "                         Default: OFF (i.e. Default = CpG context only).\n",
      "\n",
      "--zero_based             Uses 0-based genomic coordinates instead of 1-based coordinates. Default: OFF.\n",
      "\n",
      "--genome_folder <path>   Enter the genome folder you wish to use to extract sequences from (full path only). Accepted\n",
      "                         formats are FastA files ending with '.fa' or '.fasta'. Specifying a genome folder path is mandatory.\n",
      "\n",
      "--split_by_chromosome    Writes the output into individual files for each chromosome instead of a single output file. Files\n",
      "                         will be named to include the input filename and the chromosome number.\n",
      "\n",
      "\n",
      "\n",
      "OUTPUT:\n",
      "\n",
      "The bismark_methylation_extractor output is in the form:\n",
      "========================================================\n",
      "<seq-ID>  <methylation state*>  <chromosome>  <start position (= end position)>  <methylation call>\n",
      "\n",
      "* Methylated cytosines receive a '+' orientation,\n",
      "* Unmethylated cytosines receive a '-' orientation.\n",
      "\n",
      "The bismark_methylation_extractor output with --yacht (optional) specified is in the form:\n",
      "==========================================================================================\n",
      "<seq-ID>  <methylation state*>  <chromosome>  <start position (= end position)>  <methylation call>  <read start>  <read end>  <read orientation>\n",
      "\n",
      "* Methylated cytosines receive a '+' orientation,\n",
      "* Unmethylated cytosines receive a '-' orientation.\n",
      "\n",
      "\n",
      "The bedGraph output (optional) looks like this (tab-delimited; 0-based start coords, 1-based end coords):\n",
      "=========================================================================================================\n",
      "\n",
      "track type=bedGraph (header line)\n",
      "\n",
      "<chromosome>  <start position>  <end position>  <methylation percentage>\n",
      "\n",
      "\n",
      "The coverage output looks like this (tab-delimited, 1-based genomic coords; zero-based half-open coordinates available with '--zero_based'):\n",
      "============================================================================================================================================\n",
      "\n",
      "<chromosome>  <start position>  <end position>  <methylation percentage>  <count methylated>  <count non-methylated>\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "The genome-wide cytosine methylation output file is tab-delimited in the following format:\n",
      "==========================================================================================\n",
      "<chromosome>  <position>  <strand>  <count methylated>  <count non-methylated>  <C-context>  <trinucleotide context>\n",
      "\n",
      "\n",
      "\n",
      "This script was last modified on 05 October 2017.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "! ../../../../../Shared/Apps/Bismark_v0.19.0/bismark_methylation_extractor -help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "1. Path to `bismark_methylation_extractor`\n",
    "2. -p: Indicate I have paired reads\n",
    "3. --bedGraph: Create `bedGraph` file with results\n",
    "4. --counts\n",
    "5. --scaffolds: Bypass limitation of open filehandles\n",
    "6. --multicore 8\n",
    "7. --gzip: Gunzip files to save space\n",
    "8. Path to .bam files generated in Step 2\n",
    "9. Path for standard error output file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Now testing Bismark result file zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam for positional sorting (which would be bad...)\tzr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8\n",
      "\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1.mbias\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2.mbias\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3.mbias\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4.mbias\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5.mbias\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6.mbias\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7.mbias\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8.mbias\n",
      "\n",
      "/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OT_zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOT_zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOB_zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OB_zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\n",
      "\n",
      "\n",
      "Now testing Bismark result file zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam for positional sorting (which would be bad...)\tzr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8\n",
      "\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1.mbias\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2.mbias\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3.mbias\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4.mbias\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5.mbias\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6.mbias\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7.mbias\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8.mbias\n",
      "\n",
      "/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OT_zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOT_zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOB_zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OB_zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\n",
      "\n",
      "\n",
      "Now testing Bismark result file zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam for positional sorting (which would be bad...)\tzr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8\n",
      "\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1.mbias\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2.mbias\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3.mbias\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4.mbias\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5.mbias\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6.mbias\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7.mbias\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8.mbias\n",
      "\n",
      "/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OT_zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOT_zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOB_zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OB_zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\n",
      "\n",
      "\n",
      "Now testing Bismark result file zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam for positional sorting (which would be bad...)\tzr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8\n",
      "\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1.mbias\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2.mbias\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3.mbias\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4.mbias\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5.mbias\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6.mbias\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7.mbias\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8.mbias\n",
      "\n",
      "/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OT_zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOT_zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOB_zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OB_zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\n",
      "\n",
      "\n",
      "Now testing Bismark result file zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam for positional sorting (which would be bad...)\tzr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8\n",
      "\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1.mbias\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2.mbias\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3.mbias\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4.mbias\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5.mbias\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6.mbias\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7.mbias\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8.mbias\n",
      "\n",
      "/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OT_zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOT_zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOB_zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OB_zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\n",
      "\n",
      "\n",
      "Now testing Bismark result file zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam for positional sorting (which would be bad...)\tzr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8\n",
      "\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1.mbias\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2.mbias\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3.mbias\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4.mbias\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5.mbias\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6.mbias\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7.mbias\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8.mbias\n",
      "\n",
      "/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OT_zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOT_zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOB_zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OB_zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\n",
      "\n",
      "\n",
      "Now testing Bismark result file zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam for positional sorting (which would be bad...)\tzr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8\n",
      "\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1.mbias\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2.mbias\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3.mbias\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4.mbias\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5.mbias\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6.mbias\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7.mbias\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8.mbias\n",
      "\n",
      "/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OT_zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOT_zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOB_zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OB_zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\n",
      "\n",
      "\n",
      "Now testing Bismark result file zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam for positional sorting (which would be bad...)\tzr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8\n",
      "\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1.mbias\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2.mbias\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3.mbias\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4.mbias\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5.mbias\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6.mbias\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7.mbias\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8.mbias\n",
      "\n",
      "/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OT_zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOT_zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOB_zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OB_zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\n",
      "\n",
      "\n",
      "Now testing Bismark result file zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam for positional sorting (which would be bad...)\tzr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8\n",
      "\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1.mbias\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2.mbias\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3.mbias\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4.mbias\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5.mbias\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6.mbias\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7.mbias\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8.mbias\n",
      "\n",
      "/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OT_zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOT_zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOB_zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OB_zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\n",
      "\n",
      "\n",
      "Now testing Bismark result file zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam for positional sorting (which would be bad...)\tzr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8\n",
      "\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.1.mbias\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.2.mbias\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.3.mbias\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.4.mbias\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.5.mbias\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.6.mbias\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.7.mbias\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt.8.mbias\n",
      "\n",
      "/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OT_zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOT_zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_CTOB_zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\t/Users/yaamini/Documents/project-virginica-oa/analyses/2018-05-22-Bismark-Full-Samples/CpG_OB_zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated.txt\n",
      "\n"
     ]
    }
   ],
   "source": [
    "! ../../../../../Shared/Apps/Bismark_v0.19.0/bismark_methylation_extractor \\\n",
    "-p \\\n",
    "--bedGraph \\\n",
    "--counts \\\n",
    "--scaffolds \\\n",
    "--multicore 8 \\\n",
    "*deduplicated.bam \\\n",
    "2> bme-full-0524.err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\r\n",
      "zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz.methylation_calls.merged\r\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\r\n",
      "zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz.methylation_calls.merged\r\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\r\n",
      "zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz.methylation_calls.merged\r\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\r\n",
      "zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz.methylation_calls.merged\r\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\r\n",
      "zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz.methylation_calls.merged\r\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\r\n",
      "zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz.methylation_calls.merged\r\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\r\n",
      "zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz.methylation_calls.merged\r\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\r\n",
      "zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz.methylation_calls.merged\r\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\r\n",
      "zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz.methylation_calls.merged\r\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\r\n",
      "zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz.methylation_calls.merged\r\n"
     ]
    }
   ],
   "source": [
    "!ls *bedGraph*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 5. HTML Processing Report"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n",
      "  SYNOPSIS:\r\n",
      "\r\n",
      "  This script uses a Bismark alignment report to generate a graphical HTML report page. Optionally, further reports of\r\n",
      "  the Bismark suite such as deduplication, methylation extractor splitting or M-bias reports can be specified as well. If several\r\n",
      "  Bismark reports are found in the same folder, a separate report will be generated for each of these, whereby the output filename\r\n",
      "  will be derived from the Bismark alignment report file. Bismark2report attempts to find optional reports automatically based\r\n",
      "  on the file basename.\r\n",
      "\r\n",
      "\r\n",
      "  USAGE: bismark2report [options]\r\n",
      "\r\n",
      "\r\n",
      "-o/--output <filename>     Name of the output file (optional). If not specified explicitly, the output filename will be derived\r\n",
      "                           from the Bismark alignment report file. Specifying an output filename only works if the HTML report is\r\n",
      "                           to be generated for a single Bismark alignment report (and potentially additional reports).\r\n",
      "\r\n",
      "--dir                      Output directory. Output is written to the current directory if not specified explicitly.\r\n",
      "\r\n",
      "\r\n",
      "--alignment_report FILE    If not specified explicitly, bismark2report attempts to find Bismark report file(s) in the current\r\n",
      "                           directory and produces a separate HTML report for each mapping report file. Based on the basename of\r\n",
      "                           the Bismark mapping report, bismark2report will also attempt to find the other Bismark reports (see below)\r\n",
      "                           for inclusion into the HTML report. Specifying a Bismark alignment report file is mandatory.\r\n",
      "\r\n",
      "--dedup_report FILE        If not specified explicitly, bismark2report attempts to find a deduplication report file with the same\r\n",
      "                           basename as the Bismark mapping report (generated by deduplicate_bismark) in the\r\n",
      "                           current working directory. Including a deduplication report is optional, and using the FILE 'none'\r\n",
      "                           will skip this step entirely.\r\n",
      "\r\n",
      "--splitting_report FILE    If not specified explicitly, bismark2report attempts to find a splitting report file with the same\r\n",
      "                           basename as the Bismark mapping report (generated by the Bismark methylation extractor) in the current\r\n",
      "                           working directory. Including a splitting report is optional, and using the FILE 'none' will skip this\r\n",
      "                           step entirely.\r\n",
      "\r\n",
      "--mbias_report FILE        If not specified explicitly, bismark2report attempts to find a single M-bias report file with the same\r\n",
      "                           basename as the Bismark mapping report (generated by the Bismark methylation extractor) in the current\r\n",
      "                           working directory. Including an M-Bias report is optional, and using the FILE 'none' will skip this step\r\n",
      "                           entirely.\r\n",
      "\r\n",
      "--nucleotide_report FILE   If not specified explicitly, bismark2report attempts to find a single nucleotide coverage report file\r\n",
      "                           with the same basename as the Bismark mapping report (generated by Bismark with the option\r\n",
      "                           '--nucleotide_coverage') in the current working directory. Including a nucleotide coverage statistics\r\n",
      "                           report is optional, and using the FILE 'none' will skip this report entirely.\r\n",
      "\r\n",
      "                           Script last modified: 13 May 2016\r\n",
      "\r\n"
     ]
    }
   ],
   "source": [
    "! ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2report -help"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 10 alignment reports in current directory. Now trying to figure out whether there are corresponding optional reports\n",
      "\n",
      "Writing Bismark HTML report to >> zr2096_10_s1_R1_val_1_bismark_bt2_PE_report.html <<\n",
      "\n",
      "Redundant argument in sprintf at ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2report line 130.\n",
      "==============================================================================================================\n",
      "Using the following alignment report:\t\t> zr2096_10_s1_R1_val_1_bismark_bt2_PE_report.txt <\n",
      "Processing alignment report zr2096_10_s1_R1_val_1_bismark_bt2_PE_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following deduplication report:\t> zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt <\n",
      "Processing deduplication report zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following splitting report:\t\t> zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt <\n",
      "Processing splitting report zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following M-bias report:\t\t> zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt <\n",
      "Processing M-bias report zr2096_10_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt ...\n",
      "Complete\n",
      "\n",
      "No nucleotide coverage report file specified, skipping this step\n",
      "==============================================================================================================\n",
      "\n",
      "\n",
      "\n",
      "Writing Bismark HTML report to >> zr2096_1_s1_R1_val_1_bismark_bt2_PE_report.html <<\n",
      "\n",
      "Redundant argument in sprintf at ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2report line 130.\n",
      "==============================================================================================================\n",
      "Using the following alignment report:\t\t> zr2096_1_s1_R1_val_1_bismark_bt2_PE_report.txt <\n",
      "Processing alignment report zr2096_1_s1_R1_val_1_bismark_bt2_PE_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following deduplication report:\t> zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt <\n",
      "Processing deduplication report zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following splitting report:\t\t> zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt <\n",
      "Processing splitting report zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following M-bias report:\t\t> zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt <\n",
      "Processing M-bias report zr2096_1_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt ...\n",
      "Complete\n",
      "\n",
      "No nucleotide coverage report file specified, skipping this step\n",
      "==============================================================================================================\n",
      "\n",
      "\n",
      "\n",
      "Writing Bismark HTML report to >> zr2096_2_s1_R1_val_1_bismark_bt2_PE_report.html <<\n",
      "\n",
      "Redundant argument in sprintf at ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2report line 130.\n",
      "==============================================================================================================\n",
      "Using the following alignment report:\t\t> zr2096_2_s1_R1_val_1_bismark_bt2_PE_report.txt <\n",
      "Processing alignment report zr2096_2_s1_R1_val_1_bismark_bt2_PE_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following deduplication report:\t> zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt <\n",
      "Processing deduplication report zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following splitting report:\t\t> zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt <\n",
      "Processing splitting report zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following M-bias report:\t\t> zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt <\n",
      "Processing M-bias report zr2096_2_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt ...\n",
      "Complete\n",
      "\n",
      "No nucleotide coverage report file specified, skipping this step\n",
      "==============================================================================================================\n",
      "\n",
      "\n",
      "\n",
      "Writing Bismark HTML report to >> zr2096_3_s1_R1_val_1_bismark_bt2_PE_report.html <<\n",
      "\n",
      "Redundant argument in sprintf at ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2report line 130.\n",
      "==============================================================================================================\n",
      "Using the following alignment report:\t\t> zr2096_3_s1_R1_val_1_bismark_bt2_PE_report.txt <\n",
      "Processing alignment report zr2096_3_s1_R1_val_1_bismark_bt2_PE_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following deduplication report:\t> zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt <\n",
      "Processing deduplication report zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following splitting report:\t\t> zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt <\n",
      "Processing splitting report zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following M-bias report:\t\t> zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt <\n",
      "Processing M-bias report zr2096_3_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt ...\n",
      "Complete\n",
      "\n",
      "No nucleotide coverage report file specified, skipping this step\n",
      "==============================================================================================================\n",
      "\n",
      "\n",
      "\n",
      "Writing Bismark HTML report to >> zr2096_4_s1_R1_val_1_bismark_bt2_PE_report.html <<\n",
      "\n",
      "Redundant argument in sprintf at ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2report line 130.\n",
      "==============================================================================================================\n",
      "Using the following alignment report:\t\t> zr2096_4_s1_R1_val_1_bismark_bt2_PE_report.txt <\n",
      "Processing alignment report zr2096_4_s1_R1_val_1_bismark_bt2_PE_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following deduplication report:\t> zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt <\n",
      "Processing deduplication report zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following splitting report:\t\t> zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt <\n",
      "Processing splitting report zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following M-bias report:\t\t> zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt <\n",
      "Processing M-bias report zr2096_4_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt ...\n",
      "Complete\n",
      "\n",
      "No nucleotide coverage report file specified, skipping this step\n",
      "==============================================================================================================\n",
      "\n",
      "\n",
      "\n",
      "Writing Bismark HTML report to >> zr2096_5_s1_R1_val_1_bismark_bt2_PE_report.html <<\n",
      "\n",
      "Redundant argument in sprintf at ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2report line 130.\n",
      "==============================================================================================================\n",
      "Using the following alignment report:\t\t> zr2096_5_s1_R1_val_1_bismark_bt2_PE_report.txt <\n",
      "Processing alignment report zr2096_5_s1_R1_val_1_bismark_bt2_PE_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following deduplication report:\t> zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt <\n",
      "Processing deduplication report zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following splitting report:\t\t> zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt <\n",
      "Processing splitting report zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following M-bias report:\t\t> zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt <\n",
      "Processing M-bias report zr2096_5_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt ...\n",
      "Complete\n",
      "\n",
      "No nucleotide coverage report file specified, skipping this step\n",
      "==============================================================================================================\n",
      "\n",
      "\n",
      "\n",
      "Writing Bismark HTML report to >> zr2096_6_s1_R1_val_1_bismark_bt2_PE_report.html <<\n",
      "\n",
      "Redundant argument in sprintf at ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2report line 130.\n",
      "==============================================================================================================\n",
      "Using the following alignment report:\t\t> zr2096_6_s1_R1_val_1_bismark_bt2_PE_report.txt <\n",
      "Processing alignment report zr2096_6_s1_R1_val_1_bismark_bt2_PE_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following deduplication report:\t> zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt <\n",
      "Processing deduplication report zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following splitting report:\t\t> zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt <\n",
      "Processing splitting report zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following M-bias report:\t\t> zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt <\n",
      "Processing M-bias report zr2096_6_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt ...\n",
      "Complete\n",
      "\n",
      "No nucleotide coverage report file specified, skipping this step\n",
      "==============================================================================================================\n",
      "\n",
      "\n",
      "\n",
      "Writing Bismark HTML report to >> zr2096_7_s1_R1_val_1_bismark_bt2_PE_report.html <<\n",
      "\n",
      "Redundant argument in sprintf at ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2report line 130.\n",
      "==============================================================================================================\n",
      "Using the following alignment report:\t\t> zr2096_7_s1_R1_val_1_bismark_bt2_PE_report.txt <\n",
      "Processing alignment report zr2096_7_s1_R1_val_1_bismark_bt2_PE_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following deduplication report:\t> zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt <\n",
      "Processing deduplication report zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following splitting report:\t\t> zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt <\n",
      "Processing splitting report zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following M-bias report:\t\t> zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt <\n",
      "Processing M-bias report zr2096_7_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt ...\n",
      "Complete\n",
      "\n",
      "No nucleotide coverage report file specified, skipping this step\n",
      "==============================================================================================================\n",
      "\n",
      "\n",
      "\n",
      "Writing Bismark HTML report to >> zr2096_8_s1_R1_val_1_bismark_bt2_PE_report.html <<\n",
      "\n",
      "Redundant argument in sprintf at ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2report line 130.\n",
      "==============================================================================================================\n",
      "Using the following alignment report:\t\t> zr2096_8_s1_R1_val_1_bismark_bt2_PE_report.txt <\n",
      "Processing alignment report zr2096_8_s1_R1_val_1_bismark_bt2_PE_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following deduplication report:\t> zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt <\n",
      "Processing deduplication report zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following splitting report:\t\t> zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt <\n",
      "Processing splitting report zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following M-bias report:\t\t> zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt <\n",
      "Processing M-bias report zr2096_8_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt ...\n",
      "Complete\n",
      "\n",
      "No nucleotide coverage report file specified, skipping this step\n",
      "==============================================================================================================\n",
      "\n",
      "\n",
      "\n",
      "Writing Bismark HTML report to >> zr2096_9_s1_R1_val_1_bismark_bt2_PE_report.html <<\n",
      "\n",
      "Redundant argument in sprintf at ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2report line 130.\n",
      "==============================================================================================================\n",
      "Using the following alignment report:\t\t> zr2096_9_s1_R1_val_1_bismark_bt2_PE_report.txt <\n",
      "Processing alignment report zr2096_9_s1_R1_val_1_bismark_bt2_PE_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following deduplication report:\t> zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt <\n",
      "Processing deduplication report zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplication_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following splitting report:\t\t> zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt <\n",
      "Processing splitting report zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated_splitting_report.txt ...\n",
      "Complete\n",
      "\n",
      "Using the following M-bias report:\t\t> zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt <\n",
      "Processing M-bias report zr2096_9_s1_R1_val_1_bismark_bt2_pe.deduplicated.M-bias.txt ...\n",
      "Complete\n",
      "\n",
      "No nucleotide coverage report file specified, skipping this step\n",
      "==============================================================================================================\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "! ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2report"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Summary Report"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n",
      "  SYNOPSIS:\r\n",
      "\r\n",
      "  This script uses Bismark report files of several samples in a run folder to generate a graphical summary HTML report as well as\r\n",
      "  a whopping big table (tab delimited text) with all relevant alignment and methylation statistics which may be used for graphing\r\n",
      "  purposes in R, Excel or the like. Unless specific BAM files are specified, bismark2summary first identifies Bismark BAM files in\r\n",
      "  a folder (they need to use the Bismark naming conventions) and then automatically detects Bismark alignment, deduplication or\r\n",
      "  methylation extractor (splitting) reports based on the input file basename. If splitting reports are found they overwrite the\r\n",
      "  methylation statistics of the initial alignment report.\r\n",
      "\r\n",
      "\r\n",
      "  USAGE: bismark2summary [options] [<BAM file(s)>]\r\n",
      "\r\n",
      "  ARGUMENTS:\r\n",
      "\r\n",
      "  BAM file(s)                 Optional. If no BAM files are specified explicitly the current working directory is scanned for \r\n",
      "                              Bismark alignment files and their associated reports. \r\n",
      "\r\n",
      "  OPTIONS:\r\n",
      "\r\n",
      "  -o/--basename <filename>    Basename of the output file (optional). Generate a text file with all relevant extracted values\r\n",
      "                              ('basename.txt') as well as an HTML report ('basename.html'). If not specified explicitly, the\r\n",
      "                              the basename is 'bismark_summary_report'.\r\n",
      "\r\n",
      "  --title <string>            Optional HTML report title; use --title \"speech marks for text with spaces\". Default: \r\n",
      "                              'Bismark Summary Report'. \t\t\t      \r\n",
      "\r\n",
      "\r\n",
      "  --version                   Displays version information and exits.\r\n",
      "\t\t\t      \r\n",
      "  --help                      Displays this help message and exits.\r\n",
      "\r\n",
      "\r\n",
      "                                                Script last modified: 06 July 2016\r\n",
      "\r\n"
     ]
    }
   ],
   "source": [
    "! ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2summary -help"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No Bismark/Bowtie2 single-end BAM files detected\r\n",
      "Found Bismark/Bowtie2 paired-end files\r\n",
      "No Bismark/Bowtie single-end BAM files detected\r\n",
      "No Bismark/Bowtie paired-end BAM files detected\r\n",
      "\r\n",
      "Generating Bismark summary report from 10 Bismark BAM file(s)...\r\n",
      ">> Reading from Bismark report: zr2096_10_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      ">> Reading from Bismark report: zr2096_1_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      ">> Reading from Bismark report: zr2096_2_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      ">> Reading from Bismark report: zr2096_3_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      ">> Reading from Bismark report: zr2096_4_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      ">> Reading from Bismark report: zr2096_5_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      ">> Reading from Bismark report: zr2096_6_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      ">> Reading from Bismark report: zr2096_7_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      ">> Reading from Bismark report: zr2096_8_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      ">> Reading from Bismark report: zr2096_9_s1_R1_val_1_bismark_bt2_PE_report.txt\r\n",
      "\r\n",
      "Wrote Bismark project summary to >> bismark_summary_report.html <<\r\n",
      "\r\n"
     ]
    }
   ],
   "source": [
    "! ../../../../../Shared/Apps/Bismark_v0.19.0/bismark2summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}