ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/tophat
(Generate patch)
# Line 6 | Line 6
6  
7   Created by Cole Trapnell on 2008-12-25.
8   Copyright (c) 2008 Cole Trapnell. All rights reserved.
9 + Updated and maintained by Daehwan Kim and Geo Pertea since Jul 2010.
10   """
11   import sys
12   try:
# Line 22 | Line 23
23   import re
24   import signal
25   from datetime import datetime, date, time
26 + from shutil import copy
27  
28   use_message = '''
29   TopHat maps short sequences from spliced transcripts to whole genomes.
# Line 33 | Line 35
35   Options:
36      -v/--version
37      -o/--output-dir                <string>    [ default: ./tophat_out     ]
38 +    --bowtie1                                  [ default: bowtie2          ]
39      -a/--min-anchor                <int>       [ default: 8                ]
40      -m/--splice-mismatches         <0-2>       [ default: 0                ]
41      -i/--min-intron-length         <int>       [ default: 50               ]
# Line 43 | Line 46
46      -M/--prefilter-multihits                   ( for -G/--GTF option, enable
47                                                   an initial bowtie search
48                                                   against the genome )
46    -F/--min-isoform-fraction      <float>     [ default: 0.15             ]
49      --max-insertion-length         <int>       [ default: 3                ]
50      --max-deletion-length          <int>       [ default: 3                ]
51      --solexa-quals
# Line 55 | Line 57
57      --color-out
58      --library-type                 <string>    (fr-unstranded, fr-firststrand,
59                                                  fr-secondstrand)
60 <    -p/--num-threads               <int>       [ default: 1                ]
61 <    -G/--GTF                       <filename>
62 <    --transcriptome-index          <bwtidx>    [pre-built bowtie index for the
63 <                                                transcripts in the -G/--GTF file]
62 <    -T/--transcriptome-only                    [map only to the transcriptome]
60 >    -p/--num-threads               <int>       [ default: 1                   ]
61 >    -G/--GTF                       <filename>  (GTF/GFF with known transcripts)
62 >    --transcriptome-index          <bwtidx>    (transcriptome bowtie index)
63 >    -T/--transcriptome-only                    (map only to the transcriptome)
64      -j/--raw-juncs                 <filename>
65      --insertions                   <filename>
66      --deletions                    <filename>
67      -r/--mate-inner-dist           <int>
68 <    --mate-std-dev                 <int>       [ default: 20               ]
68 >    --mate-std-dev                 <int>       [ default: 20                  ]
69      --no-novel-juncs
70      --no-novel-indels
71 +    --fusion-search
72 +    --fusion-anchor-length         <int>       [ default: 20               ]
73 +    --fusion-min-dist              <int>       [ default: 10000000         ]
74 +    --fusion-read-mismatches       <int>       [ default: 2                ]
75 +    --fusion-multireads            <int>       [ default: 2                ]
76 +    --fusion-multipairs            <int>       [ default: 2                ]
77      --no-gtf-juncs
78      --no-coverage-search
79      --coverage-search
# Line 82 | Line 89
89                                                   files for color space reads]
90  
91   Advanced Options:
92 <    -N/--initial-read-mismatches   <int>       [ default: 2                ]
92 >    --genome-read-mismatches       <int>       [ default: 2                ]
93 >    -N/--read-mismatches           <int>       [ default: 2                ]
94      --segment-mismatches           <int>       [ default: 2                ]
95      --segment-length               <int>       [ default: 25               ]
96      --bowtie-n                                 [ default: bowtie -v        ]
# Line 93 | Line 101
101      --max-coverage-intron          <int>       [ default: 20000            ]
102      --min-segment-intron           <int>       [ default: 50               ]
103      --max-segment-intron           <int>       [ default: 500000           ]
104 <    --no-sort-bam                              [Output BAM is not coordinate-sorted]
105 <    --no-convert-bam                           [Do not convert to bam format.
104 >    --no-sort-bam                              (Output BAM is not coordinate-sorted)
105 >    --no-convert-bam                           (Do not convert to bam format.
106                                                  Output is <output_dir>accepted_hit.sam.
107 <                                                Implies --no-sort-bam.]
107 >                                                Implies --no-sort-bam)
108 >
109 > Bowtie2 related options:
110 >  Preset options in --end-to-end mode (local alignment is not used in TopHat2)
111 >    --b2-very-fast
112 >    --b2-fast
113 >    --b2-sensitive
114 >    --b2-very-sensitive
115 >
116 >  Alignment options
117 >    --b2-N                         <int>       [ default: 0                ]
118 >    --b2-L                         <int>       [ default: 20               ]
119 >    --b2-i                         <func>      [ default: S,1,1.25         ]
120 >    --b2-n-ceil                    <func>      [ default: L,0,0.15         ]
121 >    --b2-gbar                      <int>       [ default: 4                ]
122 >
123 >  Scoring options
124 >    --b2-mp                        <int>,<int> [ default: 6,2              ]
125 >    --b2-np                        <int>       [ default: 1                ]
126 >    --b2-rdg                       <int>,<int> [ default: 5,3              ]
127 >    --b2-rfg                       <int>,<int> [ default: 5,3              ]
128 >    --b2-score-min                 <func>      [ default: L,-0.6,-0.6      ]
129 >
130 >  Effort options
131 >    --b2-D                         <int>       [ default: 15               ]
132 >    --b2-R                         <int>       [ default: 2                ]
133 >
134  
135   SAM Header Options (for embedding sequencing run metadata in output):
136      --rg-id                        <string>    (read group ID)
# Line 109 | Line 143
143      --rg-platform                  <string>    (Sequencing platform descriptor)
144   '''
145  
146 + #    -F/--min-isoform-fraction      <float>     [ default: 0.15             ]
147 + #    (obsolete)
148 +
149   class Usage(Exception):
150      def __init__(self, msg):
151          self.msg = msg
# Line 144 | Line 181
181   # and option validation via the member function check()
182  
183  
184 < class BowtieOutputFiles():
184 > class BowtieOutputFiles:
185      def __init__(self,
186            mappings=None,
187            unmapped_reads=None,
# Line 203 | Line 240
240  
241      class SystemParams:
242          def __init__(self,
243 <                     num_cpus,
243 >                     num_threads,
244                       keep_tmp):
245 <            self.num_cpus = num_cpus
245 >            self.num_threads = num_threads
246              self.keep_tmp = keep_tmp
247              self.zipper = "gzip"
248              self.zipper_opts= []
# Line 215 | Line 252
252              global use_BWT_FIFO
253              for option, value in opts:
254                  if option in ("-p", "--num-threads"):
255 <                    self.num_cpus = int(value)
255 >                    self.num_threads = int(value)
256                  elif option == "--keep-tmp":
257                      self.keep_tmp = True
258                  elif option in ("-z","--zpacker"):
# Line 228 | Line 265
265                      use_BWT_FIFO=True
266              if self.zipper:
267                  use_zpacker=True
268 <                if self.num_cpus>1 and not self.zipper_opts:
268 >                if self.num_threads>1 and not self.zipper_opts:
269                      if self.zipper.endswith('pbzip2') or self.zipper.endswith('pigz'):
270 <                         self.zipper_opts.append('-p'+str(self.num_cpus))
270 >                         self.zipper_opts.append('-p'+str(self.num_threads))
271              else:
272                  use_zpacker=False
273                  if use_BWT_FIFO: use_BWT_FIFO=False
# Line 238 | Line 275
275              cmdline=[]
276              if self.zipper:
277                   cmdline.extend(['-z',self.zipper])
278 <            if self.num_cpus>1:
279 <                 cmdline.extend(['-p'+str(self.num_cpus)])
278 >            if self.num_threads>1:
279 >                 cmdline.extend(['-p'+str(self.num_threads)])
280              return cmdline
281  
282          def check(self):
283 <            if self.num_cpus<1 :
283 >            if self.num_threads<1 :
284                   die("Error: arg to --num-threads must be greater than 0")
285              if self.zipper:
286                  xzip=which(self.zipper)
# Line 336 | Line 373
373          def check(self):
374              if self.seed_length and self.seed_length < 20:
375                  die("Error: arg to --seed-length must be at least 20")
376 +
377              if self.mate_inner_dist_std_dev != None and self.mate_inner_dist_std_dev < 0:
378                  die("Error: arg to --mate-std-dev must at least 0")
379              if (not self.read_group_id and self.sample_id) or (self.read_group_id and not self.sample_id):
# Line 408 | Line 446
446                      self.convert_bam = False
447                      self.sort_bam = False
448  
449 +    class Bowtie2Params:
450 +        def __init__(self):
451 +            self.very_fast = False
452 +            self.fast = False
453 +            self.sensitive = False
454 +            self.very_sensitive = False
455 +
456 +            self.N = 0
457 +            self.L = 20
458 +            self.i = "S,1,1.25"
459 +            self.n_ceil = "L,0,0.15"
460 +            self.gbar = 4
461 +
462 +            self.mp = "6,2"
463 +            self.np = 1
464 +            self.rdg = "5,3"
465 +            self.rfg = "5,3"
466 +            self.score_min = "L,-0.6,-0.6"
467 +
468 +            self.D = 15
469 +            self.R = 2
470 +
471 +        def parse_options(self, opts):
472 +            for option, value in opts:
473 +                if option == "--b2-very-fast":
474 +                    self.very_fast = True
475 +                if option == "--b2-fast":
476 +                    self.fast = True
477 +                if option == "--b2-sensitive":
478 +                    self.sensitive = True
479 +                if option == "--b2-very-sensitive":
480 +                    self.very_sensitive = True
481 +
482 +                if option == "--b2-N":
483 +                    self.N = int(value)
484 +                if option == "--b2-L":
485 +                    self.L = 20
486 +                if option == "--b2-i":
487 +                    self.i = value
488 +                if option == "--b2-n-ceil":
489 +                    self.n_ceil = value
490 +                if option == "--b2-gbar":
491 +                    self.gbar = 4
492 +
493 +                if option == "--b2-mp":
494 +                    self.mp = value
495 +                if option == "--b2-np":
496 +                    self.np = int(value)
497 +                if option == "--b2-rdg":
498 +                    self.rdg = value
499 +                if option == "--b2-rfg":
500 +                    self.rfg = value
501 +                if option == "--b2-score-min":
502 +                    self.score_min = value
503 +
504 +                if option == "--b2-D":
505 +                    self.D = int(value)
506 +                if option == "--b2-R":
507 +                    self.R = int(value)
508 +
509 +        def check(self):
510 +            more_than_once = False
511 +            if self.very_fast:
512 +                if self.fast or self.sensitive or self.very_sensitive:
513 +                    more_than_once = True
514 +            else:
515 +                if self.fast:
516 +                    if self.sensitive or self.very_sensitive:
517 +                        more_than_once = True
518 +                else:
519 +                    if self.sensitive and self.very_sensitive:
520 +                        more_than_once = True
521 +
522 +            if more_than_once:
523 +                die("Error: use only one of --b2-very-fast, --b2-fast, --b2-sensitive, --b2-very-sensitive")
524 +
525 +            if not self.N in [0, 1]:
526 +                die("Error: arg to --b2-N must be either 0 or 1")
527 +
528 +            function_re = r'^[CLSG],-?\d+(\.\d+)?,-?\d+(\.\d+)?$'
529 +            function_match = re.search(function_re, self.i)
530 +
531 +            if not function_match:
532 +                die("Error: arg to --b2-i must be <func> (e.g. --b2-i S,1,1.25)")
533 +
534 +            function_match = re.search(function_re, self.n_ceil)
535 +            if not function_match:
536 +                die("Error: arg to --b2-n-ceil must be <func> (e.g. --b2-n-ceil L,0,0.15)")
537 +
538 +            function_match = re.search(function_re, self.score_min)
539 +            if not function_match:
540 +                die("Error: arg to --b2-score-min must be <func> (e.g. --b2-score-min L,-0.6,-0.6)")
541 +
542 +            pair_re = r'^\d+,\d+$'
543 +            pair_match = re.search(pair_re, self.mp)
544 +            if not pair_match:
545 +                die("Error: arg to --b2-mp must be <int>,<int> (e.g. --b2-mp 6,2)")
546 +
547 +            pair_match = re.search(pair_re, self.rdg)
548 +            if not pair_match:
549 +                die("Error: arg to --b2-rdg must be <int>,<int> (e.g. --b2-mp 5,3)")
550 +
551 +            pair_match = re.search(pair_re, self.rfg)
552 +            if not pair_match:
553 +                die("Error: arg to --b2-rfg must be <int>,<int> (e.g. --b2-mp 5,3)")
554 +
555  
556      def __init__(self):
557          self.splice_constraints = self.SpliceConstraints(8,     # min_anchor
# Line 438 | Line 582
582                                             None,                # run date
583                                             None)                # sequencing platform
584  
585 <        self.system_params = self.SystemParams(1,               # bowtie_threads (num_cpus)
585 >        self.system_params = self.SystemParams(1,               # bowtie_threads (num_threads)
586                                                 False)           # keep_tmp
587  
588          self.search_params = self.SearchParams(100,             # min_closure_exon_length
# Line 450 | Line 594
594                                                 500000)          # max_segment_intron_length
595  
596          self.report_params = self.ReportParams()
597 +
598 +        self.bowtie2_params = self.Bowtie2Params()
599 +
600 +        self.bowtie2 = True
601          self.gff_annotation = None
602          self.transcriptome_only = False
603          self.transcriptome_index = None
604 +        self.transcriptome_outdir = None
605          self.raw_junctions = None
606          self.find_novel_juncs = True
607          self.find_novel_indels = True
608 +        self.find_novel_fusions = True
609          self.find_GFF_juncs = True
460        #self.skip_check_reads = False
610          self.max_hits = 20
611          self.t_max_hits = 60
612          self.prefilter_multi = False
613 <        self.initial_read_mismatches = 2
613 >        self.genome_read_mismatches = 2
614 >        self.max_read_mismatches = 2
615          self.t_mismatches = 1
616          self.segment_length = 25
617          self.segment_mismatches = 2
# Line 474 | Line 624
624          self.closure_search = False
625          self.microexon_search = False
626          self.butterfly_search = None
627 +        self.fusion_search = False
628 +        self.fusion_anchor_length = 20
629 +        self.fusion_min_dist = 10000000
630 +        self.fusion_read_mismatches = 2
631 +        self.fusion_multireads = 2
632 +        self.fusion_multipairs = 2
633  
634      def check(self):
635          self.splice_constraints.check()
636          self.read_params.check()
637          self.system_params.check()
638 <
638 >        self.bowtie2_params.check()
639          if self.segment_length < 10:
640              die("Error: arg to --segment-length must at least 10")
641          if self.segment_mismatches < 0 or self.segment_mismatches > 3:
642              die("Error: arg to --segment-mismatches must in [0, 3]")
643 +        if self.t_mismatches < 0 or self.t_mismatches > 3:
644 +            die("Error: arg to -n or --transcriptome-mismatches must in [0, 3]")
645  
646 <        if self.read_params.color and self.butterfly_search:
647 <            die("Error: butterfly-search in colorspace is not yet supported")
646 >        if self.read_params.color:
647 >            if self.bowtie2:
648 >                die("Error: bowtie2 in colorspace is not yet supported\nuse --bowtie1 option")
649 >            if self.fusion_search:
650 >                die("Error: fusion-search in colorspace is not yet supported")
651 >            if self.butterfly_search:
652 >                die("Error: butterfly-search in colorspace is not yet supported")
653  
654          library_types = ["fr-unstranded", "fr-firststrand", "fr-secondstrand"]
655  
# Line 531 | Line 694
694                 "--max-coverage-intron", str(self.search_params.max_coverage_intron_length),
695                 "--min-segment-intron", str(self.search_params.min_segment_intron_length),
696                 "--max-segment-intron", str(self.search_params.max_segment_intron_length),
697 <               "--max-mismatches", str(self.initial_read_mismatches),
697 >               "--max-mismatches", str(self.max_read_mismatches),
698                 "--sam-header", sam_header,
699                 "--max-insertion-length", str(self.max_insertion_length),
700                 "--max-deletion-length", str(self.max_deletion_length)]
701 +
702 +        if not self.bowtie2:
703 +            cmd.extend(["--bowtie1"])
704 +
705 +        if self.fusion_search:
706 +            cmd.extend(["--fusion-search",
707 +                        "--fusion-anchor-length", str(self.fusion_anchor_length),
708 +                        "--fusion-min-dist", str(self.fusion_min_dist),
709 +                        "--fusion-read-mismatches", str(self.fusion_read_mismatches),
710 +                        "--fusion-multireads", str(self.fusion_multireads),
711 +                        "--fusion-multipairs", str(self.fusion_multipairs)])
712 +
713          cmd.extend(self.system_params.cmd())
714  
715          if self.read_params.mate_inner_dist != None:
# Line 546 | Line 721
721                 cmd.extend(["--gtf-juncs", gtf_juncs])
722          if self.closure_search == False:
723              cmd.append("--no-closure-search")
724 <        if self.coverage_search == False:
724 >        if not self.coverage_search:
725              cmd.append("--no-coverage-search")
726 <        if self.microexon_search == False:
726 >        if not self.microexon_search:
727              cmd.append("--no-microexon-search")
728          if self.butterfly_search:
729              cmd.append("--butterfly-search")
# Line 579 | Line 754
754                                          ["version",
755                                           "help",
756                                           "output-dir=",
757 +                                         "bowtie1",
758                                           "solexa-quals",
759                                           "solexa1.3-quals",
760                                           "phred64-quals",
# Line 600 | Line 776
776                                           "transcriptome-index=",
777                                           "raw-juncs=",
778                                           "no-novel-juncs",
779 +                                         "allow-fusions",
780 +                                         "fusion-search",
781 +                                         "fusion-anchor-length=",
782 +                                         "fusion-min-dist=",
783 +                                         "fusion-read-mismatches=",
784 +                                         "fusion-multireads=",
785 +                                         "fusion-multipairs=",
786                                           "no-novel-indels",
787                                           "no-gtf-juncs",
605                                         # "skip-check-reads",
788                                           "mate-inner-dist=",
789                                           "mate-std-dev=",
790                                           "no-closure-search",
# Line 619 | Line 801
801                                           "min-segment-intron=",
802                                           "max-segment-intron=",
803                                           "seed-length=",
804 <                                         "initial-read-mismatches=",
804 >                                         "genome-read-mismatches=",
805                                           "read-mismatches=",
806 <                                         "max-mismatches=",
806 >                                         "max-read-mismatches=",
807                                           "transcriptome-mismatches=",
808                                           "segment-length=",
809                                           "segment-mismatches=",
# Line 645 | Line 827
827                                           "insertions=",
828                                           "deletions=",
829                                           "no-sort-bam",
830 <                                         "no-convert-bam"])
830 >                                         "no-convert-bam",
831 >                                         "b2-very-fast",
832 >                                         "b2-fast",
833 >                                         "b2-sensitive",
834 >                                         "b2-very-sensitive",
835 >                                         "b2-N=",
836 >                                         "b2-L=",
837 >                                         "b2-i=",
838 >                                         "b2-n-ceil=",
839 >                                         "b2-gbar=",
840 >                                         "b2-ma=",
841 >                                         "b2-mp=",
842 >                                         "b2-np=",
843 >                                         "b2-rdg=",
844 >                                         "b2-rfg=",
845 >                                         "b2-score-min=",
846 >                                         "b2-D=",
847 >                                         "b2-R="])
848          except getopt.error, msg:
849              raise Usage(msg)
850  
# Line 654 | Line 853
853          self.read_params.parse_options(opts)
854          self.search_params.parse_options(opts)
855          self.report_params.parse_options(opts)
856 +        self.bowtie2_params.parse_options(opts)
857          global use_BWT_FIFO
858          global use_BAM_Unmapped
859          if not self.read_params.color:
# Line 673 | Line 873
873                  sys.exit(0)
874              if option in ("-h", "--help"):
875                  raise Usage(use_message)
876 +            if option == "--bowtie1":
877 +                self.bowtie2 = False
878              if option in ("-g", "--max-multihits"):
879                  self.max_hits = int(value)
880              if option in ("-x", "--transcriptome-max-hits"):
# Line 681 | Line 883
883                  self.gff_annotation = value
884              if option in ("-T", "--transcriptome"):
885                  self.transcriptome_only = True
886 <            if option in ("--transcriptome-index"):
886 >            if option == "--transcriptome-index":
887                  self.transcriptome_index = value
888              if option in("-M", "--prefilter-multihits"):
889                  self.prefilter_multi = True
# Line 691 | Line 893
893                  self.find_novel_juncs = False
894              if option == "--no-novel-indels":
895                  self.find_novel_indels = False
896 +            if option == "--fusion-search":
897 +                self.fusion_search = True
898 +            if option == "--fusion-anchor-length":
899 +                self.fusion_anchor_length = int(value)
900 +            if option == "--fusion-min-dist":
901 +                self.fusion_min_dist = int(value)
902 +            if option == "--fusion-read-mismatches":
903 +                self.fusion_read_mismatches = int(value)
904 +            if option == "--fusion-multireads":
905 +                self.fusion_multireads = int(value)
906 +            if option == "--fusion-multipairs":
907 +                self.fusion_multipairs = int(value)
908              if option == "--no-gtf-juncs":
909                  self.find_GFF_juncs = False
910              #if option == "--skip-check-reads":
# Line 709 | Line 923
923                  self.butterfly_search = True
924              if option == "--no-butterfly-search":
925                  self.butterfly_search = False
926 <            if option in ("-N", "--initial-read-mismatches", "--read-mismatches", "--max-mismatches"):
927 <                self.initial_read_mismatches = int(value)
926 >            if option == "--genome-read-mismatches":
927 >                self.genome_read_mismatches = int(value)
928 >            if option in ("-N", "--read-mismatches", "--max-read-mismatches"):
929 >                self.max_read_mismatches = int(value)
930              if option in ("-n", "--transcriptome-mismatches"):
931                  self.t_mismatches = int(value)
932              if option == "--segment-length":
# Line 747 | Line 963
963          return args
964  
965  
966 + # check if a file exists and has non-zero (or minimum) size
967 + def fileExists(filepath, minfsize=1):
968 +  if os.path.exists(filepath) and os.path.getsize(filepath)>=minfsize:
969 +     return True
970 +  else:
971 +     return False
972 +
973 +
974   def getFileDir(filepath):
975 <   #if fullpath, returns path including the ending /
975 >   #if fullpath given, returns path including the ending /
976     fpath, fname=os.path.split(filepath)
977     if fpath: fpath+='/'
978     return fpath
# Line 767 | Line 991
991        else:
992           return fbase
993     else:
994 <     if len(fbase)>len(fext):
994 >     if len(fbase)>0:
995          return fbase
996       else:
997          return fname
# Line 832 | Line 1056
1056          if bowtie_idx_env_var == None:
1057              die("Error: Could not find Bowtie index files " + idx_prefix + ".*")
1058          idx_prefix = bowtie_idx_env_var + idx_prefix
1059 +
1060          idx_fwd_1 = idx_prefix + ".1.ebwt"
1061          idx_fwd_2 = idx_prefix + ".2.ebwt"
1062          idx_rev_1 = idx_prefix + ".rev.1.ebwt"
# Line 903 | Line 1128
1128  
1129   # Retrive a tuple containing the system's version of Bowtie.  Parsed from
1130   # `bowtie --version`
1131 < def get_bowtie_version():
1131 > def get_bowtie_version(is_bowtie2):
1132      try:
1133          # Launch Bowtie to capture its version info
1134          proc = subprocess.Popen([bowtie_path, "--version"],
1135                            stdout=subprocess.PIPE)
1136 +
1137          stdout_value = proc.communicate()[0]
1138 +
1139          bowtie_version = None
1140          bowtie_out = repr(stdout_value)
1141  
1142          # Find the version identifier
1143 <        version_str = "bowtie version "
1143 >        if is_bowtie2:
1144 >            version_str = "bowtie2-align version "
1145 >        else:
1146 >            version_str = "bowtie version "
1147 >
1148          ver_str_idx = bowtie_out.find(version_str)
1149          if ver_str_idx != -1:
1150              nl = bowtie_out.find("\\n", ver_str_idx)
1151              version_val = bowtie_out[ver_str_idx + len(version_str):nl]
1152 <            bowtie_version = [int(x) for x in version_val.split('.')]
1152 >
1153 >            if is_bowtie2:
1154 >                ver_numbers = version_val.split('.')
1155 >                ver_numbers[2:] = ver_numbers[2].split('-')
1156 >                bowtie_version = [int(x) for x in ver_numbers[:3]] + [int(ver_numbers[3][4:])]
1157 >            else:
1158 >                bowtie_version = [int(x) for x in version_val.split('.')]
1159 >
1160          if len(bowtie_version) == 3:
1161              bowtie_version.append(0)
1162  
# Line 929 | Line 1167
1167             errmsg+="Error: bowtie not found on this system"
1168         die(errmsg)
1169  
1170 < def get_index_sam_header(read_params, idx_prefix):
1170 > def get_index_sam_header(is_bowtie2, read_params, idx_prefix, sort = True):
1171      try:
1172          bowtie_sam_header_filename = tmp_dir+idx_prefix.split('/')[-1]+".bwt.samheader.sam"
1173          bowtie_sam_header_file = open(bowtie_sam_header_filename,"w")
1174  
1175 <        bowtie_header_cmd = [bowtie_path, "--sam"]
1175 >        bowtie_header_cmd = [bowtie_path]
1176 >        if not is_bowtie2:
1177 >            bowtie_header_cmd += ["--sam"]
1178 >
1179          if read_params.color:
1180              bowtie_header_cmd.append('-C')
1181 +
1182          bowtie_header_cmd.extend([idx_prefix, '/dev/null'])
1183          subprocess.call(bowtie_header_cmd,
1184                     stdout=bowtie_sam_header_file,
1185                     stderr=open('/dev/null'))
1186  
1187          bowtie_sam_header_file.close()
946        bowtie_sam_header_file = open(bowtie_sam_header_filename,"r")
1188  
1189 +        if not sort:
1190 +            return bowtie_sam_header_filename
1191 +
1192 +        bowtie_sam_header_file = open(bowtie_sam_header_filename, "r")
1193          sam_header_file = open(sam_header, "w")
1194  
1195          preamble = []
# Line 971 | Line 1216
1216          #for line in preamble:
1217          #    print >> sam_header_file, line
1218          print >> sam_header_file, "@HD\tVN:1.0\tSO:coordinate"
974
1219          if read_params.read_group_id and read_params.sample_id:
1220              rg_str = "@RG\tID:%s\tSM:%s" % (read_params.read_group_id,
1221                                              read_params.sample_id)
# Line 999 | Line 1243
1243  
1244          sam_header_file.close()
1245  
1246 +        copy(sam_header, bowtie_sam_header_filename)
1247          return sam_header
1248  
1249      except OSError, o:
# Line 1008 | Line 1253
1253         die(errmsg)
1254  
1255   # Make sure Bowtie is installed and is recent enough to be useful
1256 < def check_bowtie():
1257 <    th_log("Checking for Bowtie")
1256 > def check_bowtie(is_bowtie2):
1257 >    log_msg = "Checking for Bowtie"
1258 >    if is_bowtie2:
1259 >        log_msg += "2"
1260 >    th_log(log_msg)
1261 >
1262 >    bowtie_bin = "bowtie"
1263 >    if is_bowtie2:
1264 >        bowtie_bin += "2"
1265 >
1266      global bowtie_path
1267 <    bowtie_path=prog_path("bowtie")
1268 <    bowtie_version = get_bowtie_version()
1267 >    bowtie_path=prog_path(bowtie_bin)
1268 >    bowtie_version = get_bowtie_version(is_bowtie2)
1269      if bowtie_version == None:
1270          die("Error: Bowtie not found on this system")
1271 <    elif bowtie_version[0] < 1 and bowtie_version[1] < 12 and bowtie_version[2] < 3:
1272 <        die("Error: TopHat requires Bowtie 0.12.3 or later")
1273 <    print >> sys.stderr, "\tBowtie version:\t\t\t %s" % ".".join([str(x) for x in bowtie_version])
1271 >    elif is_bowtie2:
1272 >        if bowtie_version[0] < 3 and bowtie_version[1] < 1 and bowtie_version[2] < 1 and bowtie_version[3] < 5:
1273 >            die("Error: TopHat requires Bowtie 2.0.0-beta5 or later")
1274 >    elif not is_bowtie2:
1275 >        if bowtie_version[0] < 1 and bowtie_version[1] < 12 and bowtie_version[2] < 3:
1276 >            die("Error: TopHat requires Bowtie 0.12.3 or later")
1277 >    print >> sys.stderr, "\t version:\t\t\t %s" % ".".join([str(x) for x in bowtie_version])
1278  
1279  
1280   # Retrive a tuple containing the system's version of samtools.  Parsed from
# Line 1056 | Line 1313
1313      print >> sys.stderr, "\tSamtools %s" % samtools_version_str
1314  
1315  
1059
1316   class FastxReader:
1317    def __init__(self, i_file, is_color=0, fname=''):
1318      self.bufline=None
# Line 1101 | Line 1357
1357      try:
1358        if line[0] != "@":
1359            raise ValueError("Records in Fastq files should start with '@' character")
1360 +
1361        seqid = line[1:].rstrip()
1362        seqstr = fline().rstrip()
1363  
# Line 1130 | Line 1387
1387            if qstrlen + self.isColor >= seq_len :
1388                 break # qv string has reached the length of seq string
1389            #loop until qv has the same length as seq
1390 +
1391        if self.isColor:
1392             # and qstrlen==seq_len :
1393             if qstrlen==seq_len:
# Line 1221 | Line 1479
1479      fhandle.write(">" + seq_id + "\n")
1480      for i in xrange(len(seq) / line_len + 1):
1481          start = i * line_len
1482 <        end = (i+1) * line_len if (i+1) * line_len < len(seq) else len(seq)
1482 >        #end = (i+1) * line_len if (i+1) * line_len < len(seq) else len(seq)
1483 >        if (i+1) * line_len < len(seq):
1484 >             end = (i+1) * line_len
1485 >        else:
1486 >             end = len(seq)
1487          fhandle.write( seq[ start:end ] + "\n")
1488  
1489   class ZReader:
# Line 1233 | Line 1495
1495          pipecmd=[]
1496          s=filename.lower()
1497          if s.endswith(".bam"):
1498 <           pipecmd=["bam2fastx", "-"]
1498 >           # daehwan - want to ouptut mapped reads as well
1499 >           pipecmd=[prog_path("bam2fastx"), "--all", "-"]
1500          else:
1501            if guess:
1502               if s.endswith(".z") or s.endswith(".gz") or s.endswith(".gzip"):
# Line 1254 | Line 1517
1517               if use_zpacker and filename.endswith(".z"):
1518                  pipecmd=[sysparams.zipper]
1519                  pipecmd.extend(sysparams.zipper_opts)
1520 +
1521            if pipecmd:
1522               pipecmd+=['-cd']
1523          if pipecmd:
# Line 1402 | Line 1666
1666   # Format a DateTime as a pretty string.
1667   # FIXME: Currently doesn't support days!
1668   def formatTD(td):
1669 <  hours = td.seconds // 3600
1670 <  minutes = (td.seconds % 3600) // 60
1671 <  seconds = td.seconds % 60
1672 <  return '%02d:%02d:%02d' % (hours, minutes, seconds)
1669 >    days = td.days
1670 >    hours = td.seconds // 3600
1671 >    minutes = (td.seconds % 3600) // 60
1672 >    seconds = td.seconds % 60
1673 >
1674 >    if days > 0:
1675 >        return '%d days %02d:%02d:%02d' % (days, hours, minutes, seconds)
1676 >    else:
1677 >        return '%02d:%02d:%02d' % (hours, minutes, seconds)
1678  
1679   class PrepReadsInfo:
1680      def __init__(self, fname, side):
# Line 1424 | Line 1693
1693                 print >> sys.stderr, "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"
1694  
1695  
1696 < def prep_reads_cmd(params, reads_list, quals_list=None, aux_file=None, filter_reads=None, hits_to_filter=None):
1696 > def prep_reads_cmd(params, reads_list, quals_list=None, aux_file=None, index_file=None, filter_reads=None, hits_to_filter=None):
1697    #generate a prep_reads cmd arguments
1698    filter_cmd = [prog_path("prep_reads")]
1699    filter_cmd.extend(params.cmd())
# Line 1437 | Line 1706
1706      filter_cmd += ["--flt-hits="+hits_to_filter]
1707    if aux_file:
1708      filter_cmd += ["--aux-outfile="+aux_file]
1709 +  if index_file:
1710 +      filter_cmd += ["--index-outfile="+index_file]
1711    if filter_reads:
1712      filter_cmd += ["--flt-reads="+filter_reads]
1713    filter_cmd.append(reads_list)
# Line 1451 | Line 1722
1722   #--> returns the tuple (internal_read_library_file, read_info)
1723   def prep_reads(params, reads_list, quals_list, side="left", prefilter_reads=None):
1724      reads_suffix = ".fq"
1725 <    if use_zpacker: reads_suffix += ".z"
1725 >
1726 >    # for parallelization, we don't compress the read files
1727 >    do_use_zpacker = use_zpacker
1728 >    if params.system_params.num_threads > 1:
1729 >        do_use_zpacker = False
1730 >
1731 >    if do_use_zpacker: reads_suffix += ".z"
1732 >
1733      output_name=side+"_kept_reads"
1734      kept_reads_filename = tmp_dir + output_name + reads_suffix
1735  
# Line 1462 | Line 1740
1740      filter_log = open(log_fname,"w")
1741  
1742      info_file=output_dir+output_name+".info"
1743 <    filter_cmd=prep_reads_cmd(params, reads_list, quals_list, info_file, prefilter_reads)
1743 >    index_file=kept_reads_filename+".index"
1744 >
1745 >    filter_cmd=prep_reads_cmd(params, reads_list, quals_list, info_file, index_file, prefilter_reads)
1746      shell_cmd = ' '.join(filter_cmd)
1747      #finally, add the compression pipe
1748      zip_cmd=[]
1749 <    if use_zpacker:
1749 >    if do_use_zpacker:
1750         zip_cmd=[ params.system_params.zipper ]
1751         zip_cmd.extend(params.system_params.zipper_opts)
1752         zip_cmd.extend(['-c','-'])
# Line 1475 | Line 1755
1755      retcode=0
1756      try:
1757          print >> run_log, shell_cmd
1758 <        if use_zpacker:
1758 >        if do_use_zpacker:
1759              filter_proc = subprocess.Popen(filter_cmd,
1760                                    stdout=subprocess.PIPE,
1761                                    stderr=filter_log)
# Line 1498 | Line 1778
1778      except OSError, o:
1779          errmsg=fail_str+str(o)
1780          die(errmsg+"\n"+log_tail(log_fname))
1501    kept_reads.close()
1781  
1782 +    kept_reads.close()
1783      return kept_reads_filename, PrepReadsInfo(info_file, side)
1784  
1785  
# Line 1512 | Line 1792
1792             mapped_reads,
1793             unmapped_reads,
1794             extra_output = "",
1795 <           t_mapping = False,
1795 >           mapping_type = "read_against_genome",
1796             multihits_out = None): #only --prefilter-multihits should activate this parameter for the initial prefilter search
1797      start_time = datetime.now()
1798      bwt_idx_name = bwt_idx_prefix.split('/')[-1]
1799      reads_file=reads_list[0]
1800      readfile_basename=getFileBaseName(reads_file)
1801 +
1802 +    g_mapping, t_mapping, seg_mapping = False, False, False
1803 +    if mapping_type == "read_against_transcriptome":
1804 +        t_mapping = True
1805 +    elif mapping_type == "seg_against_genome" or mapping_type == "seg_against_junctions":
1806 +        seg_mapping = True
1807 +    else:
1808 +        g_mapping = True
1809 +
1810 +    bowtie_str = "Bowtie"
1811 +    if params.bowtie2:
1812 +        bowtie_str += "2"
1813 +
1814 +    if seg_mapping:
1815 +        # we multiply max_hits by 2 when mapping segments
1816 +        # this multiplication factor (2) is also assumed in segment_juncs and long_spanning_reads,
1817 +        # which means "2" should not be changed without modifying those sub-programs.
1818 +        params.max_hits *= 2
1819 +        if not params.bowtie2:
1820 +            backup_bowtie_alignment_option = params.bowtie_alignment_option
1821 +            params.bowtie_alignment_option = "-v"
1822 +
1823      if t_mapping:
1824 <       print >> sys.stderr, "[%s] Mapping %s against transcriptome %s with Bowtie %s" % (start_time.strftime("%c"),
1825 <                         readfile_basename, bwt_idx_name, extra_output)
1824 >       print >> sys.stderr, "[%s] Mapping %s against transcriptome %s with %s %s" % (start_time.strftime("%c"),
1825 >                         readfile_basename, bwt_idx_name, bowtie_str, extra_output)
1826      else:
1827 <       print >> sys.stderr, "[%s] Mapping %s against %s with Bowtie %s" % (start_time.strftime("%c"),
1828 <                         readfile_basename, bwt_idx_name, extra_output)
1827 >       print >> sys.stderr, "[%s] Mapping %s against %s with %s %s" % (start_time.strftime("%c"),
1828 >                         readfile_basename, bwt_idx_name, bowtie_str, extra_output)
1829      bwt_log = open(logging_dir + 'bowtie.'+readfile_basename+'.fixmap.log', "w")
1830      #bwt_mapped=mapped_reads
1831      unmapped_reads_out=None
# Line 1544 | Line 1846
1846      # Launch Bowtie
1847      try:
1848          bowtie_cmd = [bowtie_path]
1547
1849          if reads_format == "fastq":
1850              bowtie_cmd += ["-q"]
1851          elif reads_format == "fasta":
# Line 1557 | Line 1858
1858  
1859          unzip_cmd=None
1860          bam_input=False
1861 <        if reads_list.endswith('.bam'):
1861 >        if len(reads_list) > 0 and reads_list[0].endswith('.bam'):
1862             bam_input=True
1863             unzip_cmd=[ prog_path('bam2fastx') ]
1864             if reads_format:
1865 <              unzip_cmd.append("--"+reads_format)
1865 >              unzip_cmd.append("--" + reads_format)
1866             unzip_cmd+=[reads_list]
1867  
1868          if use_zpacker and (unzip_cmd is None):
1869            unzip_cmd=[ params.system_params.zipper ]
1870            unzip_cmd.extend(params.system_params.zipper_opts)
1871 +          unzip_cmd+=['-cd']
1872  
1873          fifo_pid=None
1874          if use_FIFO:
# Line 1585 | Line 1887
1887                   os._exit(os.EX_OK)
1888  
1889          fix_map_cmd = [prog_path('fix_map_ordering')]
1890 +
1891 +        if params.bowtie2:
1892 +            if t_mapping or g_mapping:
1893 +                max_penalty, min_penalty = params.bowtie2_params.mp.split(',')
1894 +                max_penalty, min_penalty = int(max_penalty), int(min_penalty)
1895 +                min_score = (max_penalty - 1) * num_mismatches
1896 +                fix_map_cmd += ["--bowtie2-min-score", str(min_score)]
1897 +            
1898          samzip_cmd=None
1899 <        if unmapped_reads:
1899 >
1900 >        # daehwan - just use bam format for multi-threading.
1901 >        if unmapped_reads or params.system_params.num_threads >= 1:
1902             #mapping on the genome
1903             #fix_map_ordering will write BAM file directly, and unmapped_reads as BAM file too
1904             mapped_reads += ".bam"
1905 <           if params.read_params.color:
1906 <              fix_map_cmd += ["--sam-header", sam_header, "-", mapped_reads]
1907 <           else:
1908 <              unmapped_reads_out = unmapped_reads+".bam"
1909 <              fix_map_cmd += ["--sam-header", sam_header, "-", mapped_reads, unmapped_reads_out]
1905 >
1906 >           temp_sam_header = tmp_dir + bwt_idx_prefix.split('/')[-1] + ".bwt.samheader.sam"
1907 >           fix_map_cmd += ["--sam-header", temp_sam_header, "-", mapped_reads]
1908 >           if not params.read_params.color and unmapped_reads:
1909 >               unmapped_reads_out = unmapped_reads+".bam"
1910 >               fix_map_cmd += [unmapped_reads_out]
1911          else:
1912             #mapping on segment_juncs (spliced potential junctions)
1913             # don't use SAM headers, could be slow -- better use compressed SAM directly
# Line 1609 | Line 1922
1922             else:
1923                fix_map_cmd += [mapped_reads]
1924  
1925 +        fix_map_cmd += ["--index-outfile", mapped_reads + ".index"]
1926 +
1927          max_hits = params.max_hits
1928          if t_mapping:
1929            max_hits = params.t_max_hits
1930 <        if params.bowtie_alignment_option == "-n" and num_mismatches>3:
1931 <           num_mismatches=3
1932 <        bowtie_cmd += [params.bowtie_alignment_option, str(num_mismatches),
1933 <                         "-p", str(params.system_params.num_cpus),
1934 <                         "-k", str(max_hits),
1935 <                         "-m", str(max_hits),
1936 <                         "-S", "--sam-nohead"] #always use headerless SAM file
1930 >        if num_mismatches > 3:
1931 >           num_mismatches = 3
1932 >
1933 >        if params.bowtie2:
1934 >            if seg_mapping:
1935 >                # since bowtie2 does not suppress reads that map to many places,
1936 >                # we suppress those in segment_juncs and long_spanning_reads.
1937 >                bowtie_cmd += ["-k", str(max_hits + 1)]
1938 >            else:
1939 >                bowtie_cmd += ["-k", str(max_hits)]
1940 >
1941 >            bowtie2_params = params.bowtie2_params
1942 >            if seg_mapping:
1943 >                # after intensive testing,
1944 >                # the following parameters seem to work faster than Bowtie1 and as sensitive as Bowtie1,
1945 >                # but a chance for further improvements remains.
1946 >                bowtie_cmd += ["-N", str(min(num_mismatches, 1))]
1947 >                bowtie_cmd += ["-i", "C,2,0"]
1948 >                bowtie_cmd += ["-L", "20"]
1949 >                # bowtie_cmd += ["-L", str(min(params.segment_length, 20))]
1950 >            else:
1951 >                bowtie2_preset = ""
1952 >                if bowtie2_params.very_fast:
1953 >                    bowtie2_preset = "--very-fast"
1954 >                elif bowtie2_params.fast:
1955 >                    bowtie2_preset = "--fast"
1956 >                elif bowtie2_params.sensitive:
1957 >                    bowtie2_preset = "--sensitive"
1958 >                elif bowtie2_params.very_sensitive:
1959 >                    bowtie2_preset = "--very-sensitive"
1960 >
1961 >                if bowtie2_preset != "":
1962 >                    bowtie_cmd += [bowtie2_preset]
1963 >                else:
1964 >                    bowtie_cmd += ["-D", str(bowtie2_params.D),
1965 >                                   "-R", str(bowtie2_params.R),
1966 >                                   "-N", str(bowtie2_params.N),
1967 >                                   "-L", str(bowtie2_params.L),
1968 >                                   "-i", bowtie2_params.i]
1969 >
1970 >                # "--n-ceil" is not correctly parsed in Bowtie2,
1971 >                # I  (daehwan) already talked to Ben who will fix the problem.
1972 >                bowtie_cmd += [# "--n-ceil", bowtie2_params.n_ceil,
1973 >                               "--gbar", str(bowtie2_params.gbar),
1974 >                               "--mp", bowtie2_params.mp,
1975 >                               "--np", str(bowtie2_params.np),
1976 >                               "--rdg", bowtie2_params.rdg,
1977 >                               "--rfg", bowtie2_params.rfg,
1978 >                               "--score-min", bowtie2_params.score_min]
1979 >
1980 >        else:
1981 >            bowtie_cmd += [params.bowtie_alignment_option, str(num_mismatches),
1982 >                           "-k", str(max_hits),
1983 >                           "-m", str(max_hits),
1984 >                           "-S"]
1985 >
1986 >        bowtie_cmd += ["-p", str(params.system_params.num_threads)]
1987 >
1988 >        if params.bowtie2: #always use headerless SAM file
1989 >            bowtie_cmd += ["--sam-no-hd"]
1990 >        else:
1991 >            bowtie_cmd += ["--sam-nohead"]
1992 >
1993 >        if not params.bowtie2:
1994 >            if multihits_out:
1995 >                bowtie_cmd += ["--max", multihits_out]
1996 >            else:
1997 >                bowtie_cmd += ["--max", "/dev/null"]
1998 >
1999 >        if params.bowtie2:
2000 >            bowtie_cmd += ["-x"]
2001  
1623        if multihits_out:
1624            bowtie_cmd += ["--max", multihits_out]
1625        else:
1626            bowtie_cmd += ["--max", "/dev/null"]
2002          bowtie_cmd += [ bwt_idx_prefix ]
2003          bowtie_proc=None
2004          shellcmd=""
# Line 1669 | Line 2044
2044          if samzip_cmd:
2045             shellcmd += "|"+ ' '.join(samzip_cmd)+" > " + mapped_reads
2046             fix_order_proc = subprocess.Popen(fix_map_cmd,
2047 <                                          stdin=bowtie_proc.stdout,
2048 <                                          stdout=subprocess.PIPE)
2047 >                                             stdin=bowtie_proc.stdout,
2048 >                                             stdout=subprocess.PIPE)
2049             bowtie_proc.stdout.close() # needed?
2050             samzip_proc = subprocess.Popen(samzip_cmd,
2051                                   preexec_fn=subprocess_setup,
# Line 1696 | Line 2071
2071                    os.kill(fifo_pid, signal.SIGTERM)
2072                  except:
2073                    pass
2074 +
2075      except OSError, o:
2076          die(fail_str+"Error: "+str(o))
2077  
# Line 1710 | Line 2086
2086            pass
2087      if multihits_out and not os.path.exists(multihits_out):
2088          open(multihits_out, "w").close()
1713    return (mapped_reads, unmapped_reads_out)
2089  
2090 < def bowtie_segment(params,
2091 <           bwt_idx_prefix,
2092 <           reads_file,
2093 <           reads_format,
2094 <           mapped_reads,
2095 <           unmapped_reads = None,
1721 <           extra_output = ""):
2090 >    if seg_mapping:
2091 >        params.max_hits /= 2
2092 >        if not params.bowtie2:
2093 >            params.bowtie_alignment_option = backup_bowtie_alignment_option
2094 >
2095 >    return (mapped_reads, unmapped_reads_out)
2096  
1723    backup_bowtie_alignment_option = params.bowtie_alignment_option
1724    params.bowtie_alignment_option = "-v"
1725    params.max_hits *= 2
1726
1727    result = bowtie(params,
1728                    bwt_idx_prefix,
1729                    [reads_file],
1730                    reads_format,
1731                    params.segment_mismatches,
1732                    mapped_reads,
1733                    unmapped_reads,
1734                    extra_output)
1735
1736    params.bowtie_alignment_option = backup_bowtie_alignment_option
1737    params.max_hits /= 2
1738    return result
2097  
2098   # Retrieve a .juncs file from a GFF file by calling the gtf_juncs executable
2099   def get_gtf_juncs(gff_annotation):
# Line 1759 | Line 2117
2117              return (False, gtf_juncs_out_name)
2118          elif retcode != 0:
2119              die(fail_str+"Error: GTF junction extraction failed with err ="+str(retcode))
2120 +
2121      # cvg_islands not found
2122      except OSError, o:
2123         errmsg=fail_str+str(o)+"\n"
# Line 1768 | Line 2127
2127      return (True, gtf_juncs_out_name)
2128  
2129   # Call bowtie-build on the FASTA file of sythetic splice junction sequences
2130 < def build_juncs_bwt_index(external_splice_prefix, color):
2130 > def build_juncs_bwt_index(is_bowtie2, external_splice_prefix, color):
2131      th_log("Indexing splices")
2132      bowtie_build_log = open(logging_dir + "bowtie_build.log", "w")
2133  
2134      #user_splices_out_prefix  = output_dir + "user_splices_idx"
2135  
2136 <    bowtie_build_cmd = [prog_path("bowtie-build")]
2136 >    if is_bowtie2:
2137 >        bowtie_build_cmd = [prog_path("bowtie2-build")]
2138 >    else:
2139 >        bowtie_build_cmd = [prog_path("bowtie-build")]
2140 >
2141      if color:
2142          bowtie_build_cmd += ["-C"]
2143  
# Line 1796 | Line 2159
2159  
2160   # Build a splice index from a .juncs file, suitable for use with specified read
2161   # (or read segment) lengths
2162 < def build_juncs_index(min_anchor_length,
2162 > def build_juncs_index(is_bowtie2,
2163 >                      min_anchor_length,
2164                        max_seg_len,
2165                        juncs_prefix,
2166                        external_juncs,
2167                        external_insertions,
2168                        external_deletions,
2169 +                      external_fusions,
2170                        reference_fasta,
2171                        color):
2172      th_log("Retrieving sequences for splices")
# Line 1809 | Line 2174
2174      juncs_file_list = ",".join(external_juncs)
2175      insertions_file_list = ",".join(external_insertions)
2176      deletions_file_list = ",".join(external_deletions)
2177 <
2177 >    fusions_file_list = ",".join(external_fusions)
2178  
2179      juncs_db_log = open(logging_dir + "juncs_db.log", "w")
2180  
# Line 1824 | Line 2189
2189                      juncs_file_list,
2190                      insertions_file_list,
2191                      deletions_file_list,
2192 +                    fusions_file_list,
2193                      reference_fasta]
2194      try:
2195          print >> run_log, " ".join(juncs_db_cmd)
# Line 1840 | Line 2206
2206             errmsg+="Error: juncs_db not found on this system"
2207         die(errmsg)
2208  
2209 <    external_splices_out_prefix = build_juncs_bwt_index(external_splices_out_prefix, color)
2209 >    external_splices_out_prefix = build_juncs_bwt_index(is_bowtie2, external_splices_out_prefix, color)
2210      return external_splices_out_prefix
2211  
2212 < def build_idx_from_fa(fasta_fname, out_dir):
2212 > def build_idx_from_fa(is_bowtie2, fasta_fname, out_dir, color):
2213      """ Build a bowtie index from a FASTA file.
2214  
2215      Arguments:
# Line 1854 | Line 2220
2220      - The path to the Bowtie index.
2221      """
2222      bwt_idx_path = out_dir + os.path.basename(fasta_fname).replace(".fa", "")
2223 <    bowtie_idx_cmd = [prog_path("bowtie-build"),
2224 <                      "-f",
2225 <                      fasta_fname,
2226 <                      bwt_idx_path]
2223 >
2224 >    if is_bowtie2:
2225 >        bowtie_idx_cmd = [prog_path("bowtie2-build")]
2226 >    else:
2227 >        bowtie_idx_cmd = [prog_path("bowtie-build")]
2228 >
2229 >    if color:
2230 >        bowtie_idx_cmd += ["-C"]
2231 >
2232 >    bowtie_idx_cmd += [fasta_fname,
2233 >                       bwt_idx_path]
2234      try:
2235          th_log("Building Bowtie index from " + os.path.basename(fasta_fname))
2236          print >> run_log, " ".join(bowtie_idx_cmd)
# Line 1901 | Line 2274
2274      print >> sam_file, "@PG\tID:TopHat\tVN:%s\tCL:%s" % (get_version(), run_cmd)
2275  
2276   # Write final TopHat output, via tophat_reports and wiggles
2277 < def compile_reports(params, sam_header_filename, mappings, readfiles, gff_annotation):
2277 > def compile_reports(params, sam_header_filename, ref_fasta, mappings, readfiles, gff_annotation):
2278      th_log("Reporting output tracks")
2279      left_maps, right_maps = mappings
2280      left_reads, right_reads = readfiles
# Line 1916 | Line 2289
2289      junctions = output_dir + "junctions.bed"
2290      insertions = output_dir + "insertions.bed"
2291      deletions = output_dir + "deletions.bed"
1919    coverage =  "coverage.wig"
2292      accepted_hits = output_dir + "accepted_hits"
2293      report_cmdpath = prog_path("tophat_reports")
2294 +    fusions = output_dir + "fusions.out"
2295      report_cmd = [report_cmdpath]
2296 +
2297 +    alignments_output_filename = tmp_dir + "accepted_hits"
2298 +
2299      report_cmd.extend(params.cmd())
2300      report_cmd.extend(["--samtools="+samtools_path])
2301 <    report_cmd.extend([junctions,
2301 >    report_cmd.extend([ref_fasta,
2302 >                       junctions,
2303                         insertions,
2304                         deletions,
2305 <                       "-",
2305 >                       fusions,
2306 >                       alignments_output_filename,
2307                         left_maps,
2308                         left_reads])
2309 +
2310      if len(right_maps) > 0 and right_reads:
2311          report_cmd.append(right_maps)
2312          report_cmd.append(right_reads)
2313 +
2314      # -- tophat_reports now produces (uncompressed) BAM stream,
1935    #    directly piped into samtools sort
2315      try:
2316          if params.report_params.convert_bam:
2317              if params.report_params.sort_bam:
2318 <                report_proc=subprocess.Popen(report_cmd,
2319 <                                             preexec_fn=subprocess_setup,
2320 <                                             stdout=subprocess.PIPE,
2321 <                                             stderr=report_log)
2322 <
2323 <                bamsort_cmd = [samtools_path, "sort", "-", accepted_hits]
2324 <                bamsort_proc= subprocess.Popen(bamsort_cmd,
2325 <                                               preexec_fn=subprocess_setup,
2326 <                                               stderr=open(logging_dir + "reports.samtools_sort.log", "w"),
2327 <                                               stdin=report_proc.stdout)
2328 <                report_proc.stdout.close()
2329 <                print >> run_log, " ".join(report_cmd)+" | " + " ".join(bamsort_cmd)
2330 <                bamsort_proc.communicate()
2331 <                retcode = report_proc.poll()
2332 <                if retcode:
2333 <                    die(fail_str+"Error running tophat_reports\n"+log_tail(log_fname))
2318 >                print >> run_log, " ".join(report_cmd)
2319 >                report_proc=subprocess.call(report_cmd,
2320 >                                            preexec_fn=subprocess_setup,
2321 >                                            stderr=report_log)
2322 >
2323 >                bam_parts = []
2324 >                for i in range(params.system_params.num_threads):
2325 >                    bam_part_filename = "%s%d.bam" % (alignments_output_filename, i)
2326 >                    if os.path.exists(bam_part_filename):
2327 >                        bam_parts.append(bam_part_filename)
2328 >                    else:
2329 >                        break
2330 >
2331 >                num_bam_parts = len(bam_parts)
2332 >                pids = [0 for i in range(num_bam_parts)]
2333 >                sorted_bam_parts = ["%s%d_sorted" % (alignments_output_filename, i) for i in range(num_bam_parts)]
2334 >                #left_um_parts = ["%s%s%d_sorted" % (alignments_output_filename, i) for i in range(num_bam_parts)]
2335 >                #right_um_parts = ["%s%d_sorted" % (alignments_output_filename, i) for i in range(num_bam_parts)]
2336 >                for i in range(num_bam_parts):
2337 >                    bamsort_cmd = [samtools_path,
2338 >                                   "sort",
2339 >                                   bam_parts[i],
2340 >                                   sorted_bam_parts[i]]
2341 >
2342 >                    sorted_bam_parts[i] += ".bam"
2343 >                    print >> run_log, " ".join(bamsort_cmd)
2344 >                    pid = os.fork()
2345 >                    if pid == 0:
2346 >                        subprocess.call(bamsort_cmd,
2347 >                                        stderr=open(logging_dir + "reports.samtools_sort.log%d" % i, "w"))
2348 >                        os._exit(os.EX_OK)
2349 >                    else:
2350 >                        pids[i] = pid
2351 >
2352 >                for i in range(len(pids)):
2353 >                    if pids[i] > 0:
2354 >                        os.waitpid(pids[i], 0)
2355 >                        pids[i] = 0
2356 >
2357 >                # for bam_part in bam_parts:
2358 >                #    os.remove(bam_part)
2359 >
2360 >                if num_bam_parts > 1:
2361 >                    bammerge_cmd = [samtools_path,
2362 >                                    "merge",
2363 >                                    "-h",
2364 >                                    sam_header_filename,
2365 >                                    "%s.bam" % accepted_hits]
2366 >                    bammerge_cmd += sorted_bam_parts
2367 >
2368 >                    print >> run_log, " ".join(bammerge_cmd)
2369 >                    subprocess.call(bammerge_cmd,
2370 >                                    stderr=open(logging_dir + "reports.samtools_merge.log", "w"))
2371 >
2372 >                    for sorted_bam_part in sorted_bam_parts:
2373 >                        os.remove(sorted_bam_part)
2374 >                else:
2375 >                    os.rename(sorted_bam_parts[0], output_dir + "accepted_hits.bam")
2376 >
2377              else:
2378                  print >> run_log, " ".join(report_cmd)
2379                  report_proc=subprocess.call(report_cmd,
# Line 1976 | Line 2398
2398                                    stderr=bam_to_sam_log)
2399              tmp_sam_file.close()
2400              os.remove(accepted_hits)
1979            #print >> run_log, "mv %s %s" % (tmp_sam, output_dir + "accepted_hits.sam")
1980            #os.rename(tmp_sam, output_dir + "accepted_hits.sam")
2401  
2402      except OSError, o:
2403            die(fail_str+"Error: "+str(o)+"\n"+log_tail(log_fname))
2404  
2405 < # FIXME: put wiggles back!
1986 < #        wig_cmd = [prog_path("wiggles"), output_dir + accepted_hits, output_dir + coverage]
1987 < #        print >> run_log, " ".join(wig_cmd)
1988 < #        subprocess.call(wig_cmd,
1989 < #                        stderr=open("/dev/null"))
1990 <    return (coverage, junctions)
2405 >    return junctions
2406  
2407  
2408   # Split up each read in a FASTQ file into multiple segments. Creates a FASTQ file
2409   # for each segment  This function needs to be fixed to support mixed read length
2410   # inputs
1996
2411   def open_output_files(prefix, num_files_prev, num_files, out_segf, extension, params):
2412         i = num_files_prev + 1
2413         while i <= num_files:
# Line 2099 | Line 2513
2513  
2514              # Bowtie's minimum read length here is 20bp, so if the last segment
2515              # is between 20 and segment_length bp long, go ahead and write it out
2516 <            if read_length % segment_length >= 20:
2516 >            if read_length % segment_length >= min(segment_length - 2, 20):
2517                  offsets.append(read_length)
2518                  tmp_num_segments += 1
2519              else:
# Line 2142 | Line 2556
2556                              ref_fasta):
2557      th_log("Searching for junctions via mate-pair closures")
2558  
2559 +
2560      #maps = [x for x in seg_maps if (os.path.exists(x) and os.path.getsize(x) > 0)]
2561      #if len(maps) == 0:
2562      #    return None
# Line 2150 | Line 2565
2565      juncs_out = ""
2566      if slash != -1:
2567          juncs_out += left_maps[0][:slash+1]
2568 +    fusions_out = juncs_out
2569 +
2570      juncs_out += "closure.juncs"
2571 +    fusions_out += "closure.fusions"
2572  
2573      juncs_log = open(logging_dir + "closure.log", "w")
2574      juncs_cmdpath=prog_path("closure_juncs")
# Line 2161 | Line 2579
2579  
2580      juncs_cmd.extend(params.cmd())
2581      juncs_cmd.extend([juncs_out,
2582 +                      fusions_out,
2583                        ref_fasta,
2584                        left_maps,
2585                        right_maps])
# Line 2199 | Line 2618
2618      juncs_out=out_path+"segment.juncs"
2619      insertions_out=out_path+"segment.insertions"
2620      deletions_out =out_path+"segment.deletions"
2621 +    fusions_out = out_path+"segment.fusions"
2622  
2623      left_maps = ','.join(left_seg_maps)
2624      log_fname = logging_dir + "segment_juncs.log"
# Line 2206 | Line 2626
2626      segj_cmd = [prog_path("segment_juncs")]
2627  
2628      segj_cmd.extend(params.cmd())
2209
2629      segj_cmd.extend(["--ium-reads", ",".join(unmapped_reads),
2630                        ref_fasta,
2631                        juncs_out,
2632                        insertions_out,
2633                        deletions_out,
2634 +                      fusions_out,
2635                        left_reads,
2636                        left_reads_map,
2637                        left_maps])
# Line 2227 | Line 2647
2647          # spanning_reads returned an error
2648          if retcode != 0:
2649             die(fail_str+"Error: segment-based junction search failed with err ="+str(retcode)+"\n"+log_tail(log_fname))
2650 +
2651      # cvg_islands not found
2652      except OSError, o:
2653          if o.errno == errno.ENOTDIR or o.errno == errno.ENOENT:
2654             print >> sys.stderr, fail_str, "Error: segment_juncs not found on this system"
2655          die(str(o))
2656  
2657 <    return [juncs_out, insertions_out, deletions_out]
2657 >    return [juncs_out, insertions_out, deletions_out, fusions_out]
2658  
2659   # Joins mapped segments into full-length read alignments via the executable
2660   # long_spanning_reads
# Line 2244 | Line 2665
2665                           possible_juncs,
2666                           possible_insertions,
2667                           possible_deletions,
2668 +                         possible_fusions,
2669                           contig_seg_maps,
2670                           spliced_seg_maps,
2671                           alignments_out_name):
# Line 2258 | Line 2680
2680      possible_juncs = ','.join(possible_juncs)
2681      possible_insertions = ",".join(possible_insertions)
2682      possible_deletions = ",".join(possible_deletions)
2683 +    possible_fusions = ",".join(possible_fusions)
2684 +
2685      log_fname=logging_dir + "long_spanning_reads"+rn+".log"
2686      align_log = open(log_fname, "w")
2687      align_cmd = [prog_path("long_spanning_reads")]
# Line 2265 | Line 2689
2689      align_cmd.extend(params.cmd())
2690  
2691      align_cmd.append(ref_fasta)
2692 <    align_cmd.extend([ reads,
2693 <                       possible_juncs,
2694 <                       possible_insertions,
2695 <                       possible_deletions,
2696 <                       contig_seg_maps])
2692 >    align_cmd.extend([reads,
2693 >                      possible_juncs,
2694 >                      possible_insertions,
2695 >                      possible_deletions,
2696 >                      possible_fusions,
2697 >                      alignments_out_name,
2698 >                      contig_seg_maps])
2699  
2700      if spliced_seg_maps:
2701          spliced_seg_maps = ','.join(spliced_seg_maps)
2702          align_cmd.append(spliced_seg_maps)
2703  
2704      try:
2705 <        print >> run_log, " ".join(align_cmd),">",alignments_out_name
2706 <        join_proc=subprocess.Popen(align_cmd,
2707 <                          preexec_fn=subprocess_setup,
2708 <                          stdout=open(alignments_out_name, "wb"),
2283 <                          stderr=align_log)
2284 <        join_proc.communicate()
2285 <        retcode=join_proc.poll()
2286 <        if retcode:
2287 <            die(fail_str+"Error at 'long_spanning_reads'\n"+log_tail(log_fname))
2705 >        print >> run_log, " ".join(align_cmd)
2706 >        join_proc=subprocess.call(align_cmd,
2707 >                                  stderr=align_log)
2708 >
2709      except OSError, o:
2710          die(fail_str+"Error: "+str(o))
2711  
2291
2712   # This class collects spliced and unspliced alignments for each of the
2713   # left and right read files provided by the user.
2714   class Maps:
# Line 2372 | Line 2792
2792      - `left_reads`: A list of reads.
2793      - `right_reads`: A list of reads (empty if single-end).
2794  
2375    Returns:
2376    TODO: Put something here
2795      """
2796      test_input_file(params.gff_annotation)
2797  
# Line 2382 | Line 2800
2800  
2801      gtf_name = getFileBaseName(params.gff_annotation)
2802      m2g_bwt_idx = None
2803 <    if params.transcriptome_index:
2803 >    t_out_dir = tmp_dir
2804 >    if params.transcriptome_index and not params.transcriptome_outdir:
2805         m2g_bwt_idx = params.transcriptome_index
2806 <       #th_log("Using pre-built transcriptome index..")
2806 >       th_log("Using pre-built transcriptome index..")
2807      else:
2808 <       th_log("Creating multi-FASTA from GTF and genome")
2809 <       m2g_ref_fasta = tmp_dir + gtf_name + ".fa"
2808 >       th_log("Creating transcriptome data files..")
2809 >       if params.transcriptome_outdir:
2810 >         t_out_dir=params.transcriptome_outdir+"/"
2811 >       m2g_ref_fasta = t_out_dir + gtf_name + ".fa"
2812         gtf_to_fasta(params, params.gff_annotation, ref_fasta, m2g_ref_fasta)
2813         # trans_hash = transcripts_to_fasta(transcripts,
2814         #                                   ref_fasta,
2815         #                                   m2g_ref_fasta)
2816 <       m2g_bwt_idx = build_idx_from_fa(m2g_ref_fasta, tmp_dir)
2816 >       m2g_bwt_idx = build_idx_from_fa(params.bowtie2, m2g_ref_fasta, t_out_dir, params.read_params.color)
2817      mapped_gtf_list = []
2818      unmapped_gtf_list = []
2819      # do the initial mapping in GTF coordinates
# Line 2415 | Line 2836
2836                                             params.t_mismatches,
2837                                             mapped_gtf_out,
2838                                             unmapped_gtf,
2839 <                                           "", True)
2839 >                                           "", "read_against_transcriptome")
2840          mapped_gtf_list.append(mapped_gtf_map)
2841          unmapped_gtf_list.append(unmapped)
2842  
# Line 2510 | Line 2931
2931      possible_insertions.extend(user_supplied_insertions)
2932      possible_deletions = []
2933      possible_deletions.extend(user_supplied_deletions)
2934 +    possible_fusions = []
2935  
2936      left_reads, right_reads = prepared_reads
2937  
2938      maps = [[], []] # maps[0] = left_reads mapping data, maps[1] = right_reads_mapping_data
2939 +
2940      #single_segments = False
2941  
2942      # Before anything, map the reads using Map2GTF (if using annotation)
# Line 2526 | Line 2949
2949  
2950          m2g_left_maps, m2g_right_maps = mapped_gtf_list
2951          m2g_maps = [m2g_left_maps, m2g_right_maps]
2952 <        if params.transcriptome_only or not os.path.exists(unmapped_gtf_list[0]):
2952 >        if params.transcriptome_only or not fileExists(unmapped_gtf_list[0]):
2953              # The case where the user doesn't want to map to anything other
2954              # than the transcriptome OR we have no unmapped reads
2955              maps[0] = [m2g_left_maps]
2956              if right_reads:
2957 <                maps[0] = [m2g_right_maps]
2957 >                maps[1] = [m2g_right_maps]
2958              return maps
2959          # Feed the unmapped reads into spliced_alignment()
2960          initial_reads = unmapped_gtf_list[:]
2538
2961          th_log("Resuming TopHat pipeline with unmapped reads")
2962  
2963      max_seg_len = segment_len #this is the ref seq span on either side of the junctions
# Line 2549 | Line 2971
2971         # the last segment is longer
2972         if num_segs>1: max_seg_len += (read_len % segment_len)
2973  
2552    # daehwan - check this out
2974      if num_segs <= 1:
2975 <         print >> sys.stderr, "Warning: you have only one segment per read\n\twe strongly recommend that you decrease --segment-length to about half the read length because TopHat will work better with multiple segments"
2975 >         print >> sys.stderr, "Warning: you have only one segment per read.\n\tIf the read length is greater than or equal to 45bp,\n\twe strongly recommend that you decrease --segment-length to about half the read length because TopHat will work better with multiple segments"
2976  
2977      # Using the num_segs value returned by check_reads(),
2978      # decide which junction discovery strategy to use
2979 <    if num_segs<3:
2979 >    if num_segs < 3:
2980         #if params.butterfly_search != False:
2981         #   params.butterfly_search = True
2982         if params.coverage_search != False:
2983 <               params.coverage_search = True
2983 >           params.coverage_search = True
2984         if num_segs == 1:
2985           segment_len = read_len
2986      else: #num_segs >= 3:
# Line 2572 | Line 2993
2993                 params.coverage_search = False
2994          if params.butterfly_search != True:
2995                  params.butterfly_search = False
2996 +
2997      # Perform the first part of the TopHat work flow on the left and right
2998      # reads of paired ends separately - we'll use the pairing information later
2999 +    have_left_IUM=False
3000      for ri in (0,1):
3001          reads=initial_reads[ri]
3002 <        if reads == None or os.path.getsize(reads)<25 :
3002 >        if reads == None or not fileExists(reads,25):
3003              continue
3004          fbasename=getFileBaseName(reads)
3005          unspliced_out = tmp_dir + fbasename + ".mapped"
# Line 2590 | Line 3013
3013          else:
3014          # Perform the initial Bowtie mapping of the full length reads
3015            (unspliced_sam, unmapped_reads) = bowtie(params,
3016 <                                             bwt_idx_prefix,
3017 <                                             [reads],
3018 <                                             "fastq",
3019 <                                             params.initial_read_mismatches,
3020 <                                             unspliced_out,
3021 <                                             unmapped_unspliced)
3022 <
3023 < #         unspliced_sam = tmp_dir+fbasename+'.unspl.bam'
3024 < #         join_mapped_segments(params,
3025 < #                              sam_header_filename,
3026 < #                              reads,
3027 < #                              ref_fasta,
3028 < #                              ["/dev/null"],
2606 < #                              ["/dev/null"],
2607 < #                              ["/dev/null"],
2608 < #                              [unspliced_map],
2609 < #                              [],
2610 < #                              unspliced_sam)
3016 >                                                   bwt_idx_prefix,
3017 >                                                   [reads],
3018 >                                                   "fastq",
3019 >                                                   params.genome_read_mismatches,
3020 >                                                   unspliced_out,
3021 >                                                   unmapped_unspliced,
3022 >                                                   "",
3023 >                                                   "read_against_genome")
3024 >
3025 >        # daehwan - check this out
3026 >        # params.bowtie2 = False
3027 >        # check_bowtie(params.bowtie2)
3028 >
3029          seg_maps = []
3030          unmapped_segs = []
3031          segs = []
3032 <        if num_segs > 1:
3032 >        have_IUM = fileExists(unmapped_reads,25)
3033 >        if ri==0 and have_IUM:
3034 >           have_left_IUM = True
3035 >        if num_segs > 1 and have_IUM:
3036              # split up the IUM reads into segments
3037              # unmapped_reads can be in BAM format
3038              read_segments = split_reads(unmapped_reads,
# Line 2630 | Line 3051
3051                  if use_BWT_FIFO:
3052                      unmapped_seg += ".z"
3053                  extra_output = "(%d/%d)" % (i+1, len(read_segments))
3054 <                (seg_map, unmapped) = bowtie_segment(params,
3055 <                                                     bwt_idx_prefix,
3056 <                                                     seg,
3057 <                                                     "fastq",
3058 <                                                     seg_out,
3059 <                                                     unmapped_seg,
3060 <                                                     extra_output)
3054 >                (seg_map, unmapped) = bowtie(params,
3055 >                                             bwt_idx_prefix,
3056 >                                             [seg],
3057 >                                             "fastq",
3058 >                                             params.segment_mismatches,
3059 >                                             seg_out,
3060 >                                             unmapped_seg,
3061 >                                             extra_output,
3062 >                                             "seg_against_genome")
3063                  seg_maps.append(seg_map)
3064                  unmapped_segs.append(unmapped)
3065                  segs.append(seg)
# Line 2655 | Line 3078
3078  
3079      # Unless the user asked not to discover new junctions, start that process
3080      # here
2658
3081      left_reads_map = maps[0].unspliced_sam
3082      left_seg_maps = maps[0].seg_maps
3083      unmapped_reads = maps[0].unmapped_segs
# Line 2667 | Line 3089
3089          right_reads_map = None
3090          right_seg_maps = None
3091  
3092 <    if params.find_novel_juncs: # or params.find_novel_indels:
2671 <
3092 >    if params.find_novel_juncs and have_left_IUM: # or params.find_novel_indels:
3093          # Call segment_juncs to infer a list of possible splice junctions from
3094          # the regions of the genome covered in the initial and segment maps
3095          #if params.find_novel_juncs:
# Line 2688 | Line 3109
3109          if os.path.getsize(juncs[0]) != 0:
3110                      possible_juncs.append(juncs[0])
3111          if params.find_novel_indels:
3112 <                if os.path.getsize(juncs[1]) != 0:
3113 <                    possible_insertions.append(juncs[1])
3114 <                if os.path.getsize(juncs[2]) != 0:
3115 <                    possible_deletions.append(juncs[2])
3112 >            if os.path.getsize(juncs[1]) != 0:
3113 >                possible_insertions.append(juncs[1])
3114 >            if os.path.getsize(juncs[2]) != 0:
3115 >                possible_deletions.append(juncs[2])
3116 >        if params.find_novel_fusions:
3117 >            if os.path.getsize(juncs[3]) != 0:
3118 >                possible_fusions.append(juncs[3])
3119          # Optionally, and for paired reads only, use a closure search to
3120          # discover addtional junctions
3121          if params.closure_search and left_reads and right_reads:
# Line 2702 | Line 3126
3126              if os.path.getsize(juncs[0]) != 0:
3127                  possible_juncs.extend(juncs)
3128  
3129 <    if len(possible_insertions) == 0 and len(possible_deletions) == 0 and len(possible_juncs) == 0:
3129 >    if len(possible_insertions) == 0 and len(possible_deletions) == 0 and len(possible_juncs) == 0 and len(possible_fusions) == 0:
3130          spliced_seg_maps = None
3131          junc_idx_prefix = None
3132      else:
# Line 2716 | Line 3140
3140      if len(possible_juncs) == 0:
3141          possible_juncs.append(os.devnull)
3142          print >> sys.stderr, "Warning: junction database is empty!"
3143 +    if len(possible_fusions) == 0:
3144 +        possible_fusions.append(os.devnull)
3145      if junc_idx_prefix:
3146 <        build_juncs_index(3,
3147 <                          #segment_len,
3148 <                          max_seg_len,
3149 <                          junc_idx_prefix,
3150 <                          possible_juncs,
3151 <                          possible_insertions,
3152 <                          possible_deletions,
3153 <                          ref_fasta,
3154 <                          params.read_params.color)
3146 >        juncs_bwt_idx = build_juncs_index(params.bowtie2,
3147 >                                          3,
3148 >                                          max_seg_len,
3149 >                                          junc_idx_prefix,
3150 >                                          possible_juncs,
3151 >                                          possible_insertions,
3152 >                                          possible_deletions,
3153 >                                          possible_fusions,
3154 >                                          ref_fasta,
3155 >                                          params.read_params.color)
3156 >
3157 >        juncs_bwt_samheader = get_index_sam_header(params.bowtie2, params.read_params, juncs_bwt_idx, False)
3158  
3159      # Now map read segments (or whole IUM reads, if num_segs == 1) to the splice
3160      # index with Bowtie
# Line 2747 | Line 3176
3176                  seg_out = tmp_dir + fsegname + ".to_spliced.bwtout"
3177                  if use_zpacker: seg_out += ".z"
3178                  extra_output = "(%d/%d)" % (i+1, len(maps[ri].segs))
3179 <                (seg_map, unmapped) = bowtie_segment(params,
3180 <                                                     tmp_dir + junc_idx_prefix,
3181 <                                                     seg,
3182 <                                                     "fastq",
3183 <                                                     seg_out,
3184 <                                                     None,
3185 <                                                     extra_output)
3179 >                (seg_map, unmapped) = bowtie(params,
3180 >                                             tmp_dir + junc_idx_prefix,
3181 >                                             [seg],
3182 >                                             "fastq",
3183 >                                             params.segment_mismatches,
3184 >                                             seg_out,
3185 >                                             None,
3186 >                                             extra_output,
3187 >                                             "seg_against_junctions")
3188                  spliced_seg_maps.append(seg_map)
3189                  i += 1
3190  
# Line 2773 | Line 3204
3204                               possible_juncs,
3205                               possible_insertions,
3206                               possible_deletions,
3207 +                             possible_fusions,
3208                               maps[ri].seg_maps,
3209                               spliced_seg_maps,
3210                               mapped_reads)
# Line 2785 | Line 3217
3217              # NOTE: We also should be able to address bug #134 here, by replacing
3218              # contiguous alignments that poke into an intron by a small amount by
3219              # the correct spliced alignment.
3220 +
3221              try:
3222 <                merge_cmd = [ prog_path("bam_merge"), merged_map, mapped_reads ]
3222 >                merge_cmd = [ prog_path("bam_merge"), "--index-outfile", merged_map + ".index",
3223 >                              merged_map ]
3224 >
3225                  if num_segs > 1:
3226                      merge_cmd += [ maps[ri].unspliced_sam ]
3227                  if m2g_map:
3228                      merge_cmd += [ m2g_map ]
3229 +
3230 +                if os.path.exists(mapped_reads):
3231 +                    merge_cmd += [ mapped_reads ]
3232 +                else:
3233 +                    for bam_i in range(0, params.system_params.num_threads):
3234 +                        temp_bam = mapped_reads[:-4] + str(bam_i) + ".bam"
3235 +                        if os.path.exists(temp_bam):
3236 +                            merge_cmd += [ temp_bam ]
3237 +                        else:
3238 +                            break
3239 +
3240                  print >> run_log, " ".join(merge_cmd)
3241                  ret = subprocess.call( merge_cmd,
3242                                     stderr=open(logging_dir + "sam_merge.log", "w") )
# Line 2809 | Line 3255
3255              maps[ri] = [merged_map]
3256          if not params.system_params.keep_tmp:
3257                  os.remove(unspl_samfile)
3258 <        #for ri
3258 >
3259      return maps
3260  
3261   def die(msg=None):
# Line 2844 | Line 3290
3290  
3291   # FIXME: this should get set during the make dist autotools phase of the build
3292   def get_version():
3293 <   return "1.4.0"
3293 >   return "2.0.0"
3294  
3295   def mlog(msg):
3296    print >> sys.stderr, "[DBGLOG]:"+msg
# Line 2852 | Line 3298
3298   def test_input_file(filename):
3299      try:
3300          test_file = open(filename, "r")
2855    # Bowtie not found
3301      except IOError, o:
3302          die("Error: Opening file %s" % filename)
3303      return
3304  
2860
3305   def main(argv=None):
3306      warnings.filterwarnings("ignore", "tmpnam is a potential security risk")
3307  
3308      # Initialize default parameter values
3309      params = TopHatParams()
2866
3310      try:
3311          if argv is None:
3312              argv = sys.argv
# Line 2902 | Line 3345
3345  
3346          # Validate all the input files, check all prereqs before committing
3347          # to the run
3348 +        if params.gff_annotation:
3349 +           if not os.path.exists(params.gff_annotation):
3350 +             die("Error: cannot find transcript file %s" % params.gff_annotation)
3351 +           if os.path.getsize(params.gff_annotation)<10:
3352 +             die("Error: invalid transcript file %s" % params.gff_annotation)
3353 +
3354          if params.transcriptome_index:
3355 <           check_bowtie_index(params.transcriptome_index)
3355 >           if params.gff_annotation:
3356 >               #gff file given, so transcriptome data will be written there
3357 >               gff_basename = getFileBaseName(params.gff_annotation)
3358 >               #just in case, check if it's not already there (-G/--GTF given again by mistake)
3359 >               tpath, tname = os.path.split(params.transcriptome_index)
3360 >               new_subdir=False
3361 >               if tpath in (".", "./") or not tpath :
3362 >                  if not os.path.exists(params.transcriptome_index):
3363 >                    os.makedirs(params.transcriptome_index)
3364 >                    new_subdir=True
3365 >               if new_subdir or (os.path.exists(params.transcriptome_index) and os.path.isdir(params.transcriptome_index)):
3366 >                   params.transcriptome_index = os.path.join(params.transcriptome_index, gff_basename)
3367 >               gff_out=params.transcriptome_index+".gff"
3368 >               if not (os.path.exists(gff_out) and os.path.getsize(gff_out)==os.path.getsize(params.gff_annotation)):
3369 >                  #generate the transcriptome data files
3370 >                  tpath, tname = os.path.split(params.transcriptome_index)
3371 >                  params.transcriptome_outdir=tpath
3372 >           t_gff=params.transcriptome_index+".gff"
3373 >           if params.transcriptome_outdir:
3374 >              #will create the transcriptome data files
3375 >              if not os.path.exists(params.transcriptome_outdir):
3376 >                os.makedirs(params.transcriptome_outdir)
3377 >              copy(params.gff_annotation, t_gff)
3378 >           else:
3379 >              #try to use existing transcriptome data files
3380 >
3381 >              if not (os.path.exists(t_gff) and os.path.getsize(t_gff)>10):
3382 >                  die("Error: GFF transcripts file not found or invalid (%s)" % t_gff)
3383 >              check_bowtie_index(params.transcriptome_index)
3384 >           params.gff_annotation=t_gff
3385 >           #end @ transcriptome_index given
3386 >
3387          (ref_fasta, ref_seq_dict) = check_index(bwt_idx_prefix)
3388  
3389 <        check_bowtie()
3389 >        check_bowtie(params.bowtie2)
3390          check_samtools()
3391          th_log("Generating SAM header for "+bwt_idx_prefix)
3392 <        sam_header_filename = get_index_sam_header(params.read_params, bwt_idx_prefix)
2913 <
3392 >        sam_header_filename = get_index_sam_header(params.bowtie2, params.read_params, bwt_idx_prefix)
3393          #if not params.skip_check_reads:
3394          reads_list = left_reads_list
3395          if right_reads_list:
# Line 2920 | Line 3399
3399          user_supplied_juncs = []
3400          user_supplied_insertions = []
3401          user_supplied_deletions = []
3402 +        user_supplied_fusions = []
3403          global gtf_juncs
3404          if params.gff_annotation and params.find_GFF_juncs:
3405              test_input_file(params.gff_annotation)
# Line 2964 | Line 3444
3444                 rdlist=reads_list.split(',')
3445                 bwt=bowtie(params, bwt_idx_prefix, rdlist,
3446                       params.read_params.reads_format,
3447 <                     params.initial_read_mismatches,
3447 >                     params.genome_read_mismatches,
3448                       side_imap, side_ium,
3449 <                     "", False, params.preflt_data[ri].multihit_reads)
3449 >                     "", "read_against_genome", params.preflt_data[ri].multihit_reads)
3450                 params.preflt_data[ri].mappings = bwt[0] # initial mappings
3451                 params.preflt_data[ri].unmapped_reads = bwt[1] # IUM reads
3452  
# Line 2998 | Line 3478
3478          # turn off integer-quals
3479          if params.read_params.integer_quals:
3480              params.read_params.integer_quals = False
3481 +
3482          input_reads = [left_kept_reads, right_kept_reads]
3483          mappings = spliced_alignment(params,
3484                                bwt_idx_prefix,
# Line 3012 | Line 3493
3493  
3494          compile_reports(params,
3495                          sam_header_filename,
3496 +                        ref_fasta,
3497                          mappings,
3498                          input_reads,
3499                          params.gff_annotation)
# Line 3033 | Line 3515
3515          finish_time = datetime.now()
3516          duration = finish_time - start_time
3517          print >> sys.stderr,"-----------------------------------------------"
3518 <        print >> sys.stderr, "Run complete [%s elapsed]" %  formatTD(duration)
3518 >        th_log("Run complete: %s elapsed" %  formatTD(duration))
3519  
3520      except Usage, err:
3521          print >> sys.stderr, sys.argv[0].split("/")[-1] + ": " + str(err.msg)

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines