support more cmap; check mcl

Author: zhangrengang

SOI-tools.md

@@ -27,13 +27,17 @@ For example, `CHR_Ae11-Ae15-Ca2-Ca7-Vv1_143_283.concat.treefile`, `Ae11-Ae15-Ca2
 
 #### Allele identification ####
 For diploid and polyploid assembly, we may aim to identify allele genes. As allele genes are not always syntenic (e.g., non-syntenic due to samll-scale inversion), we combine synteny and orthology to identify a full set of allele genes.
+
 After `soi filter` has identified orthologous synteny in the [example pipeline](https://github.com/zhangrengang/evolution_example), we run
-`
+```
 soi-syn retrieve_allele collinearity.ortho ../OrthoFinder/OrthoFinder/Results_*/ ../all_species_gene.gff sps="Panax_ginseng Panax_notoginseng" min_block=10 win_size=10 > allele.txt
-`
-`sps="Panax_ginseng Panax_notoginseng"` sets one or more target species or pseudo-species; species labels should be different for different subgenomes (label gene ID like `SP1|sgA.g1, SP1|sgB.g1`);
-`min_block=10` sets minimum length of syntenic blocks;
-`win_size=10` sets upstream and downstream 10 genes to retrieve orthologs using syntenic genes as anchors.
+```
+`sps="Panax_ginseng Panax_notoginseng"` sets one or more target species or pseudo-species (space seperated); 
+species labels need to be different for different subgenomes (label gene ID like `SP1|sgA.g1`, `SP1|sgB.g1`);
+
+`min_block=10` sets minimum length of syntenic blocks; increasing the value will result in more reliable alleles;
+
+`win_size=10` sets upstream and downstream 10 genes to retrieve orthologs using syntenic genes as anchors; increasing the value will retrieve more alleles.
 
 The output file is like:
 ```
@@ -55,7 +59,7 @@ chrom   idx     Panax_ginseng   Panax_notoginseng       Panax_ginseng   Panax_no
 ```
 The 3rd and 4th columns indicate allelic gene pairs. 
 The last column indicates the source of alleles; for example, `orthology:15422` means ortholog #15422, 
-`synteny-1:51` means from syntenic block #1 and ortholog #51; `synteny-1:None` means that it is not pre-inferred as an ortholog.
+`synteny-1:51` means it is from syntenic block #1 and ortholog #51; `synteny-1:None` means that it is not pre-inferred as an ortholog.
 The column number will extend for polyploids.
 
 #### Orthology format conversion ####

soi/dot_plotter.py

@@ -17,13 +17,33 @@
 from .RunCmdsMP import logger
 from .WGDI import AK
 from .ploidy_plotter import add_ploidy_opts, get_ploidy, plot_bars
-from .mcscan import Collinearity, Gff, XCollinearity
+from .mcscan import Collinearity, XGff, XCollinearity
+
+cmaps = {
+            'viridis', 'plasma', 'inferno', 'magma', 'cividis',
+            'Greys', 'Purples', 'Blues', 'Greens', 'Oranges', 'Reds',
+            'YlOrBr', 'YlOrRd', 'OrRd', 'PuRd', 'RdPu', 'BuPu',
+            'GnBu', 'PuBu', 'YlGnBu', 'PuBuGn', 'BuGn', 'YlGn',
+            'binary', 'gist_yarg', 'gist_gray', 'gray', 'bone', 'pink',
+            'spring', 'summer', 'autumn', 'winter', 'cool', 'Wistia',
+            'hot', 'afmhot', 'gist_heat', 'copper',
+            'PiYG', 'PRGn', 'BrBG', 'PuOr', 'RdGy', 'RdBu',
+            'RdYlBu', 'RdYlGn', 'Spectral', 'coolwarm', 'bwr', 'seismic',
+            'berlin', 'managua', 'vanimo',
+            'twilight', 'twilight_shifted', 'hsv',
+            'Pastel1', 'Pastel2', 'Paired', 'Accent',
+            'Dark2', 'Set1', 'Set2', 'Set3',
+            'tab10', 'tab20', 'tab20b', 'tab20c',
+            'flag', 'prism', 'ocean', 'gist_earth', 'terrain', 'gist_stern',
+            'gnuplot', 'gnuplot2', 'CMRmap', 'cubehelix', 'brg',
+            'gist_rainbow', 'rainbow', 'jet', 'turbo', 'nipy_spectral',
+            'gist_ncar'}
 
 
 def dotplot_args(parser):
 	parser.add_argument('-s', metavar='FILE', type=str, required=True, nargs='+',
 						help="syntenic block file (*.collinearity, output of MCSCANX/WGDI)[required]")
-	parser.add_argument('-g', metavar='FILE', type=str, required=True,
+	parser.add_argument('-g', metavar='FILE', type=str, required=True, nargs='+',
 						help="gene annotation gff file (*.gff, one of MCSCANX/WGDI input)[required]")
 	parser.add_argument('-c', metavar='FILE', type=str, required=True,
 						help="chromosomes config file (*.ctl, same format as MCSCANX dotplotter)[required]")
@@ -99,8 +119,8 @@ def dotplot_args(parser):
 						  help="plot histogram or not [default=%(default)s]")
 	group_ks.add_argument('--max-ks', metavar='Ks', type=float, default=1,
 						  help="max Ks (x limit) [default=%(default)s]")
-	group_ks.add_argument('--ks-cmap', metavar='Ks', type=float, nargs='+', default=None,
-						  help="color map for Ks. [default=%(default)s]")
+	group_ks.add_argument('--ks-cmap', metavar='Ks', nargs='+', default=None, # type=float, 
+						  help="color map for Ks, format: `jet` or `0.2 0.6 1 ...`. [default=%(default)s]")
 	group_ks.add_argument('--ks-step', metavar='Ks', type=float, default=0.02,
 						  help="Ks step of histogram [default=%(default)s]")
 	group_ks.add_argument('--use-median', action='store_true', default=False,
@@ -327,7 +347,9 @@ def plot_blocks(blocks, outplots, ks=None, max_ks=None, ks_hist=False, ks_cmap=N
 			min_ks = min([v for v in Ks if v >= 0])
 		except ValueError:  # ValueError: min() arg is an empty sequence
 			min_ks = 0
-		if ks_cmap:
+		if ks_cmap and ks_cmap[0] in cmaps:
+			cmap = ks_cmap[0]
+		elif ks_cmap:
 			cmap = create_ks_map(ks_cmap, min_ks, max_ks)
 		else:
 			cmap = cm.jet
@@ -649,6 +671,7 @@ def create_ks_map(ks_map, min_ks, max_ks):
 	import numpy as np
 	from matplotlib import cm
 	from matplotlib.colors import ListedColormap, LinearSegmentedColormap
+	ks_map = list(map(float, ks_map))
 	length = 256
 	maps = _norm_map(ks_map, min_ks, max_ks, length)
 	print(ks_map, min_ks, max_ks, maps)
@@ -683,7 +706,7 @@ def parse_gff(gff, chrs1, chrs2):
 	coord_graph1 = nx.Graph()
 	coord_graph2 = nx.Graph()
 	d_gff = {}
-	for line in Gff(gff):
+	for line in XGff(gff):
 		if not line.chrom in chrs:
 			continue
 		key = (line.species, line.chrom)

soi/mcscan.py

@@ -240,7 +240,7 @@ def _parse_list(self, _collinearities):
 	def _parse(self):
 		if self.orthologs is not None:
 			ortholog_pairs = set(XOrthology(self.orthologs, **self.kargs))
-			logger.info('\t{} homologous pairs'.format(len(ortholog_pairs)))
+		#	logger.info('\t{} homologous pairs'.format(len(ortholog_pairs)))
 		logger.info('parsing synteny from {} collinearity files: {} ...'.format(
 			len(self.collinearities), self.collinearities[:3]))
 		nblock, ngene = 0, 0
@@ -328,8 +328,8 @@ def _parse_list(self, _orthologs):
 	def _parse(self):
 		'''yield Pair object'''
 		i = 0
+		logger.info('parsing orthology: {} ...'.format(self.orthologs))
 		for ortholog in self.orthologs:
-			logger.info('parsing orthology: {} ...'.format(ortholog))
 			if os.path.isdir(ortholog):
 				# SonicParanoid / OrthoFinder
 				parser = lazy_orthologs(ortholog)
@@ -968,21 +968,23 @@ def parse_gff(self):
 		genes = set([])
 		d_chr = {}
 		d_length = {}
-		for line in open(self.gff):
-			line = lazy_decode(line)
-			temp = line.rstrip().split('\t')
-			if not temp or line.startswith('#'):
-				continue
-			chr, gene, start, end = temp[:4]
-			if gene in genes:  # remove repeat
-				continue
-			genes.add(gene)
-			try:
-				strand = temp[4]
-			except IndexError:
-				strand = None
-			start, end = list(map(int, [start, end]))
-			g = Gene((gene, chr, start, end, strand))
+		for line in XGff(self.gff):
+			# line = lazy_decode(line)
+			# temp = line.rstrip().split('\t')
+			# if not temp or line.startswith('#'):
+				# continue
+			# chr, gene, start, end = temp[:4]
+			# if gene in genes:  # remove repeat
+				# continue
+			# genes.add(gene)
+			# try:
+				# strand = temp[4]
+			# except IndexError:
+				# strand = None
+			# start, end = list(map(int, [start, end]))
+			chr, gene, start, end, strand = \
+				line.chrom, line.gene, line.start, line.end, line.strand
+			g = line.Gene #((gene, chr, start, end, strand))
 			try:
 				d_chr[chr] += [g]
 			except KeyError:
@@ -1088,6 +1090,18 @@ def anchors2bed(collinearity, gff, chrmap, left_anchors, right_anchors, outbed=s
 	print('\t'.join(line), file=outbed)
 
 
+class XGff(XOrthology):
+	def __init__(self, gffs):
+		self.gffs = self._parse_list(gffs)
+	def __iter__(self):
+		return self._parse()
+	def _parse(self):
+		logger.info('parsing gff files: {} ...'.format(self.gffs))
+		for gff in self.gffs:
+			for line in Gff(gff):
+				yield line
+
+
 class Gff:
 	'''Gff parser'''
 
@@ -1853,6 +1867,7 @@ def parse_group(groups):
 
 def cluster_by_mcl(collinearities, orthologs=None, inflation=2,
 				   outgroup=None, ingroup=None, outpre='cluster'):
+	check_cmd('mcl')
 	ingroup = set(parse_group(ingroup))
 	outgroup = set(parse_group(outgroup))
 	logger.info('outgroup: {}'.format(outgroup))
@@ -1890,6 +1905,9 @@ def cluster_by_mcl(collinearities, orthologs=None, inflation=2,
 	nc = len([1 for line in open(cluster)])
 	logger.info('{} syntenic gene pairs reslut in {} SOGs'.format(np, nc))
 
+def check_cmd(cmd):
+	logger.info('checking `{}`'.format(cmd))
+	run_cmd(cmd, log=True, )
 
 def test_closest(collinearity, kaks, spsd, min_size=0):
 	ColinearGroups(collinearity, spsd, kaks=kaks,