Merge branch 'roman' into roman2

CIAlign/argP.py

@@ -516,11 +516,30 @@ def getParser():
                  action="store_true",
                  help="""Plot a mini alignment showing positions which are \
                  identical to or differ from the consensus. \
-                 Also has the option to plot mini alignment showing \
-                 positions based on their numerical values when \
-                 when compared via a matrix to the consensus. \
                  Default: %(default)s""")
 
+    optional.add("--plot_consensus_similarity",
+                  dest='plot_consensus_similarity', action='store_true',
+                 help="""Plot mini alignment showing \
+                 positions based on their score when \
+                 when compared via a substitution matrix to the consensus. \
+                 Default: %(default)s""")
+
+    optional.add("--plot_substitution_matrix", dest='plot_substitution_matrix',
+                 type=str, default='B',
+                 help="""Substitution matrix to use for similarity plots. \
+                 Default: %(default)s""")
+
+    optional.add("--plot_identity_palette", dest='plot_identity_palette',
+                 type=str, default='terrain_r',
+                 help="""Matplotlib palette name for identity mini alignments. \
+                 Default: %(default)s""")
+
+    optional.add("--plot_similarity_palette", dest='plot_similarity_palette',
+                 type=str, default='summer_r',
+                 help="""Matplotlib palette name for similarity mini \
+                 alignments. Default: %(default)s""")
+
     optional.add("--plot_dpi", dest="plot_dpi",
                  type=int, default=300, metavar="(int)",
                  help="DPI for mini alignments. Default: %(default)s")

CIAlign/consensusSeq.py

@@ -18,6 +18,7 @@
 import scipy.stats
 import copy
 import operator
+import pandas as pd
 matplotlib.use('Agg')
 
 
@@ -694,81 +695,127 @@ def calcConservationAli(alignment, typ):
     return (heights, ents)
 
 
-def compareAlignmentConsensus(arr, typ, booleanOrSimilarity="Boolean", MatrixName="B"):
-    consensus, _ = np.array(findConsensus(arr, '', consensus_type='majority_nongap')  )
-    if booleanOrSimilarity == "Boolean":
-        '''
-        Compares the alignment of the inputted array to the consensus of that array, and outputs a boolean array.
-
-        alignment: arr
-          The alignment stored as a numpy array
+def compareAlignmentConsensus(arr, typ, booleanOrSimilarity="Boolean",
+                              MatrixName="B"):
+    '''
+    Compares the alignment of the inputted array to the consensus of that
+    array, and will either output a boolean array, or will use a matrix 
+    (can be specified) to output an array containing the scores of the 
+    sequence compared to the consensus.
+
+
+
+        Parameters
+        ----------
+        arr: np.array
+            the sequence thats is to be aligned
+
+        typ: str
+            nt or aa
+        
+        booleanOrSimilarity: str
+            boolean or similarity (default = 'Boolean')
+
+        MatrixName: str
+            the specified matrix name (default = 'B')
+            
+        Returns
+        -------
+        new_arr: np.array
+            A boolean array containing the values of the 
+                        sequence compared to the consensus
+        new_Sarr: np.array
+            A integer array containing the values of the 
+            sequence compared to the consensus via a matrix
+    '''
+    consensus, _ = np.array(findConsensus(arr, '',
+                                          consensus_type='majority_nongap')  )
+    ci_dir = os.path.dirname(utilityFunctions.__file__)
+    matrix_dir = "%s/similarity_matrices" % (ci_dir)
 
-        return:
-        a numpy array stored as new_arr, which is a boolean array comparing the arr to the consensus of it.
-        '''
+    if booleanOrSimilarity == "Boolean":
         bool_array = np.array([])
         bool_arrL = np.empty(dtype=bool, shape=(0, len(consensus)))
         # declares the numpy arrays
         for e in range(1, (len(arr[:,0])+1)):
             # iterates over the rows of the sequences
             z = e-1
-        for i in range(1, (len(arr[0,:])+1)):
-            # iterates over the columns of the sequences
-            x = i-1
-            if arr[z,x] == consensus[x]:
-                # verifies if the current value being iterated is equal to the equivalent value inline with the consensus
-                bool_array = np.append(bool_array, [True], axis=None)
-            else:
-                bool_array = np.append(bool_array, [False], axis=None)
-        bool_arrL = np.vstack([bool_arrL, bool_array])
-        bool_array = np.array([])
+            for i in range(1, (len(arr[0,:])+1)):
+                # iterates over the columns of the sequences
+                x = i-1
+                if arr[z,x] == consensus[x]:
+                    # verifies if the current value being iterated is equal to the
+                    #equivalent value inline with the consensus
+                    bool_array = np.append(bool_array, [True], axis=None)
+                else:
+                    bool_array = np.append(bool_array, [False], axis=None)
+            bool_arrL = np.vstack([bool_arrL, bool_array])
+            bool_array = np.array([])
+
         new_arr = copy.deepcopy(bool_arrL)
         new_arr = bool_arrL.astype(bool)
-        # returns the new boolean array containing the verified alignment to the consensus
+        # returns the new boolean array containing the verified alignment to
+        #the consensus
         return new_arr
     else:
         # generates the consensus
         Sarray = np.array([])
         SarrL = np.empty(dtype=int, shape=(0, len(consensus)))
         # declares the numpy arrays
-        tab = pd.read_csv("roman_work_experience/matrices.txt", sep="\t", index_col=0)
+        tab = pd.read_csv("%s/matrices.txt" % ci_dir, sep="\t", index_col=0)
         if typ == "aa":
           # verifies if the typ is amino acid or nucleotide
           if MatrixName != "B":
             if tab.loc[MatrixName][0] != typ:
               raise RuntimeError("This matrix is not valid")
               # verifies if the matrix is valid
             else:
-              # verifies if the user would like to use the default matrix or their own
-              mat = pd.read_csv(("%s/similarity_matrices/"+MatrixName) % mydir, comment="#", sep="\s+")
+              # verifies if the user would like to use the default matrix or
+              #their own
+              mat = pd.read_csv("%s/%s" % (matrix_dir, MatrixName),
+                                comment="#", sep="\s+")
           elif MatrixName == "B":
-            mat = pd.read_csv("%s/similarity_matrices/BLOSUM62" % mydir, comment="#", sep="\s+")
+            mat = pd.read_csv("%s/BLOSUM62" % (matrix_dir),
+                              comment="#", sep="\s+")
         elif typ == "nt":
             if MatrixName != "B":
                 if tab.loc[MatrixName][0] != typ:
                     raise RuntimeError("This matrix is not valid")
                     # verifies if the matrix is valid
                 else:
-                    # verifies if the user would like to use the default matrix or their own
-                    mat = pd.read_csv(("%s/similarity_matrices/"+MatrixName) % mydir, comment="#", sep="\s+")
+                    # verifies if the user would like to use the default
+                    # matrix or their own
+                    mat = pd.read_csv("%s/%s" % (matrix_dir, MatrixName),
+                                      comment="#", sep="\s+")
             elif MatrixName == "B":
-                mat = pd.read_csv("%s/similarity_matrices/NUC.4.4" % mydir, comment="#", sep="\s+")
+                mat = pd.read_csv("%s/NUC.4.4" % (matrix_dir), comment="#",
+                                  sep="\s+")
         for e in range(1, (len(arr[:,0])+1)):
             # iterates over the rows of the sequences
             z = e-1
             for i in range(1, (len(arr[0,:])+1)):
                 #  iterates over the columns of the sequences
                 x = i-1
-                if not arr[z,x] == "-":
-                      Sarray = np.append(Sarray,[int(mat.loc[arr[z,x],consensus[x]])])
+                if not arr[z, x] == "-":
+                    if arr[z, x] == "U" and typ == 'nt':
+                        thischar = "T"
+                    else:
+                        thischar = arr[z, x]
+                    if consensus[x] == "U" and typ == 'nt':
+                        conschar = "T"
+                    else:
+                        conschar = consensus[x]
+                    score = mat.loc[thischar, conschar]
+                    Sarray = np.append(Sarray, [score])
                 elif arr[z,x] == "-":
                       # sets the value of '-' as 0
                       Sarray = np.append(Sarray, 0)
             SarrL = np.vstack([SarrL, Sarray])
             Sarray = np.array([])
         new_Sarr = copy.deepcopy(SarrL)
         new_Sarr = SarrL.astype(int)
-        # returns the new similarity array containing the verified alignment to the consensus
+        # returns the new similarity array containing the verified alignment
+        # to the consensus
         return new_Sarr
 
 

CIAlign/matrices.txt

