Answers.Body: <p>There are <a href="https://hicexplorer.readthedocs.io/en/latest/content/example_usage.html" rel="nofollow noreferrer">some interesting suggestions from the HiCExplorer developers</a>:</p>
  
  <blockquote>
    <p>Mates have to be mapped individually to avoid mapper specific
    heuristics designed for standard paired-end libraries.</p>
  </blockquote>
  
  <pre><code># map the reads, each mate individually using
  # for example bwa
  #
  # bwa mem mapping options:
  #       -A INT        score for a sequence match, which scales options -TdBOELU unless overridden [1]
  #       -B INT        penalty for a mismatch [4]
  #       -O INT[,INT]  gap open penalties for deletions and insertions [6,6]
  #       -E INT[,INT]  gap extension penalty; a gap of size k cost '{-O} + {-E}*k' [1,1] # this is set very high to avoid gaps
  #                                  at restriction sites. Setting the gap extension penalty high, produces better results as
  #                                  the sequences left and right of a restriction site are mapped independently.
  #       -L INT[,INT]  penalty for 5'- and 3'-end clipping [5,5] # this is set to no penalty.
  
  $ bwa mem -A1 -B4  -E50 -L0  index_path \
      -U mate_R1.fastq.gz 2&gt;&gt;mate_R1.log | samtools view -Shb - &gt; mate_R1.bam
  
  $ bwa mem -A1 -B4  -E50 -L0  index_path \
      -U mate_R2.fastq.gz 2&gt;&gt;mate_R2.log | samtools view -Shb - &gt; mate_R2.bam
  
  # build matrix from independently mated read pairs
  # the restriction sequence GATC is recognized by the DpnII restriction enzyme
  
  $ hicBuildMatrix --samFiles mate_R1.bam mate_R2.bam \
                   --binSize 10000 \
                   --restrictionSequence GATC \
                   --threads 4
                   --inputBufferSize 100000
                   --outBam hic.bam \
                   -o hic_matrix.h5
                   --QCfolder ./hicQC
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-09-27T01:57:00.487Z
Answers.Text: Option -SP exactly aims to disable "mapper specific heuristics designed for standard paired-end libraries".
Answers.UserDisplayName: user172818
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-09-26T01:17:45.883Z
Answers.LastActivityDate: 2018-09-26T01:17:45.883Z
Answers.UserDisplayName: burger
Answers.Body: <p>See <a href="https://github.com/hms-dbmi/hic-data-analysis-bootcamp" rel="nofollow noreferrer">here</a>, in particular <a href="https://hms-dbmi.github.io/hic-data-analysis-bootcamp/#10" rel="nofollow noreferrer">slide #10</a> from the tutorial:</p>
  
  <pre class="lang-bash prettyprint-override"><code>bwa mem -5SPM ref.fa read1.fq read2.fq &gt; out.sam
  </code></pre>
  
  <p>Here, <code>-SP</code> disables pairing. The <a href="http://aidenlab.org" rel="nofollow noreferrer">Aiden lab</a> and <a href="https://www.ncbi.nlm.nih.gov/pubmed/30166492" rel="nofollow noreferrer">this paper</a> also use a similar command line. Beware that there are many Hi-C pipelines, but not many are using bwa-mem.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-03T23:18:14.290Z
