ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/tophat
(Generate patch)
# Line 176 | Line 176
176   fail_str = "\t[FAILED]\n"
177   gtf_juncs = None #file name with junctions extracted from given GFF file
178  
179 + #mapping types:
180  
181  
182 + _reads_vs_G, _reads_vs_T, _segs_vs_G, _segs_vs_J = range(1,5)
183 +
184   def init_logger(log_fname):
185      global tophat_logger
186      tophat_logger = logging.getLogger('project')
# Line 1573 | Line 1576
1576                self.popen=subprocess.Popen(pipecmd,
1577                      preexec_fn=subprocess_setup,
1578                      stdin=self.fsrc,
1579 <                    stdout=subprocess.PIPE, close_fds=True)
1579 >                    stdout=subprocess.PIPE, stderr=tophat_log, close_fds=True)
1580             except Exception:
1581                die("Error: could not open pipe "+' '.join(pipecmd)+' < '+ self.fname)
1582             self.file=self.popen.stdout
# Line 1596 | Line 1599
1599               self.popen=subprocess.Popen(pipecmd,
1600                     preexec_fn=subprocess_setup,
1601                     stdin=subprocess.PIPE,
1602 <                   stdout=self.ftarget, close_fds=True)
1602 >                   stderr=tophat_log, stdout=self.ftarget, close_fds=True)
1603            except Exception:
1604                die("Error: could not open writer pipe "+' '.join(pipecmd)+' < '+ self.fname)
1605            self.file=self.popen.stdin # client writes to this end of the pipe
# Line 1734 | Line 1737
1737                 raise Exception()
1738             except Exception, e:
1739               die(fail_str+"Error retrieving prep_reads info.")
1740 <           th_logp("\t%s reads: min. length=%s, count=%s" %  (side, self.min_len, self.out_count))
1740 >           th_logp("\t%5s reads: min. length=%s, count=%s" %  (side, self.min_len, self.out_count))
1741             if self.min_len < 20:
1742                 th_logp("Warning: short reads (<20bp) will make TopHat quite slow and take large amount of memory because they are likely to be mapped to too many places")
1743  
# Line 1808 | Line 1811
1811              zip_proc=subprocess.Popen(zip_cmd,
1812                                    preexec_fn=subprocess_setup,
1813                                    stdin=filter_proc.stdout,
1814 <                                  stdout=kept_reads)
1814 >                                  stderr=tophat_log, stdout=kept_reads)
1815              filter_proc.stdout.close() #as per http://bugs.python.org/issue7678
1816              zip_proc.communicate()
1817              retcode=filter_proc.poll()
# Line 1838 | Line 1841
1841             mapped_reads,
1842             unmapped_reads,
1843             extra_output = "",
1844 <           mapping_type = "read_against_genome",
1844 >           mapping_type = _reads_vs_G,
1845             multihits_out = None): #only --prefilter-multihits should activate this parameter for the initial prefilter search
1846      start_time = datetime.now()
1847      bwt_idx_name = bwt_idx_prefix.split('/')[-1]
# Line 1846 | Line 1849
1849      readfile_basename=getFileBaseName(reads_file)
1850  
1851      g_mapping, t_mapping, seg_mapping = False, False, False
1852 <    if mapping_type == "read_against_transcriptome":
1852 >    if mapping_type == _reads_vs_T:
1853          t_mapping = True
1854 <    elif mapping_type == "seg_against_genome" or mapping_type == "seg_against_junctions":
1854 >    elif mapping_type >= _segs_vs_G:
1855          seg_mapping = True
1856      else:
1857          g_mapping = True
# Line 1925 | Line 1928
1928                   signal.signal(signal.SIGTERM, on_sig_exit)
1929                   subprocess.call(unm_zipcmd,
1930                                   stdin=open(unmapped_reads_fifo, "r"),
1931 +                                 stderr=tophat_log,
1932                                   stdout=open(unmapped_reads_out, "wb"))
1933                   os._exit(os.EX_OK)
1934  
# Line 2066 | Line 2070
2070             if z_input:
2071                unzip_proc = subprocess.Popen(unzip_cmd,
2072                                       stdin=open(reads_file, "rb"),
2073 <                                     stdout=subprocess.PIPE)
2073 >                                     stderr=tophat_log, stdout=subprocess.PIPE)
2074                shellcmd=' '.join(unzip_cmd) + "< " +reads_file +"|"
2075             else:
2076                 #must be uncompressed fastq input (unmapped reads from a previous run)
2077                 #or a BAM file with unmapped reads
2078                 if bam_input:
2079 <                   unzip_proc = subprocess.Popen(unzip_cmd, stdout=subprocess.PIPE)
2079 >                   unzip_proc = subprocess.Popen(unzip_cmd, stderr=tophat_log, stdout=subprocess.PIPE)
2080                     shellcmd=' '.join(unzip_cmd) + "|"
2081                 else:
2082                     bowtie_cmd += [reads_file]
# Line 2093 | Line 2097
2097             shellcmd += "|"+ ' '.join(samzip_cmd)+" > " + mapped_reads
2098             fix_order_proc = subprocess.Popen(fix_map_cmd,
2099                                               stdin=bowtie_proc.stdout,
2100 +                                             stderr=tophat_log,
2101                                               stdout=subprocess.PIPE)
2102             bowtie_proc.stdout.close() # needed?
2103             samzip_proc = subprocess.Popen(samzip_cmd,
2104                                   preexec_fn=subprocess_setup,
2105                                   stdin=fix_order_proc.stdout,
2106 +                                 stderr=tophat_log,
2107                                   stdout=open(mapped_reads, "wb"))
2108             fix_order_proc.stdout.close()
2109          else:
2110             #write BAM output directly
2111             fix_order_proc = subprocess.Popen(fix_map_cmd,
2112 <                                          stdin=bowtie_proc.stdout)
2112 >                                          stdin=bowtie_proc.stdout,
2113 >                                          stderr=tophat_log)
2114                                          # stdout=open(mapped_reads, "w"))
2115             bowtie_proc.stdout.close()
2116  
# Line 2244 | Line 2251
2251                      fusions_file_list,
2252                      reference_fasta]
2253      try:
2254 <        print >> run_log, " ".join(juncs_db_cmd)
2254 >        print >> run_log, " ".join(juncs_db_cmd) + " > " + external_splices_out_name
2255          retcode = subprocess.call(juncs_db_cmd,
2256                                   stderr=juncs_db_log,
2257                                   stdout=external_splices_out)
# Line 2912 | Line 2919
2919                                             params.t_mismatches,
2920                                             mapped_gtf_out,
2921                                             unmapped_gtf,
2922 <                                           "", "read_against_transcriptome")
2922 >                                           "", _reads_vs_T)
2923          mapped_gtf_list.append(mapped_gtf_map)
2924          unmapped_gtf_list.append(unmapped)
2925  
# Line 2966 | Line 2973
2973           zip_proc = subprocess.Popen(zip_cmd,
2974                                 preexec_fn=subprocess_setup,
2975                                 stdin=prep_proc.stdout,
2976 <                               stdout=um_reads)
2976 >                               stderr=tophat_log, stdout=um_reads)
2977           prep_proc.stdout.close() #as per http://bugs.python.org/issue7678
2978           zip_proc.communicate()
2979           retcode=prep_proc.poll()
# Line 3097 | Line 3104
3104                                                     unspliced_out,
3105                                                     unmapped_unspliced,
3106                                                     "",
3107 <                                                   "read_against_genome")
3107 >                                                   _reads_vs_G)
3108  
3109          # daehwan - check this out
3110          # params.bowtie2 = False
# Line 3122 | Line 3129
3129              for i in range(len(read_segments)):
3130                  seg = read_segments[i]
3131                  fbasename=getFileBaseName(seg)
3132 <                seg_out =  tmp_dir + fbasename + ".bwtout"
3133 <                if use_zpacker: seg_out += ".z"
3132 >                seg_out =  tmp_dir + fbasename
3133 >                #if use_zpacker: seg_out += ".z"
3134                  unmapped_seg = tmp_dir + fbasename + "_unmapped.fq"
3135                  if use_BWT_FIFO:
3136                      unmapped_seg += ".z"
# Line 3136 | Line 3143
3143                                               seg_out,
3144                                               unmapped_seg,
3145                                               extra_output,
3146 <                                             "seg_against_genome")
3146 >                                             _segs_vs_G)
3147                  seg_maps.append(seg_map)
3148                  unmapped_segs.append(unmapped)
3149                  segs.append(seg)
# Line 3252 | Line 3259
3259              for seg in maps[ri].segs:
3260                  #search each segment
3261                  fsegname = getFileBaseName(seg)
3262 <                seg_out = tmp_dir + fsegname + ".to_spliced.bwtout"
3263 <                if use_zpacker: seg_out += ".z"
3262 >                seg_out = tmp_dir + fsegname + ".to_spliced"
3263 >                #if use_zpacker: seg_out += ".z"
3264                  extra_output = "(%d/%d)" % (i+1, len(maps[ri].segs))
3265                  (seg_map, unmapped) = bowtie(params,
3266                                               tmp_dir + junc_idx_prefix,
# Line 3263 | Line 3270
3270                                               seg_out,
3271                                               None,
3272                                               extra_output,
3273 <                                             "seg_against_junctions")
3273 >                                             _segs_vs_J)
3274                  spliced_seg_maps.append(seg_map)
3275                  i += 1
3276  
# Line 3320 | Line 3327
3327  
3328                  print >> run_log, " ".join(merge_cmd)
3329                  ret = subprocess.call( merge_cmd,
3330 <                                   stderr=open(logging_dir + "sam_merge.log", "w") )
3330 >                                   stderr=open(logging_dir + "bam_merge.log", "w") )
3331                  if ret != 0:
3332                      die(fail_str+"Error executing: "+" ".join(merge_cmd))
3333              except OSError, o:
# Line 3513 | Line 3520
3520                 if not reads_list:
3521                    continue
3522                 params.preflt_data[ri].multihit_reads = tmp_dir + sides[ri]+"_multimapped.fq"
3523 <               side_imap = tmp_dir + sides[ri]+"_im.bwtout"
3524 <               if use_zpacker:
3518 <                   side_imap+=".z"
3523 >               side_imap = tmp_dir + sides[ri]+"_im"
3524 >               #if use_zpacker: side_imap+=".z"
3525                 side_ium = tmp_dir + sides[ri]+"_ium.fq"
3526                 if use_BWT_FIFO:
3527                    side_ium += ".z"
# Line 3525 | Line 3531
3531                       params.read_params.reads_format,
3532                       params.genome_read_mismatches,
3533                       side_imap, side_ium,
3534 <                     "", "read_against_genome", params.preflt_data[ri].multihit_reads)
3534 >                     "", _reads_vs_G, params.preflt_data[ri].multihit_reads)
3535                 params.preflt_data[ri].mappings = bwt[0] # initial mappings
3536                 params.preflt_data[ri].unmapped_reads = bwt[1] # IUM reads
3537  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines