What is the correct way to map Hi-C data with bwa mem?

Question

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.

Answer 1

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.

Answer 2

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