Change polyA-based and canonical splice site transcript filtration

README.md

@@ -105,6 +105,9 @@ IsoQuant support all kinds of long RNA data:
 
 Reads must be provided in FASTQ or FASTA format (can be gzipped). If you have already aligned your reads to the reference genome, simply provide sorted and indexed BAM files.
 
+IsoQuant expect reads to contain polyA tails. For more reliable transcript model construction do not trim polyA tails.
+
+
 <a name="sec1.2"></a>
 ## Supported reference data
 
@@ -519,6 +522,26 @@ The main explanation that some aligner report a lot of false unspliced alignment
 for ONT reads.
 
 
+`--report_canonical`
+    Strategy for reporting novel transcripts based on canonical splice sites, should be one of:
+
+* `auto` - automatic selection based on the data type and model construction strategy (default); 
+* `only_canonical` - report novel transcripts, which contain only canonical splice sites;
+* `only_stranded` - report novel transcripts, for which the strand can be unambiguously derived using splice sites and 
+presence of a polyA tail, allowing some splice sites to be non-canonical
+* `all` -- report all transcript model regardless of their splice sites.
+
+
+`--polya_requirement` Strategy for using polyA tails during transcript model construction, should be one of:
+
+* `auto` - default behaviour: polyA tails are required if at least 70% of the reads have polyA tail; 
+polyA tails are always required for 1/2-exon transcripts when using ONT data (this is caused by elevated number of false 1/2-exonic alignments reported by minimap2); 
+* `never` - polyA tails are never required; use this option **at your own risk** as it may noticeably increase false discovery rate, especially for ONT data;
+* `always` - reported transcripts are always required to have polyA support in the reads.
+
+Note, that polyA tails are always required for reporting novel unspliced isoforms. 
+
+
 
 ### Hidden options
 <a name="hidden"></a>
@@ -528,9 +551,6 @@ We recommend _not_ to modify these options unless you are clearly aware of their
 `--no_secondary`
     Ignore secondary alignments.
 
-`--report_unstranded`
-    Report transcripts for which the strand cannot be detected using canonical splice sites.
-
 `--aligner`
     Force to use this alignment method, can be `starlong` or `minimap2`; `minimap2` is currently used as default. Make sure the specified aligner is in the `$PATH` variable.
 
@@ -548,11 +568,9 @@ We recommend _not_ to modify these options unless you are clearly aware of their
     for storing the annotation database, because SQLite database cannot be created on a shared disks.
     The folder will be created automatically.
 
-`--low_memory`
-    Deprecated, default behaviour since 3.2.
-
 `--high_memory`
-    Cache read alignments instead for making several passes over a BAM file, noticeably increases RAM usage.
+    Cache read alignments instead for making several passes over a BAM file, noticeably increases RAM usage, 
+but may improve running time when disk I/O is relatively slow.
 
 `--min_mapq`
     Filers out all alignments with MAPQ less than this value (will also filter all secondary alignments, as they typically have MAPQ = 0).

isoquant.py

@@ -29,7 +29,8 @@
     NANOPORE_DATA,
     DataSetReadMapper
 )
-from src.dataset_processor import DatasetProcessor
+from src.dataset_processor import DatasetProcessor, PolyAUsageStrategies
+from src.graph_based_model_construction import StrandnessReportingLevel
 from src.long_read_assigner import AmbiguityResolvingMethod
 from src.long_read_counter import COUNTING_STRATEGIES
 from src.input_data_storage import InputDataStorage
@@ -48,10 +49,10 @@ def bool_str(s):
 
 def parse_args(cmd_args=None, namespace=None):
     parser = argparse.ArgumentParser(formatter_class=argparse.RawDescriptionHelpFormatter)
-    ref_args_group = parser.add_argument_group('Reference data:')
-    input_args_group = parser.add_argument_group('Input data:')
-    output_args_group = parser.add_argument_group('Output naming:')
-    pipeline_args_group = parser.add_argument_group('Pipeline options:')
+    ref_args_group = parser.add_argument_group('Reference data')
+    input_args_group = parser.add_argument_group('Input data')
+    output_args_group = parser.add_argument_group('Output naming')
+    pipeline_args_group = parser.add_argument_group('Pipeline options')
 
     other_options = parser.add_argument_group("Additional options:")
     show_full_help = '--full_help' in cmd_args