Answers.Text: Thanks for your post. It looks like the -SP flags are sufficient to map the trimmed data as I discussed above, while -5SPM is appropriate for Hi-C data that do not have the linker sequence trimmed.
Answers.UserDisplayName: conchoecia
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-09-23T04:37:45.580Z
Answers.LastActivityDate: 2018-09-23T04:37:45.580Z
Answers.UserDisplayName: user172818
Body: <h2>Library Prep</h2>
  <ul>
  <li>I have a <a href="https://www.ncbi.nlm.nih.gov/pubmed/22652625" rel="noreferrer">Hi-C</a> library prepped using an enzyme that cuts at <code>GATC</code>, so it leaves <code>GATCGATC</code> as the junction sequences. This library was sequenced on a 2x150 PE Illumina run.</li>
  </ul>
  <h2>Data Pre-processing</h2>
  <ul>
  <li>The reads were adapter and quality trimming with <a href="http://www.usadellab.org/cms/?page=trimmomatic" rel="noreferrer">trimmomatic</a> paired-end mode, keeping reads that overlap.</li>
  <li>Finally, the forward and reverse reads were trimmed such that all bases 3' of the junction sequence (<code>GATCGATC</code>). For example, the read:
  <ul>
  <li><code>5' AAAATT**GATCGATC**AATTATGGATAA 3'</code> becomes: <code>5' AAAATT**GATC 3'</code></li>
  </ul>
  </li>
  <li>In this process, sequences that don't contain the junction sequences were left alone. Trimming junction sequences in this way makes it such that each read pair maps to only two places in the genome. No three-way or four-way chimeric mappings are possible since everything <code>3'</code> of the junction sequence was removed.</li>
  </ul>
  <h2>My Question</h2>
  <ul>
  <li>What is the correct way to now map these trimmed Hi-C reads to a genome using <code>bwa mem</code>?</li>
  </ul>
  <h2>Progress so far/things I have tried</h2>
  <ul>
  <li>In the past I have just tried mapping using <code>bwa mem</code> with no special parameters using this command: <code>bwa mem ref.fa hi-c_f.fq.gz hi-c_r.fq.gz</code>. One concern is that the computations that <code>bwa mem</code> performs on insert sizes in this mode will cause some mappings to be faulty or that read pairs with chromosome-length information will be removed from the alignment because they fall on the tail of the insert size distribution (non-Gaussian, but log decay).</li>
  <li>Another option to consider is mapping each file, <code>hi-c_f.fq.gz</code> and <code>hi-c_r.fq.gz</code>, independently then adding mate information with <code>samtools fixmate</code>. However, I do not know how to do this.</li>
  <li>A last option that I have considered is to use <code>bwa mem</code> in paired-end mode, but to use a special flag. In the past I believe that <code>-5</code> was the correct flag, but in <a href="https://sourceforge.net/p/bio-bwa/mailman/message/36050442/" rel="noreferrer">this post</a> Heng Li discusses that the <code>-5</code> flag behavior has changed, as well as a new <code>-q</code> flag. In <a href="http://lh3.github.io/2018/04/02/minimap2-and-the-future-of-bwa" rel="noreferrer">this post</a> Heng Li discusses that he thinks <code>bwa mem</code> is still a good option for mapping Hi-C data but does not discuss the proper flags.</li>
  </ul>
  <h2>Limitations</h2>
  <ul>
  <li>The answer should provide the correct way to map Hi-C using <code>bwa mem</code>, but other software suggestions will also be considered.</li>
  </ul>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-09-22T19:39:04.807Z
Comments.Text: This ends up depending on the software you want to use downstream of mapping. For example, HiCExplorer still needs reads to be aligned as single-end. I imagine it's not alone in that regard.
Comments.UserDisplayName: Devon Ryan
Tags: <read-mapping><bwa><illumina><hi-c>
Title: What is the correct way to map Hi-C data with bwa mem?
ThreadText: # What is the correct way to map Hi-C data with bwa mem?
  
  (Asked by: username_0 on 2018-09-21)
  
  ## Library Prep
  
  * I have a [Hi-C](https://www.ncbi.nlm.nih.gov/pubmed/22652625) library prepped using an enzyme that cuts at `GATC`, so it leaves `GATCGATC` as the junction sequences. This library was sequenced on a 2x150 PE Illumina run.
  
  ## Data Pre-processing
  
  * The reads were adapter and quality trimming with [trimmomatic](http://www.usadellab.org/cms/?page=trimmomatic) paired-end mode, keeping reads that overlap.
  * Finally, the forward and reverse reads were trimmed such that all bases 3' of the junction sequence (`GATCGATC`). For example, the read:
    + `5' AAAATT**GATCGATC**AATTATGGATAA 3'` becomes: `5' AAAATT**GATC 3'`
  * In this process, sequences that don't contain the junction sequences were left alone. Trimming junction sequences in this way makes it such that each read pair maps to only two places in the genome. No three-way or four-way chimeric mappings are possible since everything `3'` of the junction sequence was removed.
  
  ## My Question
  
  * What is the correct way to now map these trimmed Hi-C reads to a genome using `bwa mem`?
  
  ## Progress so far/things I have tried
  
  * In the past I have just tried mapping using `bwa mem` with no special parameters using this command: `bwa mem ref.fa hi-c_f.fq.gz hi-c_r.fq.gz`. One concern is that the computations that `bwa mem` performs on insert sizes in this mode will cause some mappings to be faulty or that read pairs with chromosome-length information will be removed from the alignment because they fall on the tail of the insert size distribution (non-Gaussian, but log decay).
  * Another option to consider is mapping each file, `hi-c_f.fq.gz` and `hi-c_r.fq.gz`, independently then adding mate information with `samtools fixmate`. However, I do not know how to do this.
  * A last option that I have considered is to use `bwa mem` in paired-end mode, but to use a special flag. In the past I believe that `-5` was the correct flag, but in [this post](https://sourceforge.net/p/bio-bwa/mailman/message/36050442/) Heng Li discusses that the `-5` flag behavior has changed, as well as a new `-q` flag. In [this post](http://lh3.github.io/2018/04/02/minimap2-and-the-future-of-bwa) Heng Li discusses that he thinks `bwa mem` is still a good option for mapping Hi-C data but does not discuss the proper flags.
  
  ## Limitations
  
  * The answer should provide the correct way to map Hi-C using `bwa mem`, but other software suggestions will also be considered.
  
  **Tags:** read-mapping, bwa, illumina, hi-c
  
  > **Comments:**
  > **username_1**: This ends up depending on the software you want to use downstream of mapping. For example, HiCExplorer still needs reads to be aligned as single-end. I imagine it's not alone in that regard.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_3 on 2018-09-23. Score: 4 | Accepted answer
  
  See [here](https://github.com/hms-dbmi/hic-data-analysis-bootcamp), in particular [slide #10](https://hms-dbmi.github.io/hic-data-analysis-bootcamp/#10) from the tutorial:
  
  ```bash
  bwa mem -5SPM ref.fa read1.fq read2.fq > out.sam
  
  ```
  
  Here, `-SP` disables pairing. The [Aiden lab](http://aidenlab.org) and [this paper](https://www.ncbi.nlm.nih.gov/pubmed/30166492) also use a similar command line. Beware that there are many Hi-C pipelines, but not many are using bwa-mem.
  
  > **Comments:**
  > **username_0**: Thanks for your post. It looks like the -SP flags are sufficient to map the trimmed data as I discussed above, while -5SPM is appropriate for Hi-C data that do not have the linker sequence trimmed.  
  
  ---
  
  ### Answer by username_2 on 2018-09-26. Score: 1
  
  There are [some interesting suggestions from the HiCExplorer developers](https://hicexplorer.readthedocs.io/en/latest/content/example_usage.html):
  
  > Mates have to be mapped individually to avoid mapper specific
  > heuristics designed for standard paired-end libraries.
  
  ```
  # map the reads, each mate individually using
  # for example bwa
  #
  # bwa mem mapping options:
  #       -A INT        score for a sequence match, which scales options -TdBOELU unless overridden [1]
  #       -B INT        penalty for a mismatch [4]
  #       -O INT[,INT]  gap open penalties for deletions and insertions [6,6]
  #       -E INT[,INT]  gap extension penalty; a gap of size k cost '{-O} + {-E}*k' [1,1] # this is set very high to avoid gaps
  #                                  at restriction sites. Setting the gap extension penalty high, produces better results as
  #                                  the sequences left and right of a restriction site are mapped independently.
  #       -L INT[,INT]  penalty for 5'- and 3'-end clipping [5,5] # this is set to no penalty.
  
  $ bwa mem -A1 -B4  -E50 -L0  index_path \
      -U mate_R1.fastq.gz 2>>mate_R1.log | samtools view -Shb - > mate_R1.bam
  
  $ bwa mem -A1 -B4  -E50 -L0  index_path \
      -U mate_R2.fastq.gz 2>>mate_R2.log | samtools view -Shb - > mate_R2.bam
  
  # build matrix from independently mated read pairs
  # the restriction sequence GATC is recognized by the DpnII restriction enzyme
  
  $ hicBuildMatrix --samFiles mate_R1.bam mate_R2.bam \
                   --binSize 10000 \
                   --restrictionSequence GATC \
                   --threads 4
                   --inputBufferSize 100000
                   --outBam hic.bam \
                   -o hic_matrix.h5
                   --QCfolder ./hicQC
  
  ```
  
  > **Comments:**
  > **username_3**: Option -SP exactly aims to disable "mapper specific heuristics designed for standard paired-end libraries".  
  
  ---
  
  