Depicting deletions with plotsr

README.md

@@ -127,6 +127,7 @@ Valid values for annotation type: SYN, INV, TRA, INVTR, DUP, INVDP. Here:
 | INVTR    | Inverted translocation |
 | DUP      | Duplication |
 | INVDP    | Inverted duplication |
+| DEL | Deletion | 
 
 <b><i>NOTE</b>: The BEDPE file must have syntenic region annotations. These are required to group homologous chromosomes from different genomes. Syntenic regions can only be between homologous chromosomes. In case, syntenic regions between homologous chromosomes are not available, then entire homologous chromosomes can be added as syntenic in the BEDPE file manually to allow clustering of homologous chromosomes by plotsr. While plotting, use the `--nosyn` option to skip plotting of these manually added syntenic regions.  </i>
 
@@ -254,4 +255,4 @@ Additional parameters (colors, spacing, legends) of the plot can be adjusted by
 ## Citation:
 If you find plotsr helpful, please [cite](https://doi.org/10.1093/bioinformatics/btac196):
 
-`Manish Goel, Korbinian Schneeberger, plotsr: visualizing structural similarities and rearrangements between multiple genomes, Bioinformatics, 2022; btac196, https://doi.org/10.1093/bioinformatics/btac196`
+`Manish Goel, Korbinian Schneeberger, plotsr: visualizing structural similarities and rearrangements between multiple genomes, Bioinformatics, 2022; btac196, https://doi.org/10.1093/bioinformatics/btac196`

config/base.cfg

@@ -3,10 +3,13 @@ syncol:#CCCCCC
 invcol:#FFA500
 tracol:#9ACD32
 dupcol:#00BBFF
+delcol:#FF0055
 synlwd:0                ## Line width for syntenic annotations
 invlwd:0.1              ## Line width for inversions
 tralwd:0.1              ## Line width for translocations
 duplwd:0.1              ## Line width for duplications
+delwd:0.1              ## Line width for deletion
+
 alpha:0.8
 
 ## Margins and dimensions:
@@ -26,4 +29,4 @@ norm:T                  ## For each chromosome, independently normalise the y-ax
 
 ## Axis
 maxl:-1                 ## Manually set maximum chromosome position. Use `-1` for automatic selection. Does not work with --itx
-genname:T               ## Write genome names adjacent to the chromosome (T) or not (F)
+genname:T               ## Write genome names adjacent to the chromosome (T) or not (F)

plotsr/scripts/func.py

@@ -39,8 +39,8 @@
            "i9": "caretright (centered at base)",
            "i10": "caretup (centered at base)",
            "i11": "caretdown"}
-VARS = ['SYN', 'INV', 'TRANS', 'INVTR', 'DUP', 'INVDP']
-COLORS = ['#DEDEDE', '#FFA500', '#9ACD32', '#00BBFF']
+VARS = ['SYN', 'INV', 'TRANS', 'INVTR', 'DUP', 'INVDP','DEL']
+COLORS = ['#DEDEDE', '#FFA500', '#9ACD32', '#00BBFF','#FF0055']
 
 FONT_NAMES = []
 for fn in matplotlib.font_manager.findSystemFonts():
@@ -194,10 +194,12 @@ def readbasecfg(f, v):
     cfg['invcol'] = '#FFA500'
     cfg['tracol'] = '#9ACD32'
     cfg['dupcol'] = '#00BBFF'
+    cfg['delcol'] = '#FF0055'
     cfg['synlwd'] = 0
     cfg['invlwd'] = 0.1
     cfg['tralwd'] = 0.1
     cfg['duplwd'] = 0.1
+    cfg['delwd'] = 0.1
     cfg['alpha'] = 0.8
 
     # Set chromosome margins
@@ -234,15 +236,15 @@ def readbasecfg(f, v):
             if line[0] not in cfgk:
                 logger.error("{} is not a valid config parameter. Using default value.".format(line[0]))
                 continue
-            if line[0] in ['syncol', 'invcol', 'tracol', 'dupcol']:
+            if line[0] in ['syncol', 'invcol', 'tracol', 'dupcol', 'delcol']:
                 try:
                     if line[1] == '#': matplotlib.colors.to_rgb(line[1])
                     else: matplotlib.colors.to_hex(line[1])
                 except ValueError:
                     logger.error("Error in using colour: {} for {}. Use correct hexadecimal colours or named colours defined in matplotlib (https://matplotlib.org/stable/gallery/color/named_colors.html). Using default value.".format(line[1], line[0]))
                     continue
                 cfg[line[0]] = line[1]
-            elif line[0] in ['synlwd', 'invlwd', 'tralwd', 'duplwd', 'alpha', 'chrmar', 'exmar', 'bboxmar', 'genlegcol', 'marginchr', 'maxl']:
+            elif line[0] in ['synlwd', 'invlwd', 'tralwd', 'duplwd', 'delwd','alpha', 'chrmar', 'exmar', 'bboxmar', 'genlegcol', 'marginchr', 'maxl']:
                 try:
                     float(line[1])
                 except ValueError:
@@ -442,7 +444,7 @@ def readsyriout(f):
     # Reads syri.out. Select: achr, astart, aend, bchr, bstart, bend, srtype
     logger = logging.getLogger("readsyriout")
     syri_regs = deque()
-    skipvartype = ['CPG', 'CPL', 'DEL', 'DUPAL', 'HDR', 'INS', 'INVAL', 'INVDPAL', 'INVTRAL', 'NOTAL', 'SNP', 'SYNAL', 'TDM', 'TRANSAL']
+    skipvartype = ['CPG', 'CPL', 'DUPAL', 'HDR', 'INS', 'INVAL', 'INVDPAL', 'INVTRAL', 'NOTAL', 'SNP', 'SYNAL', 'TDM', 'TRANSAL']
     logger.info('Reading input files generated by syri.')
     with open(f, 'r') as fin:
         for line in fin:
@@ -1082,6 +1084,8 @@ def filterinput(args, df, chrid, itx=False):
         df = df.loc[~df['type'].isin(['TRANS', 'INVTR'])]
     if args.nodup:
         df = df.loc[~df['type'].isin(['DUP', 'INVDP'])]
+    if args.nodel:
+        df = df.loc[df['type'] != 'DEL']
     df.sort_values(['bchr', 'bstart', 'bend'], inplace=True)
     df.sort_values(['achr', 'astart', 'aend'], inplace=True)
     return df
@@ -1583,22 +1587,25 @@ def annotodict(anno):
     adinvlab = False
     adtralab = False
     adduplab = False
+    addelab = False
     svlabels = dict()
-    legenddict = {'SYN': adsynlab, 'INV': adinvlab, 'TRANS': adtralab, 'DUP': adduplab}
+    legenddict = {'SYN': adsynlab, 'INV': adinvlab, 'TRANS': adtralab, 'DUP': adduplab, 'DEL': addelab}
     for s in range(len(alignments)):
         df = deepcopy(alignments[s][1])
         df.loc[df['type'] == 'INVTR', 'type'] = 'TRANS'
         df.loc[df['type'] == 'INVDP', 'type'] = 'DUP'
         coldict = {'SYN': cfg['syncol'],
                    'INV': cfg['invcol'],
                    'TRANS': cfg['tracol'],
-                   'DUP': cfg['dupcol']}
+                   'DUP': cfg['dupcol'],
+                   'DEL': cfg['delcol']}
         lwddict = {'SYN': cfg['synlwd'],
                    'INV': cfg['invlwd'],
                    'TRANS': cfg['tralwd'],
-                   'DUP': cfg['duplwd']}
+                   'DUP': cfg['duplwd'],
+                   'DEL': cfg['delwd']}
         # df['col'] = [coldict[c] for c in df['type']]
-        labdict = {'SYN': 'Syntenic', 'INV': 'Inversion', 'TRANS': 'Translocation', 'DUP': 'Duplication'}
+        labdict = {'SYN': 'Syntenic', 'INV': 'Inversion', 'TRANS': 'Translocation', 'DUP': 'Duplication', 'DEL': 'Deletion'}
         df['lab'] = [labdict[c] for c in df['type']]
         df.loc[df.duplicated(['lab']), 'lab'] = ''
         # df['lw'] = 0
@@ -1666,7 +1673,7 @@ def annotodict(anno):
                     if not legenddict[row.type]:
                         svlabels[row.type] = l
                         legenddict[row.type] = True
-    return ax, [svlabels[i] for i in ['SYN', 'INV', 'TRANS', 'DUP'] if i in svlabels]
+    return ax, [svlabels[i] for i in ['SYN', 'INV', 'TRANS', 'DUP', 'DEL'] if i in svlabels]
 # END
 
 

plotsr/scripts/plotsr.py

@@ -47,26 +47,26 @@ def plotsr(args):
     ## Validate input
     if args.sr is None and args.bp is None:
         logger.error("No structural annotations provided. Use --sr or -bp to provide path to input files")
-        sys.exit()
+        sys.exit(1)
 
     if args.sr is not None and args.bp is not None:
         logger.error("Both --sr and --bp cannot be used. Use single file type for all input structural annotations files. User converter to reformat BEDPE/syri.out files")
-        sys.exit()
+        sys.exit(1)
 
     # Check if both --chr and --reg are defined
     if args.chr is not None and args.reg is not None:
         logger.error("Both --chr and --reg are provided. Only one parameter can be provided at a time. Exiting.")
-        sys.exit()
+        sys.exit(1)
 
     # Check if both --chr and --chrord are defined
     if args.chr is not None and args.chrord is not None:
         logger.error("Both --chr and --chrord are provided. Only one parameter can be provided at a time. Exiting.")
-        sys.exit()
+        sys.exit(1)
 
     # Check if --rtr is used without --reg
     if args.rtr and args.reg is None:
         logger.error("Cannot use --rtr without --reg. Exiting.")
-        sys.exit()
+        sys.exit(1)
 
 
     ###################################################################
@@ -151,7 +151,7 @@ def plotsr(args):
         # Check that the chrorder file contains all chromosomes
         if len(chrs) != len(cs):
             logger.error("Number of chromsomes in {} is less than the number of chromsomes in the alignment file {}. Either list the order of all chromosomes or use --chr if chromosome selection is requires. Exiting.".format(args.chrord.name, alignments[0][0]))
-            sys.exit()
+            sys.exit(1)
 
     chrgrps = OrderedDict()
     for c in chrs:
@@ -192,7 +192,7 @@ def plotsr(args):
         g = set(df.loc[invindex, 'bstart'] < df.loc[invindex, 'bend'])
         if len(g) == 2:
             logger.error("Inconsistent coordinates in input file {}. For INV, INVTR, INVDUP annotations, either bstart < bend for all annotations or bstart > bend for all annotations. Mixing is not permitted.".format(alignments[i][0]))
-            sys.exit()
+            sys.exit(1)
         elif False in g:
             continue
         df.loc[invindex, 'bstart'] = df.loc[invindex, 'bstart'] + df.loc[invindex, 'bend']
@@ -217,7 +217,7 @@ def plotsr(args):
             fig = plt.figure(figsize=[W, H])
     except Exception as e:
         logger.error("Error in initiliazing figure. Try using a different backend.\n{}".format(e.with_traceback()))
-        sys.exit()
+        sys.exit(1)
     ax = fig.add_subplot(111, frameon=False)
 
     allal = pdconcat([alignments[i][1] for i in range(len(alignments))])
@@ -242,6 +242,8 @@ def plotsr(args):
         labelcnt += 1
     if 'DUP' in allal['type'].array or 'INVDP' in allal['type'].array:
         labelcnt += 1
+    if 'DEL' in allal['type'].array:
+        labelcnt += 1
 
     ## Draw Axes
     ax = drawax(ax, chrgrps, chrlengths, V, S, cfg, ITX, minl=minl, maxl=maxl, chrname=CHRNAME)
@@ -311,6 +313,7 @@ def main():
     filtering.add_argument('--noinv', help='Do not plot inversions', default=False, action='store_true')
     filtering.add_argument('--notr', help='Do not plot translocations regions', default=False, action='store_true')
     filtering.add_argument('--nodup', help='Do not plot duplications regions', default=False, action='store_true')
+    filtering.add_argument('--nodel', help='Do not plot deletions', default=False, action='store_true')
     filtering.add_argument('-s', help='minimum size of a SR to be plotted', type=int, default=10000)
 
     plotting = parser.add_argument_group("Plot adjustment")
@@ -330,8 +333,6 @@ def main():
     other.add_argument('--version', action='version', version='{version}'.format(version=__version__))
     parser._action_groups.append(other)
 
-    # args = parser.parse_args([]) # TODO: Delete this line
     args = parser.parse_args()
-    # args = parser.parse_args('--sr col_lersyri.out --sr ler_cvisyri.out --sr cvi_erisyri.out --sr eri_shasyri.out --sr sha_kyosyri.out --sr kyo_an1syri.out --sr an1_c24syri.out --genomes genomes.txt  --chr Chr3 -S 1 -o ampril_col0_chr3.png -W 5 -H 3 -f 8 --cfg base.cfg'.split())
     plotsr(args)
 # END