@@ -150,16 +151,20 @@ def add_hidden_option(*args, **kwargs):  # show command only with --full-help
                           help="gene quantification strategy", type=str, default="with_inconsistent")
     add_additional_option("--report_novel_unspliced", "-u", type=bool_str,
                           help="report novel monoexonic transcripts (true/false), "
-                               "default: False for ONT, True for other data types")
+                               "default: false for ONT, true for other data types")
+    add_additional_option("--report_canonical",  type=str, choices=[e.name for e in StrandnessReportingLevel],
+                          help="reporting level for novel transcripts based on canonical splice sites;"
+                               " default: " + StrandnessReportingLevel.auto.name,
+                          default=StrandnessReportingLevel.only_stranded.name)
+    add_additional_option("--polya_requirement", type=str, choices=[e.name for e in PolyAUsageStrategies],
+                          help="require polyA tails to be present when reporting transcripts; "
+                               "default: auto (requires polyA only when polyA percentage is >= 70%%)",
+                          default=PolyAUsageStrategies.auto.name)
     # OUTPUT PROPERTIES
     pipeline_args_group.add_argument("--threads", "-t", help="number of threads to use", type=int,
                                      default="16")
     pipeline_args_group.add_argument('--check_canonical', action='store_true', default=False,
                                      help="report whether splice junctions are canonical")
-    add_additional_option_to_group(pipeline_args_group, "--report_unstranded",
-                                   help="report transcripts for which the strand cannot be detected "
-                                        "using canonical splice sites",
-                                   action='store_true', default=False)
     pipeline_args_group.add_argument("--sqanti_output", help="produce SQANTI-like TSV output",
                                      action='store_true', default=False)
     pipeline_args_group.add_argument("--count_exons", help="perform exon and intron counting",
@@ -180,8 +185,6 @@ def add_hidden_option(*args, **kwargs):  # show command only with --full-help
                                    help="align reads to reference without running further analysis")
 
     # ADDITIONAL
-    add_additional_option("--low_memory", help="decrease RAM consumption (deprecated, set by default)",
-                          action='store_true', default=True)
     add_additional_option("--high_memory", help="increase RAM consumption (store alignment and the genome in RAM)",
                           action='store_true', default=False)
     add_additional_option("--no_junc_bed", action="store_true", default=False,
@@ -234,8 +237,12 @@ def add_hidden_option(*args, **kwargs):  # show command only with --full-help
                                    type=str, required=True)
         resume_parser.add_argument('--debug', action='store_true', default=argparse.SUPPRESS,
                                    help='Debug log output.')
-        resume_parser.add_argument("--threads", "-t", help="number of threads to use", type=int, default=argparse.SUPPRESS)
-        resume_parser.add_argument("--low_memory", help="decrease RAM consumption (deprecated, set by default)",
+        resume_parser.add_argument("--threads", "-t", help="number of threads to use",
+                                   type=int, default=argparse.SUPPRESS)
+        resume_parser.add_argument("--high_memory",
+                                   help="increase RAM consumption (store alignment and the genome in RAM)",
+                                   action='store_true', default=False)
+        resume_parser.add_argument("--keep_tmp", help="do not remove temporary files in the end",
                                    action='store_true', default=argparse.SUPPRESS)
         resume_parser.add_argument("--keep_tmp", help="do not remove temporary files in the end", action='store_true',
                                    default=argparse.SUPPRESS)
@@ -271,11 +278,12 @@ def check_and_load_args(args, parser):
                 logger.warning("Output folder already contains a previous run, will be overwritten.")
             else:
                 logger.warning("Output folder already contains a previous run, some files may be overwritten. "
-                               "Use --resume to resume a failed run. Use --force to avoid this message. "
-                               "Press Ctrl+C to interrupt the run now.")
+                               "Use --resume to resume a failed run. Use --force to avoid this message.")
+                logger.warning("Press Ctrl+C to interrupt the run now.")
                 delay = 9
                 for i in range(delay):
