6
$\begingroup$

Library Prep

  • I have a Hi-C 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 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 Heng Li discusses that the -5 flag behavior has changed, as well as a new -q flag. In this post 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.
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Devon Ryan Sep 22 '18 at 19:39
4
$\begingroup$

See here, in particular slide #10 from the tutorial:

bwa mem -5SPM ref.fa read1.fq read2.fq > out.sam

Here, -SP disables pairing. The Aiden lab and this paper also use a similar command line. Beware that there are many Hi-C pipelines, but not many are using bwa-mem.

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – conchoecia Oct 3 '18 at 23:18
1
$\begingroup$

There are some interesting suggestions from the HiCExplorer developers:

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
| improve this answer | |
$\endgroup$
  • 2
    $\begingroup$ Option -SP exactly aims to disable "mapper specific heuristics designed for standard paired-end libraries". $\endgroup$ – user172818 Sep 27 '18 at 1:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.