add an option to tell IsoQuant polyA tails were trimmed

Author: andrewprzh

isoquant.py

@@ -35,6 +35,7 @@
     DataSetReadMapper
 )
 from src.dataset_processor import DatasetProcessor, PolyAUsageStrategies
+from src.alignment_processor import PolyATrimmed
 from src.graph_based_model_construction import StrandnessReportingLevel
 from src.long_read_assigner import AmbiguityResolvingMethod
 from src.long_read_counter import COUNTING_STRATEGIES, CountingStrategy, GroupedOutputFormat, NormalizationMethod
@@ -138,6 +139,9 @@ def add_hidden_option(*args, **kwargs):  # show command only with --full-help
                         help="type of data to process, supported types are: " + ", ".join(DATA_TYPE_ALIASES.keys()))
     input_args_group.add_argument('--stranded',  type=str, help="reads strandness type, supported values are: " +
                         ", ".join(SUPPORTED_STRANDEDNESS), default="none")
+    input_args_group.add_argument('--polya_trimmed', default=None, type=str,
+                                  choices=[e.name for e in PolyATrimmed],
+                                  help="define reads which had polyA tail trimmed")
     input_args_group.add_argument('--fl_data', action='store_true', default=False,
                         help="reads represent FL transcripts; both ends of the read are considered to be reliable")
 

src/alignment_info.py

@@ -8,6 +8,7 @@
 import logging
 
 from .common import get_read_blocks
+from .polya_finder import PolyAInfo
 from .polya_verification import shift_polya, shift_polyt
 
 logger = logging.getLogger('IsoQuant')
@@ -28,7 +29,7 @@ def __init__(self, alignment):
             return
         self.read_start = self.read_exons[0][0]
         self.read_end = self.read_exons[-1][1]
-        self.polya_info = None
+        self.polya_info = PolyAInfo(-1, -1, -1, -1)
         self.exons_changed = False
         self.cage_hits = []
         self.combined_profile = None

src/alignment_processor.py

@@ -33,6 +33,14 @@
 logger = logging.getLogger('IsoQuant')
 
 
+@unique
+class PolyATrimmed(Enum):
+    none = 1
+    all = 2
+    stranded = 3
+    # file = 3
+
+
 @unique
 class AlignmentType(Enum):
     primary = 1
@@ -360,32 +368,34 @@ def process_intergenic(self, alignment_storage, region, skip_read_fraction=1):
             if skip_read_fraction > 1 and counter % skip_read_fraction != 0:
                 continue
 
-            alignment_info.add_polya_info(self.polya_finder, self.polya_fixer)
-            # if self.params.cage:
-            #    alignment_info.add_cage_info(self.cage_finder)
+            if self.params.polya_trimmed == PolyATrimmed.none:
+                alignment_info.add_polya_info(self.polya_finder, self.polya_fixer)
 
             read_assignment = ReadAssignment(read_id, ReadAssignmentType.intergenic,
                                              IsoformMatch(MatchClassification.intergenic))
 
             if alignment_info.exons_changed:
                 read_assignment.add_match_attribute(MatchEvent(MatchEventSubtype.aligned_polya_tail))
-            read_assignment.polyA_found = (alignment_info.polya_info.external_polya_pos != -1 or
-                                           alignment_info.polya_info.external_polyt_pos != -1 or
-                                           alignment_info.polya_info.internal_polya_pos != -1 or
-                                           alignment_info.polya_info.internal_polyt_pos != -1)
+
             read_assignment.polya_info = alignment_info.polya_info
             read_assignment.cage_found = len(alignment_info.cage_hits) > 0
             read_assignment.genomic_region = region
             read_assignment.exons = alignment_info.read_exons
             read_assignment.corrected_exons = corrector.correct_read(alignment_info)
             read_assignment.corrected_introns = junctions_from_blocks(read_assignment.corrected_exons)
-
             read_assignment.read_group = self.read_groupper.get_group_id(alignment, self.bam_merger.bam_pairs[bam_index][1])
             read_assignment.mapped_strand = "-" if alignment.is_reverse else "+"
             read_assignment.strand = self.get_assignment_strand(read_assignment)
             read_assignment.chr_id = self.chr_id
             read_assignment.multimapper = alignment.is_secondary
             read_assignment.mapping_quality = alignment.mapping_quality
+
+            self.add_artificial_polya(read_assignment)
+            read_assignment.polyA_found = (alignment_info.polya_info.external_polya_pos != -1 or
+                                           alignment_info.polya_info.external_polyt_pos != -1 or
+                                           alignment_info.polya_info.internal_polya_pos != -1 or
+                                           alignment_info.polya_info.internal_polyt_pos != -1)
+
             AlignmentCollector.import_bam_tags(alignment, read_assignment, self.params.bam_tags)
             yield read_assignment
 
@@ -429,10 +439,7 @@ def process_genic(self, alignment_storage, gene_info, region, skip_read_fraction
 
             if alignment_info.exons_changed:
                 read_assignment.add_match_attribute(MatchEvent(MatchEventSubtype.aligned_polya_tail))
-            read_assignment.polyA_found = (alignment_info.polya_info.external_polya_pos != -1 or
-                                           alignment_info.polya_info.external_polyt_pos != -1 or
-                                           alignment_info.polya_info.internal_polya_pos != -1 or
-                                           alignment_info.polya_info.internal_polyt_pos != -1)
+
             read_assignment.polya_info = alignment_info.polya_info
             read_assignment.cage_found = len(alignment_info.cage_hits) > 0
             read_assignment.genomic_region = region
@@ -446,11 +453,16 @@ def process_genic(self, alignment_storage, gene_info, region, skip_read_fraction
             read_assignment.strand = self.get_assignment_strand(read_assignment)
             AlignmentCollector.check_antisense(read_assignment)
             AlignmentCollector.import_bam_tags(alignment, read_assignment, self.params.bam_tags)
-
             read_assignment.chr_id = gene_info.chr_id
             read_assignment.multimapper = alignment.is_secondary
             read_assignment.mapping_quality = alignment.mapping_quality
 
+            self.add_artificial_polya(read_assignment)
+            read_assignment.polyA_found = (alignment_info.polya_info.external_polya_pos != -1 or
+                                           alignment_info.polya_info.external_polyt_pos != -1 or
+                                           alignment_info.polya_info.internal_polya_pos != -1 or
+                                           alignment_info.polya_info.internal_polyt_pos != -1)
+
             if self.params.count_exons:
                 read_assignment.exon_gene_profile = alignment_info.combined_profile.read_exon_profile.gene_profile
                 read_assignment.intron_gene_profile = alignment_info.combined_profile.read_intron_profile.gene_profile
@@ -519,6 +531,18 @@ def get_gene_info_for_region(self, current_region):
             gene_info.set_reference_sequence(current_region[0], current_region[1], self.chr_record)
         return gene_info
 
+    def add_artificial_polya(self, read_assignment):
+        if self.params.polya_trimmed == PolyATrimmed.stranded:
+            if read_assignment.strand == '+':
+                read_assignment.polya_info.external_polya_pos = read_assignment.corrected_exons[-1][1] + 1
+            elif read_assignment.strand == '-':
+                read_assignment.polya_info.external_polyt_pos = read_assignment.corrected_exons[0][0] - 1
+        elif self.params.polya_trimmed == PolyATrimmed.all:
+            if read_assignment.mapped_strand == '+':
+                read_assignment.polya_info.external_polya_pos = read_assignment.corrected_exons[-1][1] + 1
+            elif read_assignment.mapped_strand == '-':
+                read_assignment.polya_info.external_polyt_pos = read_assignment.corrected_exons[0][0] - 1
+
     @staticmethod
     def split_coverage_regions(genomic_region, alignment_storage):
         if interval_len(genomic_region) < AlignmentCollector.MIN_REGION_LEN or \