-                    sys.stdout.write("Resuming the run in %d seconds\r" % (delay - i))
+                    countdown = delay - i
+                    sys.stdout.write("Resuming the run in %d second%s\r" % (countdown, "s" if countdown > 1 else ""))
                     time.sleep(1)
                 logger.info("Overwriting the previous run")
                 time.sleep(1)
@@ -288,9 +296,6 @@ def check_and_load_args(args, parser):
     elif not os.path.exists(args.genedb_output):
         os.makedirs(args.genedb_output)
 
-    if args.high_memory:
-        args.low_memory = False
-
     if not check_input_params(args):
         parser.print_usage()
         exit(-1)
@@ -501,7 +506,7 @@ def set_data_dependent_options(args):
         args.splice_correction_strategy = splice_correction_strategies[args.data_type]
 
     args.resolve_ambiguous = 'monoexon_and_fsm' if args.fl_data else 'default'
-    args.needs_polya_for_construction = False
+    args.requires_polya_for_construction = False
     if args.read_group is None and args.input_data.has_replicas():
         args.read_group = "file_name"
     args.use_technical_replicas = args.read_group == "file_name"
@@ -591,16 +596,26 @@ def set_model_construction_options(args):
                                             'min_known_count', 'min_nonfl_count',
                                             'min_novel_count', 'min_novel_count_rel',
                                             'min_mono_count_rel', 'singleton_adjacent_cov',
-                                            'fl_only', 'novel_monoexonic'))
+                                            'fl_only', 'novel_monoexonic',
+                                            'require_monointronic_polya', 'require_monoexonic_polya',
+                                            'report_canonical'))
     strategies = {
-        'reliable':        ModelConstructionStrategy(2, 0.5, 20,  5, 0.05,  1, 0.1,  0.1,  2, 4, 8, 0.05, 0.05, 50, True, False),
-        'default_pacbio':  ModelConstructionStrategy(1, 0.5, 10,  2, 0.02,  1, 0.05,  0.05,  1, 2, 2, 0.02, 0.005, 100, False, True),
-        'sensitive_pacbio':ModelConstructionStrategy(1, 0.5, 5,   2, 0.005,  1, 0.01,  0.02,  1, 2, 2, 0.005, 0.001, 100, False, True),
-        'default_ont':     ModelConstructionStrategy(1, 0.5, 20,  3, 0.02,  1, 0.05,  0.05,  1, 3, 3, 0.02, 0.02, 10, False, False),
-        'sensitive_ont':   ModelConstructionStrategy(1, 0.5, 20,  3, 0.005,  1, 0.01,  0.02,  1, 2, 3, 0.005, 0.005, 10, False, True),
-        'fl_pacbio':       ModelConstructionStrategy(1, 0.5, 10,  2, 0.02,  1, 0.05,  0.01,  1, 2, 3, 0.02, 0.005, 100, True, True),
-        'all':             ModelConstructionStrategy(0, 0.3, 5,   1, 0.002,  1, 0.01, 0.01, 1, 1, 1, 0.002, 0.001, 500, False, True),
-        'assembly':        ModelConstructionStrategy(0, 0.3, 5,   1, 0.05,  1, 0.01, 0.02,  1, 1, 1, 0.05, 0.01, 50, False, True)
+        'reliable':        ModelConstructionStrategy(2, 0.5, 20,  5, 0.05,  1, 0.1,  0.1,  2, 4, 8, 0.05, 0.05, 50,
+                                                     True, False, True, True, StrandnessReportingLevel.only_canonical),
+        'default_pacbio':  ModelConstructionStrategy(1, 0.5, 10,  2, 0.02,  1, 0.05,  0.05,  1, 2, 2, 0.02, 0.005, 100,
+                                                     False, True, False, True, StrandnessReportingLevel.only_canonical),
+        'sensitive_pacbio':ModelConstructionStrategy(1, 0.5, 5,   2, 0.005,  1, 0.01,  0.02,  1, 2, 2, 0.005, 0.001, 100,
+                                                     False, True, False, False, StrandnessReportingLevel.only_stranded),
+        'default_ont':     ModelConstructionStrategy(1, 0.5, 20,  3, 0.02,  1, 0.05,  0.05,  1, 3, 3, 0.02, 0.02, 10,
+                                                     False, False, True, True, StrandnessReportingLevel.only_canonical),
+        'sensitive_ont':   ModelConstructionStrategy(1, 0.5, 20,  3, 0.005,  1, 0.01,  0.02,  1, 2, 3, 0.005, 0.005, 10,
+                                                     False, True, False, False, StrandnessReportingLevel.only_stranded),
+        'fl_pacbio':       ModelConstructionStrategy(1, 0.5, 10,  2, 0.02,  1, 0.05,  0.01,  1, 2, 3, 0.02, 0.005, 100,
+                                                     True, True, False, False, StrandnessReportingLevel.only_canonical),
+        'all':             ModelConstructionStrategy(0, 0.3, 5,   1, 0.002,  1, 0.01, 0.01, 1, 1, 1, 0.002, 0.001, 500,
+                                                     False, True, False, False, StrandnessReportingLevel.all),
+        'assembly':        ModelConstructionStrategy(0, 0.3, 5,   1, 0.05,  1, 0.01, 0.02,  1, 1, 1, 0.05, 0.01, 50,
+                                                     False, True, False, False, StrandnessReportingLevel.only_stranded)
     }
     strategy = strategies[args.model_construction_strategy]
 
@@ -633,6 +648,13 @@ def set_model_construction_options(args):
         logger.info("Novel unspliced transcripts will not be reported, "
                     "set --report_novel_unspliced true to discover them")
 
+    args.require_monointronic_polya = strategy.require_monointronic_polya
+    args.require_monoexonic_polya = strategy.require_monoexonic_polya
+    args.polya_requirement_strategy = PolyAUsageStrategies[args.polya_requirement]
+    args.report_canonical_strategy = StrandnessReportingLevel[args.report_canonical]
+    if args.report_canonical_strategy == StrandnessReportingLevel.auto:
+        args.report_canonical_strategy = strategy.report_canonical
+
 
 def set_configs_directory(args):
     config_dir = os.path.join(os.environ['HOME'], '.config', 'IsoQuant')
@@ -656,8 +678,6 @@ def set_additional_params(args):
     set_splice_correction_options(args)
 
     args.print_additional_info = True
-    args.no_polya = False
-
     args.indel_near_splice_site_dist = 10
     args.upstream_region_len = 20
 
@@ -673,6 +693,8 @@ def set_additional_params(args):
         args.needs_reference = False
 
     args.simple_models_mapq_cutoff = 30
+    args.polya_percentage_threshold = 0.7
+    args.low_polya_percentage_threshold = 0.1
 
 
 def run_pipeline(args):

src/alignment_processor.py

@@ -229,7 +229,7 @@ def __init__(self, chr_id, bam_pairs, params, illumina_bam, genedb=None, chr_rec
 
         self.bam_merger = BAMOnlineMerger(self.bam_pairs, self.chr_id, 0,
                                           self.bam_pairs[0][0].get_reference_length(self.chr_id),
-                                          multiple_iterators=self.params.low_memory)
+                                          multiple_iterators=not self.params.high_memory)
         self.strand_detector = StrandDetector(self.chr_record)
         self.read_groupper = read_groupper
         self.polya_finder = PolyAFinder(self.params.polya_window, self.params.polya_fraction)
@@ -238,7 +238,7 @@ def __init__(self, chr_id, bam_pairs, params, illumina_bam, genedb=None, chr_rec
 
     def process(self):
         current_region = None
-        alignment_storage = BAMAlignmentStorage(self.bam_merger) if self.params.low_memory else InMemoryAlignmentStorage()
+        alignment_storage = BAMAlignmentStorage(self.bam_merger) if not self.params.high_memory else InMemoryAlignmentStorage()
 
         for bam_index, alignment in self.bam_merger.get():
             if alignment_storage.alignment_is_not_adjacent(alignment):