@@ -0,0 +1,73 @@
+Matrix	Type
+BLOSUM100	aa
+BLOSUM30	aa
+BLOSUM35	aa
+BLOSUM40	aa
+BLOSUM45	aa
+BLOSUM50	aa
+BLOSUM55	aa
+BLOSUM60	aa
+BLOSUM62	aa
+BLOSUM65	aa
+BLOSUM70	aa
+BLOSUM75	aa
+BLOSUM80	aa
+BLOSUM85	aa
+BLOSUM90	aa
+BLOSUMN	aa
+DAYHOFF	aa
+GONNET	aa
+IDENTITY	aa
+MATCH	aa
+NUC.4.4	nt
+NUC.4.2	nt
+PAM10	aa
+PAM100	aa
+PAM110	aa
+PAM120	aa
+PAM130	aa
+PAM140	aa
+PAM150	aa
+PAM160	aa
+PAM170	aa
+PAM180	aa
+PAM190	aa
+PAM20	aa
+PAM200	aa
+PAM210	aa
+PAM220	aa
+PAM230	aa
+PAM240	aa
+PAM250	aa
+PAM260	aa
+PAM270	aa
+PAM280	aa
+PAM290	aa
+PAM30	aa
+PAM300	aa
+PAM310	aa
+PAM320	aa
+PAM330	aa
+PAM340	aa
+PAM350	aa
+PAM360	aa
+PAM370	aa
+PAM380	aa
+PAM390	aa
+PAM40	aa
+PAM400	aa
+PAM410	aa
+PAM420	aa
+PAM430	aa
+PAM440	aa
+PAM450	aa
+PAM460	aa
+PAM470	aa
+PAM480	aa
+PAM490	aa
+PAM50	aa
+PAM500	aa
+PAM60	aa
+PAM70	aa
+PAM80	aa
+PAM90	aa

CIAlign/miniAlignments.py

@@ -250,7 +250,10 @@ def drawMarkUpLegend(outfile, palette="CBS"):
 def drawMiniAlignment(arr, nams, log, outfile, typ, plot_type='standard',
                       dpi=300, title=None, width=5, height=3, markup=False,
                       markupdict=None, ret=False, orig_nams=[],
-                      keep_numbers=False, force_numbers=False, palette="CBS"):
+                      keep_numbers=False, force_numbers=False, palette="CBS",
+                      sub_matrix_name='B', plot_identity_palette='terrain_r',
+                      plot_similarity_palette='summer_r',
+                      plot_substitution_matrix='B'):
     '''
     Draws a "mini alignment" image showing a small representation of the
     whole alignment so that gaps and poorly aligned regions are visible.
@@ -322,8 +325,15 @@ def drawMiniAlignment(arr, nams, log, outfile, typ, plot_type='standard',
     if plot_type == 'standard':
         arr2, cm = arrNumeric(arr, typ, palette)
     elif plot_type == 'boolean':
-        arr2 = consensusSeq.compareAlignmentConsensus(arr)
-        cm = matplotlib.colormaps['jet']
+        arr2 = consensusSeq.compareAlignmentConsensus(
+            arr, typ=typ, booleanOrSimilarity="Boolean")
+        cm = matplotlib.colormaps[plot_identity_palette]
+    elif plot_type == 'similarity':
+        arr2 = consensusSeq.compareAlignmentConsensus(
+            arr, typ=typ, booleanOrSimilarity="similarity",
+            MatrixName=sub_matrix_name)
+        cm = matplotlib.colormaps[plot_similarity_palette]
+
     # display it on the axis
     a.imshow(arr2, cmap=cm, aspect='auto', interpolation='nearest')
 

CIAlign/runCIAlign.py

@@ -162,6 +162,7 @@ def whichFunctions(args):
             args.plot_output,
             args.plot_markup,
             args.plot_consensus_identity,
+            args.plot_consensus_similarity,
             args.visualise,
             args.all_options]):
         which_functions.append("mini_alignments")
@@ -892,14 +893,33 @@ def runMiniAlignments(args, log, orig_arr, orig_nams, arr, nams,
             print("Plotting identity to consensus")
         outf = "%s_consensus_identity.%s" % (args.outfile_stem,
                                              args.plot_format)
-        miniAlignments.drawMiniAlignment(arr, nams, log,
-                                         outf, typ,
-                                         plot_type='boolean',
-                                         dpi=args.plot_dpi,
-                                         width=args.plot_width,
-                                         height=args.plot_height,
-                                         force_numbers=fn,
-                                         palette=args.palette)
+        miniAlignments.drawMiniAlignment(
+            arr, nams, log,
+            outf, typ,
+            plot_type='boolean',
+            dpi=args.plot_dpi,
+            width=args.plot_width,
+            height=args.plot_height,
+            force_numbers=fn,
+            palette=args.palette,
+            plot_identity_palette=args.plot_identity_palette)
+    if args.plot_consensus_similarity:
+        log.info("Plotting similarity to consensus")
+        if not args.silent:
+            print("Plotting similarity to consensus")
+        outf = "%s_consensus_similarity.%s" % (args.outfile_stem,
+                                             args.plot_format)
+        miniAlignments.drawMiniAlignment(
+            arr, nams, log,
+            outf, typ,
+            plot_type='similarity',
+            dpi=args.plot_dpi,
+            width=args.plot_width,
+            height=args.plot_height,
+            force_numbers=fn,
+            palette=args.palette,
+            plot_similarity_palette=args.plot_similarity_palette,
+            plot_substitution_matrix=args.plot_substitution_matrix)
 
 
 def runConsensus(args, log, orig_arr, orig_nams, arr, nams, removed_seqs):

CIAlign/similarity_matrices/BLOSUM100

@@ -0,0 +1,31 @@
+#  Matrix made by matblas from blosum100_3.iij
+#  * column uses minimum score
+#  BLOSUM Clustered Scoring Matrix in 1/3 Bit Units
+#  Blocks Database = /data/blocks_5.0/blocks.dat
+#  Cluster Percentage: >= 100
+#  Entropy =   1.4516, Expected =  -1.0948
+   A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  B  Z  X  *
+A  8 -3 -4 -5 -2 -2 -3 -1 -4 -4 -4 -2 -3 -5 -2  1 -1 -6 -5 -2 -4 -2 -2 -10 
+R -3 10 -2 -5 -8  0 -2 -6 -1 -7 -6  3 -4 -6 -5 -3 -3 -7 -5 -6 -4 -1 -3 -10 
+N -4 -2 11  1 -5 -1 -2 -2  0 -7 -7 -1 -5 -7 -5  0 -1 -8 -5 -7  5 -2 -3 -10 
+D -5 -5  1 10 -8 -2  2 -4 -3 -8 -8 -3 -8 -8 -5 -2 -4 -10 -7 -8  6  0 -4 -10 
+C -2 -8 -5 -8 14 -7 -9 -7 -8 -3 -5 -8 -4 -4 -8 -3 -3 -7 -6 -3 -7 -8 -5 -10 
+Q -2  0 -1 -2 -7 11  2 -5  1 -6 -5  2 -2 -6 -4 -2 -3 -5 -4 -5 -2  5 -2 -10 
+E -3 -2 -2  2 -9  2 10 -6 -2 -7 -7  0 -5 -8 -4 -2 -3 -8 -7 -5  0  7 -3 -10 
+G -1 -6 -2 -4 -7 -5 -6  9 -6 -9 -8 -5 -7 -8 -6 -2 -5 -7 -8 -8 -3 -5 -4 -10 
+H -4 -1  0 -3 -8  1 -2 -6 13 -7 -6 -3 -5 -4 -5 -3 -4 -5  1 -7 -2 -1 -4 -10 
+I -4 -7 -7 -8 -3 -6 -7 -9 -7  8  2 -6  1 -2 -7 -5 -3 -6 -4  4 -8 -7 -3 -10 
+L -4 -6 -7 -8 -5 -5 -7 -8 -6  2  8 -6  3  0 -7 -6 -4 -5 -4  0 -8 -6 -3 -10 
+K -2  3 -1 -3 -8  2  0 -5 -3 -6 -6 10 -4 -6 -3 -2 -3 -8 -5 -5 -2  0 -3 -10 
+M -3 -4 -5 -8 -4 -2 -5 -7 -5  1  3 -4 12 -1 -5 -4 -2 -4 -5  0 -7 -4 -3 -10 
+F -5 -6 -7 -8 -4 -6 -8 -8 -4 -2  0 -6 -1 11 -7 -5 -5  0  4 -3 -7 -7 -4 -10 
+P -2 -5 -5 -5 -8 -4 -4 -6 -5 -7 -7 -3 -5 -7 12 -3 -4 -8 -7 -6 -5 -4 -4 -10 
+S  1 -3  0 -2 -3 -2 -2 -2 -3 -5 -6 -2 -4 -5 -3  9  2 -7 -5 -4 -1 -2 -2 -10 
+T -1 -3 -1 -4 -3 -3 -3 -5 -4 -3 -4 -3 -2 -5 -4  2  9 -7 -5 -1 -2 -3 -2 -10 
+W -6 -7 -8 -10 -7 -5 -8 -7 -5 -6 -5 -8 -4  0 -8 -7 -7 17  2 -5 -9 -7 -6 -10 
+Y -5 -5 -5 -7 -6 -4 -7 -8  1 -4 -4 -5 -5  4 -7 -5 -5  2 12 -5 -6 -6 -4 -10 
+V -2 -6 -7 -8 -3 -5 -5 -8 -7  4  0 -5  0 -3 -6 -4 -1 -5 -5  8 -7 -5 -3 -10 
+B -4 -4  5  6 -7 -2  0 -3 -2 -8 -8 -2 -7 -7 -5 -1 -2 -9 -6 -7  6  0 -4 -10 
+Z -2 -1 -2  0 -8  5  7 -5 -1 -7 -6  0 -4 -7 -4 -2 -3 -7 -6 -5  0  6 -2 -10 
+X -2 -3 -3 -4 -5 -2 -3 -4 -4 -3 -3 -3 -3 -4 -4 -2 -2 -6 -4 -3 -4 -2 -3 -10 
+* -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 -10  1