Iterative vs fragment-based mapping

Iterative mapping first proposed by (Imakaev et al., 2012), allows to map usually a high number of reads. However other methodologies, less "brute-force" can be used to take into account the chimeric nature of Hi-C reads.

A simple alternative is to allow split mapping, just as with RNA-seq data.

Another way consists in pre-truncating (Ay and Noble, 2015) reads that contain a ligation site and map only the longest part of the read (Wingett et al., 2015).

Finally, an intermediate approach, fragment-based, consists in mapping full length reads first, and than splitting unmapped reads at the ligation sites (Serra, Ba`{u, Filion and Marti-Renom, 2016).

Advantages of iterative mapping

  • It's the only solution when no restriction enzyme has been used (i.e. micro-C)
  • Can be faster when few windows (2 or 3) are used

Advantages of fragment-based mapping

  • Generally faster
  • Safer: mapped reads are generally larger than 25-30 nm (the largest window used in iterative mapping). Less reads are mapped, but the difference is usually canceled or reversed when looking for "valid-pairs".

Note: We use GEM (Marco-Sola, Sammeth, Guig\'{o and Ribeca, 2012), performance are very similar to Bowtie2, perhaps a bit better.

For now TADbit is only compatible with GEM.

Mapping


In [1]:
from pytadbit.mapping.full_mapper import full_mapping

The full mapping function can be used to perform either iterative or fragment-based mapping, or a combination of both.

Iterative mapping

Here an example of use as iterative mapping:


In [11]:
r_enz = 'HindIII'

In [ ]:
! mkdir -p results/iterativ/$r_enz
! mkdir -p results/iterativ/$r_enz/01_mapping

In [6]:
# for the first side of the reads
full_mapping(gem_index_path='/media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem',
             out_map_dir='results/iterativ/{0}/01_mapping/mapped_{0}_r1/'.format(r_enz),
             fastq_path='/media/storage/FASTQs/K562_%s_1.fastq' % (r_enz),
             r_enz='hindIII', frag_map=False, clean=True, nthreads=20,
             windows=((1,25),(1,30),(1,35),(1,40),(1,45),(1,50),(1,55),(1,60),(1,65),(1,70),(1,75)), 
             temp_dir='results/iterativ/{0}/01_mapping/mapped_{0}_r1_tmp/'.format(r_enz))


Preparing FASTQ file
  - conversion to MAP format
  - trimming reads 1-25
Mapping reads in window 1-25...
TO GEM /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_DtQxPB
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_DtQxPB -o /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_DtQxPB_full_1-25
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_DtQxPB
   x removing map /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_DtQxPB_full_1-25.map
Preparing MAP file
  - trimming reads 1-30
   x removing original input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_DtQxPB_filt_1-25.map
Mapping reads in window 1-30...
TO GEM /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_mhZLf6
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_mhZLf6 -o /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_mhZLf6_full_1-30
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_mhZLf6
   x removing map /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_mhZLf6_full_1-30.map
Preparing MAP file
  - trimming reads 1-35
   x removing original input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_mhZLf6_filt_1-30.map
Mapping reads in window 1-35...
TO GEM /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_wezLHr
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_wezLHr -o /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_wezLHr_full_1-35
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_wezLHr
   x removing map /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_wezLHr_full_1-35.map
Preparing MAP file
  - trimming reads 1-40
   x removing original input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_wezLHr_filt_1-35.map
Mapping reads in window 1-40...
TO GEM /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_1o4Mbp
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_1o4Mbp -o /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_1o4Mbp_full_1-40
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_1o4Mbp
   x removing map /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_1o4Mbp_full_1-40.map
Preparing MAP file
  - trimming reads 1-45
   x removing original input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_1o4Mbp_filt_1-40.map
Mapping reads in window 1-45...
TO GEM /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_jm52Qw
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_jm52Qw -o /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_jm52Qw_full_1-45
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_jm52Qw
   x removing map /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_jm52Qw_full_1-45.map
Preparing MAP file
  - trimming reads 1-50
   x removing original input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_jm52Qw_filt_1-45.map
Mapping reads in window 1-50...
TO GEM /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_n2STO3
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_n2STO3 -o /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_n2STO3_full_1-50
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_n2STO3
   x removing map /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_n2STO3_full_1-50.map
Preparing MAP file
  - trimming reads 1-55
   x removing original input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_n2STO3_filt_1-50.map
Mapping reads in window 1-55...
TO GEM /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_oTOw3K
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_oTOw3K -o /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_oTOw3K_full_1-55
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_oTOw3K
   x removing map /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_oTOw3K_full_1-55.map
Preparing MAP file
  - trimming reads 1-60
   x removing original input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_oTOw3K_filt_1-55.map
Mapping reads in window 1-60...
TO GEM /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_13gsmu
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_13gsmu -o /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_13gsmu_full_1-60
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_13gsmu
   x removing map /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_13gsmu_full_1-60.map
Preparing MAP file
  - trimming reads 1-65
   x removing original input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_13gsmu_filt_1-60.map
Mapping reads in window 1-65...
TO GEM /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_jcqagN
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_jcqagN -o /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_jcqagN_full_1-65
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_jcqagN
   x removing map /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_jcqagN_full_1-65.map
Preparing MAP file
  - trimming reads 1-70
   x removing original input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_jcqagN_filt_1-65.map
Mapping reads in window 1-70...
TO GEM /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_MVCVk5
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_MVCVk5 -o /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_MVCVk5_full_1-70
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_MVCVk5
   x removing map /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_MVCVk5_full_1-70.map
Preparing MAP file
  - trimming reads 1-75
   x removing original input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_MVCVk5_filt_1-70.map
Mapping reads in window 1-75...
TO GEM /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_w13Ngf
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_w13Ngf -o /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_w13Ngf_full_1-75
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_w13Ngf
   x removing map /home/student/Course/3DAROC/Notebooks/results/iterativ/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_w13Ngf_full_1-75.map
Out[6]:
['results/iterativ/HindIII/01_mapping/mapped_HindIII_r1/K562_HindIII_1_full_1-25.map',
 'results/iterativ/HindIII/01_mapping/mapped_HindIII_r1/K562_HindIII_1_full_1-30.map',
 'results/iterativ/HindIII/01_mapping/mapped_HindIII_r1/K562_HindIII_1_full_1-35.map',
 'results/iterativ/HindIII/01_mapping/mapped_HindIII_r1/K562_HindIII_1_full_1-40.map',
 'results/iterativ/HindIII/01_mapping/mapped_HindIII_r1/K562_HindIII_1_full_1-45.map',
 'results/iterativ/HindIII/01_mapping/mapped_HindIII_r1/K562_HindIII_1_full_1-50.map',
 'results/iterativ/HindIII/01_mapping/mapped_HindIII_r1/K562_HindIII_1_full_1-55.map',
 'results/iterativ/HindIII/01_mapping/mapped_HindIII_r1/K562_HindIII_1_full_1-60.map',
 'results/iterativ/HindIII/01_mapping/mapped_HindIII_r1/K562_HindIII_1_full_1-65.map',
 'results/iterativ/HindIII/01_mapping/mapped_HindIII_r1/K562_HindIII_1_full_1-70.map',
 'results/iterativ/HindIII/01_mapping/mapped_HindIII_r1/K562_HindIII_1_full_1-75.map']

And for the second side of the read:


In [8]:
# for the second side of the reads
full_mapping(gem_index_path='/media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem',
             out_map_dir='results/iterativ/{0}/01_mapping/mapped_{0}_r2/'.format(r_enz),
             fastq_path='/media/storage/FASTQs/K562_%s_2.fastq' % (r_enz),
             r_enz=r_enz, frag_map=False,  clean=True, nthreads=20,
             windows=((1,25),(1,30),(1,35),(1,40),(1,45),(1,50),(1,55),(1,60),(1,65),(1,70),(1,75)),
             temp_dir='results/iterativ/{0}/01_mapping/mapped_{0}_r2_tmp/'.format(r_enz))


Preparing FASTQ file
  - conversion to MAP format
  - trimming reads 1-25
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-8-0c7087896b84> in <module>()
      5              r_enz=r_enz, frag_map=False,  clean=True, nthreads=20,
      6              windows=((1,25),(1,30),(1,35),(1,40),(1,45),(1,50),(1,55),(1,60),(1,65),(1,70),(1,75)),
----> 7              temp_dir='results/iterativ/{0}/01_mapping/mapped_{0}_r2_tmp/'.format(r_enz))

/home/student/.miniconda2/lib/python2.7/site-packages/pytadbit/mapping/full_mapper.pyc in full_mapping(gem_index_path, fastq_path, out_map_dir, r_enz, frag_map, min_seq_len, windows, add_site, clean, get_nread, **kwargs)
    418                    or input_reads.endswith('.dsrc'    )),
    419             min_seq_len=min_seq_len, trim=win, skip=skip, nthreads=nthreads,
--> 420             light_storage=light_storage)
    421         # clean
    422         if input_reads != fastq_path and clean:

/home/student/.miniconda2/lib/python2.7/site-packages/pytadbit/mapping/full_mapper.pyc in transform_fastq(fastq_path, out_fastq, trim, r_enz, add_site, min_seq_len, fastq, verbose, light_storage, **kwargs)
    183             except StopIteration:
    184                 continue
--> 185         out.write(_map2fastq('\t'.join((insert_mark(header, cnt),
    186                                         seq, qal, '0', '-\n'))))
    187         # the next fragments should be preceded by the RE site

KeyboardInterrupt: 

Fragment-based mapping

With fragment based mapping it would be:


In [12]:
! mkdir -p results/fragment/$r_enz
! mkdir -p results/fragment/$r_enz/01_mapping

In [13]:
# for the first side of the reads 
full_mapping(gem_index_path='/media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem',
             out_map_dir='results/fragment/{0}/01_mapping/mapped_{0}_r1/'.format(r_enz),
             fastq_path='/media/storage/FASTQs/K562_%s_1.fastq' % (r_enz),
             r_enz=r_enz, frag_map=True, clean=True, nthreads=20, 
             temp_dir='results/fragment/{0}/01_mapping/mapped_{0}_r1_tmp/'.format(r_enz))


/home/student/.miniconda2/lib/python2.7/site-packages/pytadbit/mapping/full_mapper.py:390: UserWarning: WARNING: only 61 Gb left on tmp_dir: /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r1_tmp

  warn('WARNING: only %d Gb left on tmp_dir: %s\n' % (fspace, temp_dir))
Preparing FASTQ file
  - conversion to MAP format
Mapping reads in window 1-end...
TO GEM /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_ivD4Lh
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_ivD4Lh -o /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_ivD4Lh_full_1-end
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_ivD4Lh
   x removing map /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_ivD4Lh_full_1-end.map
  - splitting into restriction enzyme (RE) fragments using ligation sites
  - ligation sites are replaced by RE sites to match the reference genome
    * enzyme: HindIII, ligation site: AAGCTAGCTT, RE site: AAGCTT
Preparing MAP file
   x removing pre-GEM input /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_ivD4Lh_filt_1-end.map
Mapping fragments of remaining reads...
TO GEM /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_hRD31s
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_hRD31s -o /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_hRD31s_frag_1-end
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_hRD31s
   x removing failed to map /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_ivD4Lh_fail.map
   x removing tmp mapped /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r1_tmp/K562_HindIII_1_hRD31s_frag_1-end.map
Out[13]:
['results/fragment/HindIII/01_mapping/mapped_HindIII_r1/K562_HindIII_1_full_1-end.map',
 'results/fragment/HindIII/01_mapping/mapped_HindIII_r1/K562_HindIII_1_frag_1-end.map']

In [14]:
# for the second side of the reads
full_mapping(gem_index_path='/media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem',
             out_map_dir='results/fragment/{0}/01_mapping/mapped_{0}_r2/'.format(r_enz),
             fastq_path='/media/storage/FASTQs/K562_%s_2.fastq' % (r_enz),
             r_enz=r_enz, frag_map=True, clean=True, nthreads=20, 
             temp_dir='results/fragment/{0}/01_mapping/mapped_{0}_r2_tmp/'.format(r_enz))


/home/student/.miniconda2/lib/python2.7/site-packages/pytadbit/mapping/full_mapper.py:390: UserWarning: WARNING: only 41 Gb left on tmp_dir: /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r2_tmp

  warn('WARNING: only %d Gb left on tmp_dir: %s\n' % (fspace, temp_dir))
Preparing FASTQ file
  - conversion to MAP format
Mapping reads in window 1-end...
TO GEM /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r2_tmp/K562_HindIII_2_BkYlKm
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r2_tmp/K562_HindIII_2_BkYlKm -o /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r2_tmp/K562_HindIII_2_BkYlKm_full_1-end
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r2_tmp/K562_HindIII_2_BkYlKm
   x removing map /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r2_tmp/K562_HindIII_2_BkYlKm_full_1-end.map
  - splitting into restriction enzyme (RE) fragments using ligation sites
  - ligation sites are replaced by RE sites to match the reference genome
    * enzyme: HindIII, ligation site: AAGCTAGCTT, RE site: AAGCTT
Preparing MAP file
   x removing pre-GEM input /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r2_tmp/K562_HindIII_2_BkYlKm_filt_1-end.map
Mapping fragments of remaining reads...
TO GEM /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r2_tmp/K562_HindIII_2_98v6S7
/usr/local/bin/gem-mapper -I /media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 20 -i /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r2_tmp/K562_HindIII_2_98v6S7 -o /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r2_tmp/K562_HindIII_2_98v6S7_frag_1-end
Parsing result...
   x removing GEM input /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r2_tmp/K562_HindIII_2_98v6S7
   x removing failed to map /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r2_tmp/K562_HindIII_2_BkYlKm_fail.map
   x removing tmp mapped /home/student/Course/3DAROC/Notebooks/results/fragment/HindIII/01_mapping/mapped_HindIII_r2_tmp/K562_HindIII_2_98v6S7_frag_1-end.map
Out[14]:
['results/fragment/HindIII/01_mapping/mapped_HindIII_r2/K562_HindIII_2_full_1-end.map',
 'results/fragment/HindIII/01_mapping/mapped_HindIII_r2/K562_HindIII_2_frag_1-end.map']

References

[^](#ref-1) Imakaev, Maxim V and Fudenberg, Geoffrey and McCord, Rachel Patton and Naumova, Natalia and Goloborodko, Anton and Lajoie, Bryan R and Dekker, Job and Mirny, Leonid A. 2012. Iterative correction of Hi-C data reveals hallmarks of chromosome organization.. URL

[^](#ref-2) Ay, Ferhat and Noble, William Stafford. 2015. Analysis methods for studying the 3D architecture of the genome. URL

[^](#ref-3) Wingett, Steven and Ewels, Philip and Furlan-Magaril, Mayra and Nagano, Takashi and Schoenfelder, Stefan and Fraser, Peter and Andrews, Simon. 2015. HiCUP: pipeline for mapping and processing Hi-C data.. URL

[^](#ref-4) Serra, Fran\c{cois and Ba`{u, Davide and Filion, Guillaume and Marti-Renom, Marc A.. 2016. Structural features of the fly chromatin colors revealed by automatic three-dimensional modeling.. URL

[^](#ref-5) Marco-Sola, Santiago and Sammeth, Michael and Guig\'{o, Roderic and Ribeca, Paolo. 2012. The GEM mapper: fast, accurate and versatile alignment by filtration.