Iterative mapping


In [1]:
r_enz ='HindIII'

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

In [7]:
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)))


---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
<ipython-input-7-1ed3f35e1fea> in <module>()
      3              out_map_dir='results/{0}/01_mapping/map{0}_r1'.format(r_enz),
      4              fastq_path='/media/storage/FASTQs/K562_{0}_1.fastq'.format(r_enz),
----> 5              windows=((1.,25), (1,30), (1,40), (1, 75)))

/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)
    384     # create directories
    385     for rep in [temp_dir, out_map_dir]:
--> 386         mkdir(rep)
    387     # check space
    388     fspace = int(get_free_space_mb(temp_dir, div=3))

/home/student/.miniconda2/lib/python2.7/site-packages/pytadbit/utils/file_handling.pyc in mkdir(dnam)
    100 def mkdir(dnam):
    101     try:
--> 102         os.mkdir(dnam)
    103     except OSError as exc:
    104         if exc.errno == errno.EEXIST and os.path.isdir(dnam):

OSError: [Errno 2] No such file or directory: 'results/HindIII/01_mapping/mapHindIII_r1'

In [ ]:


In [ ]:


In [21]:
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='/media/storage/Notebooks2/Debora/results/HindIII/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)))


Preparing FASTQ file
  - conversion to MAP format
  - trimming reads 1-25
Mapping reads in window 1-25...
TO GEM /home/student/tmp/K562_HindIII_1_493SaP
/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_493SaP -o /home/student/tmp/K562_HindIII_1_493SaP_full_1-25
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-21-f45631558063> in <module>()
      3              out_map_dir='/media/storage/Notebooks2/Debora/results/HindIII/01_mapping/map{0}_r1'.format(r_enz),
      4              fastq_path='/media/storage/FASTQs/K562_{0}_1.fastq'.format(r_enz),
----> 5              windows=((1,25), (1,30), (1,40), (1, 75)))

/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)
    435 
    436         if not skip:
--> 437             gem_mapping(gem_index_path, curr_map, out_map_path, **kwargs)
    438             # parse map file to extract not uniquely mapped reads
    439             print 'Parsing result...'

/home/student/.miniconda2/lib/python2.7/site-packages/pytadbit/mapping/full_mapper.pyc in gem_mapping(gem_index_path, fastq_path, out_map_path, gem_binary, **kwargs)
    328     try:
    329         # check_call(gem_cmd, stdout=PIPE, stderr=PIPE)
--> 330         out, err = Popen(gem_cmd, stdout=PIPE, stderr=PIPE).communicate()
    331     except CalledProcessError as e:
    332         print out

/home/student/.miniconda2/lib/python2.7/subprocess.pyc in communicate(self, input)
    477             return (stdout, stderr)
    478 
--> 479         return self._communicate(input)
    480 
    481 

/home/student/.miniconda2/lib/python2.7/subprocess.pyc in _communicate(self, input)
   1096 
   1097             if _has_poll:
-> 1098                 stdout, stderr = self._communicate_with_poll(input)
   1099             else:
   1100                 stdout, stderr = self._communicate_with_select(input)

/home/student/.miniconda2/lib/python2.7/subprocess.pyc in _communicate_with_poll(self, input)
   1150             while fd2file:
   1151                 try:
-> 1152                     ready = poller.poll()
   1153                 except select.error, e:
   1154                     if e.args[0] == errno.EINTR:

KeyboardInterrupt: 

Fragment based mapping


In [8]:
! mkdir -p results

In [9]:
! mkdir -p results/HindIII

In [10]:
! mkdir -p results/MboI

In [11]:
ls


00_fastq_QC.ipynb  01_mapping.ipynb  learning_python.ipynb  results/

In [12]:
! mkdir -p results/HindIII/01_mapping

In [13]:
! mkdir -p results/MboI/01_mapping

In [19]:
cd results/HindIII/01_mapping/


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

In [34]:
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='/media/storage/Notebooks2/Debora/results/HindIII/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 150 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_LLm74P
/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_LLm74P -o /home/student/tmp/K562_HindIII_2_LLm74P_full_1-end
Parsing result...
   x removing GEM input /home/student/tmp/K562_HindIII_2_LLm74P
   x removing map /home/student/tmp/K562_HindIII_2_LLm74P_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_LLm74P_filt_1-end.map
Mapping fragments of remaining reads...
TO GEM /home/student/tmp/K562_HindIII_2_uGlYse
/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_uGlYse -o /home/student/tmp/K562_HindIII_2_uGlYse_frag_1-end
Parsing result...
   x removing GEM input /home/student/tmp/K562_HindIII_2_uGlYse
   x removing failed to map /home/student/tmp/K562_HindIII_2_LLm74P_fail.map
   x removing tmp mapped /home/student/tmp/K562_HindIII_2_uGlYse_frag_1-end.map
Out[34]:
['/media/storage/Notebooks2/Debora/results/HindIII/01_mapping/mapHindIII_r2/K562_HindIII_2_full_1-end.map',
 '/media/storage/Notebooks2/Debora/results/HindIII/01_mapping/mapHindIII_r2/K562_HindIII_2_frag_1-end.map']

In [40]:
! head /media/storage/Notebooks2/Debora/results/HindIII/01_mapping/mapHindIII_r2/K562_HindIII_2_full_1-end.map


NS500645:59:HCL32BGXY:1:11101:14163:1054	CAACCAATTGTCAAACAGAAGAATTTTAAATCTACCTATAAGCTGGAAGCCCCACCCCACTTCAAGTTGTCCCAT	AAAAAEEE/EEEEEEEEEEEEEEEEEEEEEEEEEEEEAEE6EEEEEAAEE<EE/<EEEEEEEEEEEEEEE<A/<<	1	chr2:-:62077857:75
NS500645:59:HCL32BGXY:1:11101:14193:1055	CAGGCTGGTCTCGAACTCCTGGGCTCAAGCGATCCACCCACCTCAGCCTCCCAAAGTGCTGGGATTACAGGCGTG	AAAAAEEEEEEEE<AEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEE/EAAEEE	1	chr16:+:14833322:75
NS500645:59:HCL32BGXY:1:11101:10953:1055	GATTCTCAAATTATCTATAACAGATCCATCTTCAGAATCTAAGCACAGAGGGAAGCGATCGGGACATACAAGGGC	AAAAAEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEE	1	chr10:-:2392912:75
NS500645:59:HCL32BGXY:1:11101:16360:1056	GTCACAGCCAGTTGGGAAAACTGGCTAGCCATATGCAGAAAACTGAAACTGGACCCCATCCTTACACCTTGCACA	AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE<EEEEEEEE6EEEEAEEEEAEEEEEEEAEEEEEEEEA	0:1	chr2:+:7141515:3G71
NS500645:59:HCL32BGXY:1:11101:21852:1056	CAGACAGACAAATACTGAAAGAATTCACAACTACCAAGCCAATACTACCAACACTGCTAAAAGGAGCTCTAAATC	AAAAAEEEEEEEEEEEEEEEEEEEEEEE<EEEEEEE<EEAAEEEEEEEEAEEEEEEEEEEEEEEAEEEEEEEEE<	1	chr6:+:23863609:75
NS500645:59:HCL32BGXY:1:11101:17157:1056	CCTCATGGGGACCAAGTCTGTGTGCGAAGTAGGTAGCAACTGGTTCCAGCCCATCTACCTCGGAGCCATGTTCTC	AAAAAEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEAEEEEEEAEEEEEEEEEEA	0:1	chr9:-:114790552:21A53
NS500645:59:HCL32BGXY:1:11101:11342:1057	CTGCTCCTGTGAGTCAGAGTGTGTCATTTCCTCACCTAAAACACTCCAGTGGCTCCACCTCGGTCTTGTGAAGCT	AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEE	1	chr6:-:29924856:75
NS500645:59:HCL32BGXY:1:11101:15358:1057	CATGTTTATGTACTCATATGTATATATGTGTATTTGTGTGTATATATACAGACACTCTGATACACATGCACAAAC	AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE	1	chr2:+:58688316:75
NS500645:59:HCL32BGXY:1:11101:8600:1058	CTATTCTTTATTGTAAGATATCAGAACAGCAGAGCTCAAAGTGCACACGATTTCCCTTCTTATGTCAGGCACTGG	AAAAAAEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEA/EEEEEEEEEEEEEAAEEEE	1	chr5:-:5504049:75
NS500645:59:HCL32BGXY:1:11101:23591:1059	CTTGGCTCATCATCTCAGAGTGTTTGAATATTGAATTTCAAGTAGGAAAGCAGTAGACTGCAACAAGACTAAAAC	//AAA/EEE//<EEEEEEEEEEEEEEEEE/EEE</<EEAE6</EEEE/EEEEE/EAEE<EEEAEEEE/A/EEEEE	0:1	chr18:-:37712209:29C45

In [41]:
!head /media/storage/Notebooks2/Debora/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 [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]: