Iterative mapping

Iterative mapping is used when no restriction enzyme is used



In [1]:

    
r_enz = 'HindIII'



In [2]:

    
from pytadbit.mapping.full_mapper import full_mapping



In [ ]:

    
full_mapping('/media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem',
            nthreads=2, clean=True, r_enz=r_enz,
            frag_map=False, 
            out_map_dir='results/{0}/01_mapping/map{0}_r1'.format(r_enz),
            fastq_path='/media/storage/FASTQs/K562_{0}_1.fastq'.format(r_enz),
            windows=((1,25),(1,30),(1,40),(1,75)))

Fragment based mapping

Fragment based mapping is more precise, faster and with less background. It first look for the ligation site, then split the reads into the two frags and map each one individually. Then looks at the 75nt from the other side and if contiguous (short [<100ntdistance] from the right side of the first read) then combine it. If from other chr or location then adds an iteraction with that other fragment.



In [3]:

    
! mkdir -p results



In [4]:

    
! mkdir -p results/HindIII



In [5]:

    
! mkdir -p results/MboI



In [8]:

    
! mkdir -p results/HindIII/01_mapping



In [9]:

    
! mkdir -p results/MboI/01_mapping



In [6]:

    
ls









    



00_fastq_QC.ipynb       learning_python2.ipynb  results/
01_Mapping_Reads.ipynb  learning_python.ipynb



In [7]:

    
! ls









    



00_fastq_QC.ipynb	learning_python2.ipynb	results
01_Mapping_Reads.ipynb	learning_python.ipynb



In [13]:

    
! ls









    



00_fastq_QC.ipynb	learning_python2.ipynb	results
01_Mapping_Reads.ipynb	learning_python.ipynb



In [10]:

    
full_mapping('/media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem',
            nthreads=2, clean=True, r_enz=r_enz,
            frag_map=True, 
            out_map_dir='results/{0}/01_mapping/map{0}_r1'.format(r_enz),
            fastq_path='/media/storage/FASTQs/K562_{0}_1.fastq'.format(r_enz))









    



Preparing FASTQ file
  - conversion to MAP format
Mapping reads in window 1-end...
TO GEM /home/student/tmp/K562_HindIII_1_fu4k2a
/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 2 -i /home/student/tmp/K562_HindIII_1_fu4k2a -o /home/student/tmp/K562_HindIII_1_fu4k2a_full_1-end
Parsing result...
   x removing GEM input /home/student/tmp/K562_HindIII_1_fu4k2a
   x removing map /home/student/tmp/K562_HindIII_1_fu4k2a_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/tmp/K562_HindIII_1_fu4k2a_filt_1-end.map
Mapping fragments of remaining reads...
TO GEM /home/student/tmp/K562_HindIII_1_09JF3D
/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 2 -i /home/student/tmp/K562_HindIII_1_09JF3D -o /home/student/tmp/K562_HindIII_1_09JF3D_frag_1-end
Parsing result...
   x removing GEM input /home/student/tmp/K562_HindIII_1_09JF3D
   x removing failed to map /home/student/tmp/K562_HindIII_1_fu4k2a_fail.map
   x removing tmp mapped /home/student/tmp/K562_HindIII_1_09JF3D_frag_1-end.map






    Out[10]:





['results/HindIII/01_mapping/mapHindIII_r1/K562_HindIII_1_full_1-end.map',
 'results/HindIII/01_mapping/mapHindIII_r1/K562_HindIII_1_frag_1-end.map']



In [14]:

    
full_mapping('/media/storage/db/reference_genome/Homo_sapiens/hg38/hg38.gem',
            nthreads=2, clean=True, r_enz=r_enz,
            frag_map=True, 
            out_map_dir='results/{0}/01_mapping/map{0}_r2'.format(r_enz),
            fastq_path='/media/storage/FASTQs/K562_{0}_2.fastq'.format(r_enz))









    



/home/student/.miniconda2/lib/python2.7/site-packages/pytadbit/mapping/full_mapper.py:390: UserWarning: WARNING: only 169 Gb left on tmp_dir: /home/student/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/tmp/K562_HindIII_2_VPCLbv
/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 2 -i /home/student/tmp/K562_HindIII_2_VPCLbv -o /home/student/tmp/K562_HindIII_2_VPCLbv_full_1-end
Parsing result...
   x removing GEM input /home/student/tmp/K562_HindIII_2_VPCLbv
   x removing map /home/student/tmp/K562_HindIII_2_VPCLbv_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/tmp/K562_HindIII_2_VPCLbv_filt_1-end.map
Mapping fragments of remaining reads...
TO GEM /home/student/tmp/K562_HindIII_2_C3SDhB
/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 2 -i /home/student/tmp/K562_HindIII_2_C3SDhB -o /home/student/tmp/K562_HindIII_2_C3SDhB_frag_1-end
Parsing result...
   x removing GEM input /home/student/tmp/K562_HindIII_2_C3SDhB
   x removing failed to map /home/student/tmp/K562_HindIII_2_VPCLbv_fail.map
   x removing tmp mapped /home/student/tmp/K562_HindIII_2_C3SDhB_frag_1-end.map






    Out[14]:





['results/HindIII/01_mapping/mapHindIII_r2/K562_HindIII_2_full_1-end.map',
 'results/HindIII/01_mapping/mapHindIII_r2/K562_HindIII_2_frag_1-end.map']



In [15]:

    
ls









    



00_fastq_QC.ipynb       learning_python2.ipynb  results/
01_Mapping_Reads.ipynb  learning_python.ipynb



In [16]:

    
cd results/









    



/media/storage/Notebooks2/JCarlos/results



In [17]:

    
ls









    



HindIII/  MboI/



In [18]:

    
cd HindIII/









    



/media/storage/Notebooks2/JCarlos/results/HindIII



In [19]:

    
ls









    



01_mapping/



In [20]:

    
cd 01_mapping/









    



/media/storage/Notebooks2/JCarlos/results/HindIII/01_mapping



In [21]:

    
ls









    



mapHindIII_r1/  mapHindIII_r2/



In [22]:

    
cd mapHindIII_r1/









    



/media/storage/Notebooks2/JCarlos/results/HindIII/01_mapping/mapHindIII_r1



In [23]:

    
ls









    



K562_HindIII_1_frag_1-end.map  K562_HindIII_1_full_1-end.map



In [24]:

    
ll









    



total 1967960
-rw-rw-r-- 1 student  703591014 Apr  4 15:55 K562_HindIII_1_frag_1-end.map
-rw-rw-r-- 1 student 1311588316 Apr  4 15:51 K562_HindIII_1_full_1-end.map



In [25]:

    
! ls -lh









    



total 1.9G
-rw-rw-r-- 1 student student 671M Apr  4 15:55 K562_HindIII_1_frag_1-end.map
-rw-rw-r-- 1 student student 1.3G Apr  4 15:51 K562_HindIII_1_full_1-end.map



In [28]:

    
! head results/HindIII/01_mapping/mapHindIII_r2/K562_HindIII_1_full_1-end.map









    



head: cannot open ‘results/HindIII/01_mapping/mapHindIII_r2/K562_HindIII_1_full_1-end.map’ for reading: No such file or directory



In [29]:

    
! head K562_HindIII_1_full_1-end.map









    



NS500645:59:HCL32BGXY:1:11101:10953:1055	CTGACNTACAGGACAGTTATTTTCCAGACTGAAAATTCTCAAACCCAGTCGTGCAAATGATTCTTCCTGCAGTTG	AAAAA#EEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEAEEEEEEEEEEEEEEEEEEEE	0:1:0	chr10:-:2688535:5G69
NS500645:59:HCL32BGXY:1:11101:24830:1055	CTTTTNTTTGCTAGTCAAAGGCCTAGGCCTTTATCTTCCTTTTTGCTACAACCCTAAAGGCAAAGAAAGAAAGGG	AAAAA#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEE	0:1:0	chr1:-:50841578:5C69
NS500645:59:HCL32BGXY:1:11101:6039:1056	AGGGTNACATAGGTGACATAGCACGTAAGAGCTTGGGTCCTGAAGTTAGATGGGGCTGGGTTTCAAATCTATCCT	AAAAA#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE	0:0:1:0	chr7:+:129557395:5G65C3
NS500645:59:HCL32BGXY:1:11101:21852:1056	CTCTGNTTTGTTTCTTTGTTTGTTATGGTTGTTTTTCAGATAGATATAGTTTCATTACAAAAGAAAATTTGTTAT	AAAAA#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEE	0:0:1:0	chr6:-:23446287:5T68G
NS500645:59:HCL32BGXY:1:11101:19311:1056	ATTTTNTACCTTTTGAGAAACACTGAACTAATAATATTTCTGAGACCTCTTCAAGCTCTGACATGCTGTGTTTAG	AAAAA#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE	0:1:0	chr19:-:41066219:5A69
NS500645:59:HCL32BGXY:1:11101:17157:1056	ATAATNTATTTGCTCTCTTCAGCCTTTTTCCAGTTAGTACTCTCATCAGTAAGGCACATGAAGTGTGAAATTTTG	AAAAA#EEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE	0:1:0	chr9:+:114790228:5A69
NS500645:59:HCL32BGXY:1:11101:26187:1057	AAACTNATTTTTTTTTTTGCTTTTTTTTCCTTTTAAAAATATGAGCAGTAATCACATTGACTAAACAAACAGCCT	AAAAA#EEEEEEEEEEEEEEEEEEEEEEE/EEEEEEEEEEEEAEEEEEEEEEEEE6EEEEEEEAEEAEEE/EEEE	0:1:0	chr10:+:63481424:5G69
NS500645:59:HCL32BGXY:1:11101:8600:1058	AGAGANCCAAAGCTGGATAAATTTTCTACCTGAGAGCTGGAGTTTCCAAGAATGACAGACAGCCATGGGAAGAAG	AAAAA#EEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEAEEEE	0:1:0	chr5:+:53943740:5A69
NS500645:59:HCL32BGXY:1:11101:8731:1058	ATGAANCTACTACAAGAAAACATTAGGGAAACTCTTCAGGATATTGGTCTGGGCAAAAATTTCTTGAGTAATACC	AAAAA#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE	0:1:0	chr3:+:17104208:5A69
NS500645:59:HCL32BGXY:1:11101:23591:1059	CTATTNGTCTTAATAGTCTGAGTTCCCAAGATAAACAGAACATATACATATATTTATATTATAAGGAATTGGCTC	AAAAA#EEE/EEEEEEEEEEEE6E<EEEEE/EEEAEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEE/EAAEEEE	0:1:0	chr18:+:41009204:5A69



In [30]:

    
cd ..









    



/media/storage/Notebooks2/JCarlos/results/HindIII/01_mapping



In [31]:

    
cd ..









    



/media/storage/Notebooks2/JCarlos/results/HindIII



In [32]:

    
cd ..









    



/media/storage/Notebooks2/JCarlos/results



In [33]:

    
cd ..









    



/media/storage/Notebooks2/JCarlos



In [36]:

    
! head results/HindIII/01_mapping/mapHindIII_r2/K562_HindIII_2_frag_1-end.map









    



NS500645:59:HCL32BGXY:1:11101:24830:1055~2~	AAGCTTTTTTTTCCCTTTCTTTCTTTGCCTTTAGGGTTGTAGCAAAAAGGAAGATAA	HHHHHHEEEEEEEAEEEEEEEEEEEEEEEEEEAEEEEAEEEEE/EEE6EEEEEEEEE	1	chr1:+:50841566:57
NS500645:59:HCL32BGXY:1:11101:6039:1056	GGCTTCTTCTCTTTAAGGCAACATAGTATAAAGCTT	AAAAAEEEEEEEEEEEEEEEEEEE/EEAE6HHHHHH	1	chr7:+:128976029:36
NS500645:59:HCL32BGXY:1:11101:6039:1056~2~	AAGCTTTCGAAGGCTAGCTAGGATAGATTTGAAACCCAGCC	HHHHHHEEEEEEEEEEEEEEEEEEEEEEEEEEE/EEEAEEE	0:1	chr7:-:129557448:22G18
NS500645:59:HCL32BGXY:1:11101:19311:1056	ACCTACTTCATAAAAAGCCAAGGATGTAATTTTTTCAGTTTAGTGTAAGCTT	/AAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEHHHHHH	1	chr19:+:41557056:52
NS500645:59:HCL32BGXY:1:11101:19311:1056~2~	AAGCTTAGGGTTCAAGTCTAAACAC	HHHHHHEEEEEEAEEEEEEEEEEEE	1	chr19:+:41066202:25
NS500645:59:HCL32BGXY:1:11101:17064:1056	CAACAATCATATATTTCCTGCAAATGGAAAGGCAAATGGTACAGAGTCCCCAATGTACAAGCTT	AAAAAEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEAHHHHHH	0:1	chr4:+:33593268:40C23
NS500645:59:HCL32BGXY:1:11101:26187:1057	CCCACTTACCACTGAAGTTAAAGCTT	AAAAAEEEA<<AEEEEEEEEHHHHHH	1	chr10:+:63502264:26
NS500645:59:HCL32BGXY:1:11101:26187:1057~2~	AAGCTTCTAAGTGAAGTCTCAAAGCCACAGAATTAAGAAAGGCTGTTTGTT	HHHHHHAEA6EEAAEEE<E/<EA/6/A6A</EA<A/A/AA//<</AA<//E	1	chr10:-:63481487:51
NS500645:59:HCL32BGXY:1:11101:8250:1059	AGCCCGCTCAAGGTCTAAACTCGGCTCGTCCCTGGGCCTTCCACCTACAAGAAGCTT	AAAAAE666EEEEEEEEEAAAEAEEEA/EEEEEAEEE6EEEE/EEEE6EAEHHHHHH	0:1	chr13:+:29695513:3A53
NS500645:59:HCL32BGXY:1:11101:17211:1059	AGTCATAGCAACAACAGATACTTAAGCTT	AAAAAEEEEEEEEEEEEEEEEEEHHHHHH	1	chr2:+:25493831:29



In [ ]: