Refactor parent gene

Dockerfile

@@ -9,7 +9,7 @@ RUN git clone git://github.com/samtools/htslib.git
 RUN cd htslib && make install
 
 # bcftools
-RUN git clone git://github.com/samtools/bcftools.git 
+RUN git clone git://github.com/samtools/bcftools.git
 RUN cd bcftools && make
 
 # samtools
@@ -33,7 +33,7 @@ RUN cd augustus && make
 # HDF5
 RUN wget -q http://www.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/hdf5-1.10.1/src/hdf5-1.10.1.tar.gz
 RUN tar xzf hdf5-1.10.1.tar.gz
-RUN cd hdf5-1.10.1 && ./configure --enable-cxx --prefix=/usr 
+RUN cd hdf5-1.10.1 && ./configure --enable-cxx --prefix=/usr
 RUN cd hdf5-1.10.1 && make && make install
 
 # sonLib
@@ -62,7 +62,7 @@ RUN tar xvjf sambamba_v0.6.7_linux.tar.bz2
 
 FROM ubuntu:18.04
 RUN apt-get update
-RUN apt-get install -y wget bedtools bamtools samtools sqlite3 libgsl0-dev libcolamd2 software-properties-common libcurl4-openssl-dev
+RUN apt-get install -y wget bedtools bamtools samtools sqlite3 libgsl0-dev libcolamd2 software-properties-common libcurl4-openssl-dev exonerate
 RUN add-apt-repository -y ppa:deadsnakes/ppa
 RUN apt-get install -y python3.7 python3-pip
 # Kent

cat/__init__.py

@@ -53,7 +53,6 @@
 from .exceptions import *
 
 logger = logging.getLogger('cat')
-logger.setLevel('INFO')
 
 
 ###
@@ -185,7 +184,7 @@ def get_pipeline_args(self):
         args.set('global_near_best', self.global_near_best, True)
         args.set('filter_overlapping_genes', self.filter_overlapping_genes, True)
         args.set('overlapping_gene_distance', self.overlapping_gene_distance, True)
-        
+
         # user specified flags for consensus finding
         args.set('intron_rnaseq_support', self.intron_rnaseq_support, False)
         args.set('exon_rnaseq_support', self.exon_rnaseq_support, False)
@@ -398,6 +397,8 @@ def load_docker(self):
                                            'Either install it or use a different option for --binary-mode.')
             os.environ['SINGULARITY_PULLFOLDER'] = os.path.abspath(self.work_dir)
             os.environ['SINGULARITY_CACHEDIR'] = os.path.abspath(self.work_dir)
+            if os.environ.get('SINGULARITY_IMAGE'):
+                return
             tools.fileOps.ensure_dir(self.work_dir)
             if not os.path.isfile(os.path.join(self.work_dir, 'cat.img')):
                 subprocess.check_call(['singularity', 'pull', '--name', 'cat.img',
@@ -565,7 +566,7 @@ def success(task):
             os.remove(os.path.abspath(task.job_store))
         except OSError:
             pass
-    else: 
+    else:
         cmd = ['toil', 'stats', '--raw', os.path.abspath(task.job_store)]
     raw = tools.procOps.call_proc(cmd)
     parsed = raw[raw.index('{'):raw.rfind('}') + 1]
@@ -1620,14 +1621,13 @@ def get_args(pipeline_args, mode):
         args = tools.misc.HashableNamespace()
         if mode == 'augPB':
             args.tablename = tools.sqlInterface.AugPbAlternativeGenes.__tablename__
+            pb_genomes = set(pipeline_args.target_genomes) & (pipeline_args.isoseq_genomes - {pipeline_args.ref_genome})
             args.gps = {genome: AugustusPb.get_args(pipeline_args, genome).augustus_pb_gp
-                        for genome in pipeline_args.isoseq_genomes}
+                        for genome in pb_genomes}
             args.filtered_tm_gps = {genome: TransMap.get_args(pipeline_args, genome).filtered_tm_gp
-                                    for genome in pipeline_args.isoseq_genomes - {pipeline_args.ref_genome}}
+                                    for genome in pb_genomes}
             args.unfiltered_tm_gps = {genome: TransMap.get_args(pipeline_args, genome).tm_gp
-                                      for genome in pipeline_args.isoseq_genomes - {pipeline_args.ref_genome}}
-            args.chrom_sizes = {genome: GenomeFiles.get_args(pipeline_args, genome).sizes
-                                for genome in pipeline_args.isoseq_genomes}
+                                      for genome in pb_genomes}
             # add the reference annotation as a pseudo-transMap to assign parents in reference
             args.filtered_tm_gps[pipeline_args.ref_genome] = ReferenceFiles.get_args(pipeline_args).annotation_gp
             args.unfiltered_tm_gps[pipeline_args.ref_genome] = ReferenceFiles.get_args(pipeline_args).annotation_gp
@@ -1643,23 +1643,20 @@ def get_args(pipeline_args, mode):
             unfiltered_tm_gp_files[pipeline_args.ref_genome] = ReferenceFiles.get_args(pipeline_args).annotation_gp
             args.filtered_tm_gps = filtered_tm_gp_files
             args.unfiltered_tm_gps = unfiltered_tm_gp_files
-            args.chrom_sizes = {genome: GenomeFiles.get_args(pipeline_args, genome).sizes
-                                for genome in list(pipeline_args.target_genomes) + [pipeline_args.ref_genome]}
         elif mode == 'exRef':
+            exref_genomes = set(pipeline_args.target_genomes) & set(pipeline_args.external_ref_genomes)
             args.tablename = tools.sqlInterface.ExRefAlternativeGenes.__tablename__
             args.gps = {genome: ExternalReferenceFiles.get_args(pipeline_args, genome).annotation_gp
-                        for genome in pipeline_args.external_ref_genomes}
+                        for genome in exref_genomes}
             filtered_tm_gp_files = {genome: TransMap.get_args(pipeline_args, genome).filtered_tm_gp
-                                    for genome in pipeline_args.external_ref_genomes}
+                                    for genome in exref_genomes}
             unfiltered_tm_gp_files = {genome: TransMap.get_args(pipeline_args, genome).tm_gp
-                                      for genome in pipeline_args.external_ref_genomes}
+                                      for genome in exref_genomes}
             # add the reference annotation as a pseudo-transMap to assign parents in reference
             filtered_tm_gp_files[pipeline_args.ref_genome] = ReferenceFiles.get_args(pipeline_args).annotation_gp
             unfiltered_tm_gp_files[pipeline_args.ref_genome] = ReferenceFiles.get_args(pipeline_args).annotation_gp
             args.filtered_tm_gps = filtered_tm_gp_files
             args.unfiltered_tm_gps = unfiltered_tm_gp_files
-            args.chrom_sizes = {genome: GenomeFiles.get_args(pipeline_args, genome).sizes
-                                for genome in list(pipeline_args.target_genomes) + [pipeline_args.ref_genome]}
         else:
             raise Exception('Invalid mode passed to FindDenovoParents')
         return args
@@ -1696,8 +1693,7 @@ def run(self):
             table_target = self.get_table_targets(genome, denovo_args.tablename, pipeline_args)
             filtered_tm_gp = denovo_args.filtered_tm_gps[genome]
             unfiltered_tm_gp = denovo_args.unfiltered_tm_gps[genome]
-            chrom_sizes = denovo_args.chrom_sizes[genome]
-            df = assign_parents(filtered_tm_gp, unfiltered_tm_gp, chrom_sizes, denovo_gp)
+            df = assign_parents(filtered_tm_gp, unfiltered_tm_gp, denovo_gp)
             db = pipeline_args.dbs[genome]
             with tools.sqlite.ExclusiveSqlConnection(db) as engine:
                 df.to_sql(denovo_args.tablename, engine, if_exists='replace')
@@ -2298,7 +2294,7 @@ def output(self):
     def run(self):
         pipeline_args = self.get_pipeline_args()
         luigi_stats = tools.sqlInterface.load_luigi_stats(pipeline_args.stats_db, 'stats')
-        
+
         try:
             toil_stats = tools.sqlInterface.load_luigi_stats(pipeline_args.stats_db, 'toil_stats')
         except ValueError:
@@ -3184,7 +3180,7 @@ def construct_consensus_gp_as(has_rna, has_pb):
 '''
 
 
-snake_composite = '''track hubCentral
+snake_composite = '''track cactus
 compositeTrack on
 shortLabel Cactus
 longLabel Cactus Alignment Tracks
@@ -3199,11 +3195,11 @@ def construct_consensus_gp_as(has_rna, has_pb):
 visibility full
 type bigBed 3
 
-    track hubCentralAlignments
+    track cactusAlignments
     shortLabel Alignments
     view Alignments
     visibility full
-    subTrack hubCentral
+    subTrack cactus
 
 '''
 
@@ -3212,7 +3208,7 @@ def construct_consensus_gp_as(has_rna, has_pb):
         shortLabel {genome}
         otherSpecies {genome}
         visibility {visibility}
-        parent hubCentralAlignments off
+        parent cactusAlignments off
         priority 3
         bigDataUrl {hal_path}
         type halSnake

cat/chaining.py

@@ -1,5 +1,8 @@
 """"
-Toil program to generate UCSC chains and nets between two genomes in a HAL file.
+Toil program to generate UCSC chains and nets between one reference genome and N target genomes in a HAL file.
+
+The reference genome is split up by contig/chromosome, and those are chained against all the things they align to
+in each genome specified.
 """
 import argparse
 import collections

cat/classify.py

@@ -42,8 +42,6 @@
 2) mRNA
 
 """
-import itertools
-
 import pandas as pd
 
 import tools.bio
@@ -244,29 +242,28 @@ def find_indels(tx, psl, aln_mode):
     """
     def convert_coordinates_to_chromosome(left_pos, right_pos, coordinate_fn, strand):
         """convert alignment coordinates to target chromosome coordinates, inverting if negative strand"""
-        left_chrom_pos = coordinate_fn(left_pos)
-        assert left_chrom_pos is not None
-        right_chrom_pos = coordinate_fn(right_pos)
-        if right_chrom_pos is None:
-            right_chrom_pos = coordinate_fn(right_pos - 1)
-            if strand == '-':
-                left_chrom_pos += 1
-            else:
-                left_chrom_pos -= 1
-        assert right_chrom_pos is not None
-        if strand == '-':
-            left_chrom_pos, right_chrom_pos = right_chrom_pos, left_chrom_pos
-        assert right_chrom_pos >= left_chrom_pos
+        # this is a deletion, and we don't need to adjust
+        if left_pos == right_pos:
+            left_chrom_pos = right_chrom_pos = coordinate_fn(left_pos)
+            # for + strand insertion, just translate
+        elif strand == '+':
+            right_chrom_pos = coordinate_fn(right_pos)
+            left_chrom_pos = coordinate_fn(left_pos)
+        # for - strand insertion, translate but move outwards
+        else:
+            right_chrom_pos = coordinate_fn(left_pos) + 1
+            left_chrom_pos = coordinate_fn(right_pos - 1)
         return left_chrom_pos, right_chrom_pos
 
     def parse_indel(left_pos, right_pos, coordinate_fn, tx, offset, gap_type):
         """Converts either an insertion or a deletion into a output transcript"""
         left_chrom_pos, right_chrom_pos = convert_coordinates_to_chromosome(left_pos, right_pos, coordinate_fn,
                                                                             tx.strand)
         if left_chrom_pos is None or right_chrom_pos is None:
-            assert aln_mode == 'CDS'
             return None
 
+        assert right_chrom_pos >= left_chrom_pos
+
         if left_chrom_pos > tx.thick_start and right_chrom_pos < tx.thick_stop:
             indel_type = 'CodingMult3' if offset % 3 == 0 else 'Coding'
         else:
@@ -295,13 +292,20 @@ def parse_indel(left_pos, right_pos, coordinate_fn, tx, offset, gap_type):
         t_offset = t_start - block_size - t_pos
         assert (q_offset >= 0 and t_offset >= 0)
         if q_offset != 0:  # query insertion -> insertion in target sequence
-            left_pos = q_start - q_offset
-            right_pos = q_start
+            if tx.strand == '+':
+                left_pos = q_start - q_offset
+                right_pos = q_start
+            else:
+                left_pos = psl.q_size - q_start - 1
+                right_pos = psl.q_size - q_start + q_offset
             row = parse_indel(left_pos, right_pos, coordinate_fn, tx, q_offset, 'Insertion')
             if row is not None:
                 r.append(row)
         if t_offset != 0:  # target insertion -> insertion in reference sequence
-            left_pos = right_pos = q_start
+            if tx.strand == '+':
+                left_pos = right_pos = q_start
+            else:
+                left_pos = right_pos = psl.q_size - q_start
             row = parse_indel(left_pos, right_pos, coordinate_fn, tx, t_offset, 'Deletion')
             if row is not None:
                 r.append(row)

cat/consensus.py

@@ -648,9 +648,11 @@ def add_exref_ids(s):
                                   'source_gene_common_name': s.CommonName}
 
         # bring in extra tags for exRef
-        if 'exRef' in denovo_tx_modes:
-            for key, val in tools.misc.parse_gff_attr_line(exref_annot.loc[aln_id].ExtraTags).items():
-                denovo_tx_dict[aln_id][key] = val
+        if tools.nameConversions.aln_id_is_exref(aln_id):
+            tags = exref_annot.loc[aln_id].ExtraTags
+            if len(tags) > 0:
+                for key, val in tools.misc.parse_gff_attr_line(tags).items():
+                    denovo_tx_dict[aln_id][key] = val
 
         # record some metrics
         metrics['denovo'][tx_mode][s.TranscriptClass.replace('_', ' ').capitalize()] += 1
@@ -907,12 +909,15 @@ def convert_attrs(attrs, id_field):
         else:
             attrs['Name'] = attrs['gene_id']
         # convert empty strings into nan
-        attrs = {x: y if y != '' else 'nan' for x, y in attrs.items()}
-        # don't include the support vectors in the string, they will be placed in their respective places
-        attrs_str = ['='.join([key, str(val).replace('=', '_')]) for key, val in sorted(attrs.items()) if 'support' not in key]
-        # explicitly escape any semicolons that may exist in the input strings
-        attrs_str = [x.replace(';', '%3B') for x in attrs_str]
-        return score, ';'.join(attrs_str)
+        attrs_str = []
+        for key, val in attrs.items():
+            val = str(val)
+            if len(val) == 0:
+                val = 'nan'
+            val = val.replace('=', '%3D').replace(';', '%3B')
+            key = key.replace('=', '%3D').replace(';', '%3B')
+            attrs_str.append(f"{key}={val}")
+        return score, ";".join(attrs_str)
 
     def find_feature_support(attrs, feature, i):
         """Extracts the boolean value from the comma delimited string"""
@@ -929,7 +934,7 @@ def find_all_tx_modes(attrs_list):
             for attrs in attrs_list:
                 tx_modes.update(attrs['transcript_modes'].split(','))
             return ','.join(tx_modes)
-        
+
         intervals = set()
         for tx in tx_objs:
             intervals.update(tx.exon_intervals)

cat/hints_db.py

@@ -307,7 +307,7 @@ def run_protein_aln(job, protein_subset, genome_fasta_file_id):
     # perform alignment
     tmp_exonerate = tools.fileOps.get_tmp_toil_file()
     cmd = ['exonerate', '--model', 'protein2genome', '--showvulgar', 'no', '--showalignment', 'no',
-           '--showquerygff', 'yes', genome_fasta, protein_fasta]
+           '--showquerygff', 'yes', protein_fasta, genome_fasta]
     tools.procOps.run_proc(cmd, stdout=tmp_exonerate)
     return job.fileStore.writeGlobalFile(tmp_exonerate)