too-many-cells 2.1.1.0 → 3.0.1.0

raw patch · 34 files changed

+2536/−648 lines, 34 filesdep +hvegadep +hvega-themedep +optparse-applicativedep −optparse-generic

Dependencies added: hvega, hvega-theme, optparse-applicative, ploterific

Dependencies removed: optparse-generic

Files

• app/Main.hs +20/−13
• src/TooManyCells/Differential/Differential.hs +7/−6
• src/TooManyCells/Differential/Types.hs +1/−1
• src/TooManyCells/Diversity/Load.hs +2/−2
• src/TooManyCells/File/Types.hs +2/−0
• src/TooManyCells/MakeTree/Load.hs +1/−1
• src/TooManyCells/MakeTree/Plot.hs +2/−2
• src/TooManyCells/MakeTree/Utility.hs +11/−0
• src/TooManyCells/Matrix/Load.hs +46/−31
• src/TooManyCells/Matrix/Preprocess.hs +23/−0
• src/TooManyCells/Matrix/Types.hs +4/−2
• src/TooManyCells/Matrix/Utility.hs +75/−13
• src/TooManyCells/Motifs/FindMotif.hs +2/−2
• src/TooManyCells/Peaks/ClusterPeaks.hs +5/−3
• src/TooManyCells/Program/Classify.hs +13/−11
• src/TooManyCells/Program/Differential.hs +26/−28
• src/TooManyCells/Program/Diversity.hs +13/−12
• src/TooManyCells/Program/Interactive.hs +13/−12
• src/TooManyCells/Program/LoadMatrix.hs +34/−27
• src/TooManyCells/Program/MakeTree.hs +82/−98
• src/TooManyCells/Program/MatrixOutput.hs +6/−5
• src/TooManyCells/Program/Motifs.hs +8/−16
• src/TooManyCells/Program/Options.hs +420/−289
• src/TooManyCells/Program/Paths.hs +13/−19
• src/TooManyCells/Program/Peaks.hs +18/−28
• src/TooManyCells/Program/Spatial.hs +245/−0
• src/TooManyCells/Program/Utility.hs +20/−20
• src/TooManyCells/Spatial/AnnoSpat.hs +108/−0
• src/TooManyCells/Spatial/ProjectionPlot.hs +264/−0
• src/TooManyCells/Spatial/Relationships.hs +248/−0
• src/TooManyCells/Spatial/SummaryPlot.hs +597/−0
• src/TooManyCells/Spatial/Types.hs +42/−0
• src/TooManyCells/Spatial/Utility.hs +146/−0
• too-many-cells.cabal +19/−7

@@ -9,7 +9,7 @@ module Main where  -- Remote-import Options.Generic+import Options.Applicative  -- Local import TooManyCells.Program.Differential@@ -22,19 +22,26 @@ import TooManyCells.Program.Options import TooManyCells.Program.Paths import TooManyCells.Program.Peaks+import TooManyCells.Program.Spatial  main :: IO () main = do-    opts <- getRecord "too-many-cells, Gregory W. Schwartz.\-                      \ Clusters and analyzes single cell data."+  let opts = info (sub <**> helper)+                ( fullDesc+              <> progDesc "Clusters and analyzes single cell data."+              <> header "too-many-cells, Gregory W. Schwartz" ) -    case opts of-        MakeTree{}     -> makeTreeMain opts-        Interactive{}  -> interactiveMain opts-        Differential{} -> differentialMain opts-        Diversity{}    -> diversityMain opts-        Paths{}        -> pathsMain opts-        Classify{}     -> classifyMain opts-        Peaks{}        -> peaksMain opts-        Motifs{}       -> motifsMain opts-        MatrixOutput{} -> matrixOutputMain opts+  config <- customExecParser (prefs (showHelpOnEmpty <> showHelpOnError)) opts++  let mainEntry (MakeTreeCommand _)     = makeTreeMain config+      mainEntry (InteractiveCommand _)  = interactiveMain config+      mainEntry (DifferentialCommand _) = differentialMain config+      mainEntry (DiversityCommand _)    = diversityMain config+      mainEntry (PathsCommand _)        = pathsMain config+      mainEntry (ClassifyCommand _)     = classifyMain config+      mainEntry (PeaksCommand _)        = peaksMain config+      mainEntry (MotifsCommand _)       = motifsMain config+      mainEntry (MatrixOutputCommand _) = matrixOutputMain config+      mainEntry (SpatialCommand _)      = spatialMain config++  mainEntry config

@@ -260,7 +260,7 @@             -> B.ByteString getDEString xs = header <> "\n" <> body   where-    header = "feature,log2FC,pVal,FDR"+    header = "feature,log2FC,pVal,qVal"  -- edgeR calls q-values FDR     body   = CSV.encode            . fmap ( L.over L._1 Diff.unName                   . L.over L._3 Diff.unPValue@@ -275,10 +275,10 @@   -> B.ByteString getAllDEStringKruskalWallis xs = header <> "\n" <> body   where-    header = "node,feature,log2FC,pVal,FDR"+    header = "node,feature,log2FC,pVal,qVal"     body   = CSV.encode-           . fmap ( L.over L._6 (maybe "NA" (showt . Diff.unQValue))-                  . L.over L._5 (maybe "NA" (showt . Diff.unFDR))+           . fmap ( (\(!a, !b, !c, !d, !e, !f) -> (a,b,c,d,f))+                  . L.over L._6 (maybe "NA" (showt . Diff.unQValue))                   . L.over L._4 (maybe "NA" (showt . Diff.unPValue))                   . L.over L._3 Diff.unLog2Diff                   . L.over L._2 unFeature@@ -292,9 +292,10 @@   -> B.ByteString getDEStringKruskalWallis xs = header <> "\n" <> body   where-    header = "feature,log2FC,pVal,FDR,qVal"+    header = "feature,log2FC,pVal,qVal"     body   = CSV.encode-           . fmap ( L.over L._5 (maybe "NA" (showt . Diff.unQValue))+           . fmap ( (\(!a, !b, !c, !d, !e) -> (a,b,c,e))+                  . L.over L._5 (maybe "NA" (showt . Diff.unQValue))                   . L.over L._4 (maybe "NA" (showt . Diff.unFDR))                   . L.over L._3 (maybe "NA" (showt . Diff.unPValue))                   . L.over L._2 Diff.unLog2Diff

@@ -17,7 +17,7 @@  -- Basic newtype TopN = TopN { unTopN :: Int }-newtype NoEdger = NoEdger { unNoEdger :: Bool }+newtype Edger = Edger { unEdger :: Bool } newtype DiffNodes = DiffNodes {unDiffNodes :: ([G.Node], [G.Node])} newtype DiffLabels =   DiffLabels { unDiffLabels :: (Maybe (Set.Set Label), Maybe (Set.Set Label)) }

@@ -27,7 +27,7 @@ import qualified Control.Lens as L import qualified Data.Aeson as A import qualified Data.ByteString.Lazy.Char8 as B-import qualified Data.ByteString.Streaming.Char8 as BS+import qualified Streaming.ByteString.Char8 as BS import qualified Data.Foldable as F import qualified Data.Map.Strict as Map import qualified Data.Sequence as Seq@@ -77,7 +77,7 @@                     . S.toList_                     . S.map getCols                     . S.decodeByName-                    $ (contents :: BS.ByteString (ExceptT S.CsvParseException Managed) ())+                    $ (contents :: BS.ByteStream (ExceptT S.CsvParseException Managed) ())          return population 

@@ -30,6 +30,7 @@ import qualified Data.Sequence as Seq import qualified Data.Sparse.Common as S import qualified Numeric.LinearAlgebra as H+import qualified Turtle as TU  -- Local @@ -45,6 +46,7 @@     { unPriorPath :: FilePath     } deriving (Eq,Ord,Read,Show) newtype OutputDirectory  = OutputDirectory { unOutputDirectory :: FilePath }+newtype TempPath  = TempPath { unTempPath :: TU.FilePath }  -- Advanced data MatrixFileFolder = MatrixFile FilePath | MatrixFolder FilePath

@@ -29,7 +29,7 @@ import qualified Control.Lens as L import qualified Data.Aeson as A import qualified Data.ByteString.Lazy.Char8 as B-import qualified Data.ByteString.Streaming.Char8 as BS+import qualified Streaming.ByteString.Char8 as BS import qualified Data.Csv as CSV import qualified Data.Foldable as F import qualified Data.Map.Strict as Map

@@ -70,7 +70,7 @@                -> Maybe (Double, Double, String) getProjections (ProjectionMap pm) (CellInfo { _barcode = !b }, Cluster !c) = do     p <- Map.lookup b pm-    return $ (unX . fst $ p, unY . snd $ p, show c)+    return $ (unX . fst . snd $ p, unY . snd . snd $ p, show c)  -- | Get the projections and label with colors for a cell. getProjectionsColor :: LabelMap@@ -90,7 +90,7 @@       l = T.unpack . unLabel . barcodeToLabel . getId $ cell       c = barcodeToColor . getId $ cell -  return $ (unX . fst $ p, unY . snd $ p, l, c)+  return $ (unX . fst . snd $ p, unY . snd . snd $ p, l, c)  -- | Plot clusters. plotClustersR :: String -> ProjectionMap -> [(CellInfo, [Cluster])] -> R s ()

@@ -10,6 +10,7 @@ module TooManyCells.MakeTree.Utility     ( updateTreeRowIndex     , updateTreeRowBool+    , projectionMapToCoordinateMap     ) where  -- Remote@@ -18,6 +19,9 @@ import Data.Tree (Tree) import qualified Control.Lens as L import qualified Data.HashMap.Strict as HMap+import qualified Data.List as List+import qualified Data.Map.Strict as Map+import qualified Data.Sparse.Common as S import qualified Data.Vector as V  -- Local@@ -49,3 +53,10 @@ updateTreeRowBool (UpdateTreeRowsFlag True) Nothing tree = tree updateTreeRowBool (UpdateTreeRowsFlag True) (Just sc) tree =   updateTreeRowIndex sc tree++-- | Convert ProjectionMap to CoordinateMap for birch-beer.+projectionMapToCoordinateMap :: ProjectionMap -> CoordinateMap+projectionMapToCoordinateMap = CoordinateMap+  . Map.map (L.over L._2 (\(X !x, Y !y) -> S.fromListDenseSV 2 [x, y]))+  . Map.mapKeys (Id . unCell)+  . unProjectionMap

@@ -44,7 +44,7 @@ import qualified Control.Lens as L import qualified Data.ByteString.Char8 as B import qualified Data.ByteString.Lazy.Char8 as BL-import qualified Data.ByteString.Streaming.Char8 as BS+import qualified Streaming.ByteString.Char8 as BS import qualified Data.Csv as CSV import qualified Data.Foldable as F import qualified Data.HashMap.Strict as HMap@@ -241,7 +241,7 @@                 . (S.store gS)                 . (S.store mS)                 . S.decode S.NoHeader-                $ (contents :: BS.ByteString (ExceptT S.CsvParseException Managed) ())+                $ (contents :: BS.ByteStream (ExceptT S.CsvParseException Managed) ())          let finalMat = MatObsRow . HS.sparsifySM . HS.fromColsL $ m -- We want observations as rows @@ -286,7 +286,10 @@       stream = preprocessStream              . fmap (T.splitOn "\t")              . (\x -> maybe x (flip filterFragments x) excludeFragments)  -- Filter out unwanted fragments by name match-             . Turtle.mfilter (not . T.null)+             . Turtle.mfilter (\ x+                              -> (not . T.null $ x)+                              && ((fmap fst . T.uncons $ x) /= Just '#')  -- No commented or empty lines+                              )              . fmap Turtle.lineToText              . (maybe id filterBlacklist blacklist)              . Turtle.inproc "bedtools" ["sort", "-i", "stdin"]@@ -358,7 +361,10 @@       stream = preprocessStream              . fmap (T.splitOn "\t")              . (\x -> maybe x (flip filterFragments x) excludeFragments)  -- Filter out unwanted fragments by name match-             . Turtle.mfilter (not . T.null)+             . Turtle.mfilter (\ x+                              -> (not . T.null $ x)+                              && ((fmap fst . T.uncons $ x) /= Just '#')  -- No commented or empty lines+                              )              . fmap Turtle.lineToText              . (maybe id filterBlacklist blacklist)              . Turtle.inproc "bedtools" ["sort", "-i", "stdin"]@@ -387,22 +393,22 @@                   , _colNames = fmap Feature features                   } --- | Load a projection file to get a vector of projections.-loadProjectionFile :: ProjectionFile -> IO (Vector (X, Y))-loadProjectionFile =-    fmap (\ x -> either-                                error-                                (fmap getProjection . V.drop 1)-                                (CSV.decode CSV.NoHeader x :: Either String (Vector [T.Text]))-         )-        . BL.readFile-        . unProjectionFile-  where-    getProjection [_, x, y] = ( X . either error fst . T.double $ x-                              , Y . either error fst . T.double $ y-                              )-    getProjection xs        =-        error $ "\nUnrecognized projection row: " <> show xs+-- -- | Load a projection file to get a vector of projections. Legacy.+-- loadProjectionFile :: ProjectionFile -> IO (Vector (X, Y))+-- loadProjectionFile =+--     fmap (\ x -> either+--                                 error+--                                 (fmap getProjection . V.drop 1)+--                                 (CSV.decode CSV.NoHeader x :: Either String (Vector [T.Text]))+--          )+--         . BL.readFile+--         . unProjectionFile+--   where+--     getProjection [_, x, y] = ( X . either error fst . T.double $ x+--                               , Y . either error fst . T.double $ y+--                               )+--     getProjection xs        =+--         error $ "\nUnrecognized projection row: " <> show xs  -- | Load a projection file as a map. loadProjectionMap :: ProjectionFile -> IO ProjectionMap@@ -410,21 +416,30 @@     fmap (\ x -> ProjectionMap                . Map.fromList                . V.toList-               . either-                  error-                  (fmap getProjection . V.drop 1)-               $ (CSV.decode CSV.NoHeader x :: Either String (Vector [T.Text]))+               . either error (fmap getProjection . snd)+               $ ( CSV.decodeByName x+                :: Either String (CSV.Header, V.Vector (Map.Map T.Text T.Text))+                 )          )         . BL.readFile         . unProjectionFile   where-    getProjection [b, x, y] = ( Cell b-                              , ( X . either error fst . T.double $ x-                                , Y . either error fst . T.double $ y-                                )-                              )-    getProjection xs        =-        error $ "\nUnrecognized projection row: " <> show xs+    getProjection m = ( Cell $ lookupErr "item" m+                      , ( fmap Sample $ Map.lookup "sample" m+                        , ( X . either error fst . T.double . lookupErr "x" $ m+                          , Y . either error fst . T.double . lookupErr "y" $ m+                          )+                        )+                      )+    lookupErr col m = Map.findWithDefault+                        ( error+                        $ "\nUnrecognized projection row "+                       <> show col+                       <> " in row: "+                       <> show m+                        )+                        col+                        m  -- | Decide to use first two values as the projection. toPoint :: HS.SpVector Double -> (X, Y)

@@ -17,6 +17,8 @@     , logCPMSparseMat     , uqScaleSparseMat     , medScaleSparseMat+    , minMaxNormSparseMat+    , transposeSparseMat     , quantileScaleSparseMat     , centerScaleSparseCell     , tfidfScaleSparseMat@@ -202,6 +204,20 @@ tfidfScaleSparseMat :: MatObsRow -> MatObsRow tfidfScaleSparseMat = MatObsRow . unB2 . b1ToB2 . B1 . unMatObsRow +-- | Scale matrix rows based on min-max normalization.+minMaxNormSparseMat :: MatObsRow -> MatObsRow+minMaxNormSparseMat (MatObsRow mat) =+  MatObsRow+    . S.sparsifySM+    . S.fromRowsL+    . fmap (\x -> fromMaybe x $ minMaxNormSV x)+    . S.toRowsL+    $ mat++-- | Transpose matrix (helper to apply normalizations to columns instead, for instance).+transposeSparseMat :: MatObsRow -> MatObsRow+transposeSparseMat (MatObsRow mat) = MatObsRow . S.transposeSM $ mat+ -- | Scale a cell by the library size. scaleDenseCell :: H.Vector H.R -> H.Vector H.R scaleDenseCell xs = H.cmap (/ total) xs@@ -277,6 +293,13 @@         . fmap snd         . S.toListSV         $ xs++-- | Min-max normalization of sparse vector.+minMaxNormSV :: S.SpVector Double -> Maybe (S.SpVector Double)+minMaxNormSV xs = do+  (mi, ma) <-+    L.sequenceOf L.both $ Fold.fold ((,) <$> Fold.minimum <*> Fold.maximum) xs+  pure $ fmap (\x -> (x - mi) / (ma - mi)) xs  -- | Filter a matrix to remove low count cells and features. filterDenseMat :: FilterThresholds -> SingleCells -> SingleCells

@@ -126,8 +126,8 @@     { unY :: Double     } deriving (Eq,Ord,Read,Show,Generic) deriving newtype (Num) deriving anyclass (A.ToJSON,A.FromJSON) newtype ProjectionMap = ProjectionMap-  { unProjectionMap :: Map.Map Cell (X, Y)-  }+  { unProjectionMap :: Map.Map Cell (Maybe Sample, (X, Y))+  } deriving (Read, Show) newtype RMat s          = RMat { unRMat :: R.SomeSEXP s } newtype RMatObsRow s    = RMatObsRow { unRMatObsRow :: R.SomeSEXP s } newtype RMatFeatRow s   = RMatFeatRow { unRMatFeatRow :: R.SomeSEXP s }@@ -355,6 +355,8 @@               | TotalNorm               | LogCPMNorm Double               | QuantileNorm+              | MinMaxNorm+              | TransposeNorm               | NoneNorm                 deriving (Read, Show, Eq) 

@@ -13,6 +13,7 @@     ( matToRMat     , scToRMat     , sparseToHMat+    , sparseMatToSparseRMat     , hToSparseMat     , matToHMat     , matToSpMat@@ -29,6 +30,7 @@     , extractCellSc     , extractCellsSc     , removeCellsSc+    , subsetCellsSc     , extractCellV     , aggSc     ) where@@ -36,6 +38,7 @@ -- Remote import BirchBeer.Types import Codec.Compression.GZip (compress)+import Control.Monad (liftM2) import Control.Monad.Managed (runManaged) import Control.Monad.State (MonadState (..), State (..), evalState, execState, modify) import Control.Monad.IO.Class (MonadIO)@@ -44,7 +47,7 @@ import Data.Function (on) import Data.Hashable (Hashable) import Data.List (maximumBy, isInfixOf, foldl')-import Data.Maybe (fromMaybe)+import Data.Maybe (fromMaybe, catMaybes) import Data.Matrix.MatrixMarket (Matrix(RMatrix, IntMatrix), Structure (..), writeMatrix') import Data.Scientific (toRealFloat, Scientific) import Data.Streaming.Zlib (WindowBits)@@ -57,9 +60,10 @@ import qualified Data.Attoparsec.Text as A import qualified Data.ByteString.Char8 as B import qualified Data.ByteString.Lazy.Char8 as BL-import qualified Data.ByteString.Streaming.Char8 as BS+import qualified Streaming.ByteString.Char8 as BS import qualified Data.Clustering.Hierarchical as HC import qualified Data.Graph.Inductive as G+import qualified Data.HashMap.Lazy as HMap import qualified Data.HashSet as HSet import qualified Data.IntMap.Strict as IMap import qualified Data.IntervalMap.Strict as IntervalMap@@ -99,6 +103,35 @@     mat <- [r| as.matrix(fromJSON(mString_hs)) |]     return . RMatObsRow $ mat +-- | Convert a mat to an RMatrix.+sparseMatToSparseRMat :: SingleCells -> R s (RMatObsRow s)+sparseMatToSparseRMat sc = do+    [r| library(Matrix) |]++    let rowNamesR = fmap (T.unpack . unCell) . V.toList . L.view rowNames $ sc+        colNamesR = fmap (T.unpack . unFeature) . V.toList . L.view colNames $ sc+        idxs = S.toListSM . unMatObsRow . L.view matrix $ sc+        d :: [Double]+        d = (\(!m, !n) -> fmap fromIntegral [m, n])+          . S.dimSM+          . unMatObsRow+          . L.view matrix+          $ sc+        is :: [Double]+        is = fmap (fromIntegral . (\x -> x + 1) . L.view L._1) idxs+        js :: [Double]+        js = fmap (fromIntegral . (\x -> x + 1) . L.view L._2) idxs+        xs :: [Double]+        xs = fmap (L.view L._3) idxs++    namedMat <- [r| mat = sparseMatrix(i = is_hs, j = js_hs, x = xs_hs, dims = d_hs)+                    rownames(mat) = rowNamesR_hs+                    colnames(mat) = colNamesR_hs+                    mat+                |]++    return . RMatObsRow $ namedMat+ -- | Convert a sc structure to an RMatrix. scToRMat :: SingleCells -> R s (RMatObsRow s) scToRMat sc = do@@ -220,7 +253,7 @@ printDenseMatrixLike :: (MatrixLike a, Monad m)                      => MatrixTranspose                      -> a-                     -> BS.ByteString m ()+                     -> BS.ByteStream m () printDenseMatrixLike (MatrixTranspose mt) mat =   Stream.encode (Just $ Stream.header header)     . Stream.zipWith (:) rowN@@ -295,7 +328,7 @@ {-# INLINABLE hashNub #-}  -- | Keep decompressing a compressed bytestream until exhaused.-decompressStreamAll :: (MonadIO m) => WindowBits -> BS.ByteString m r -> BS.ByteString m r+decompressStreamAll :: (MonadIO m) => WindowBits -> BS.ByteStream m r -> BS.ByteStream m r decompressStreamAll w bs = S.decompress' w bs >>= go   where     go (Left bs) = S.decompress' w bs >>= go@@ -341,6 +374,14 @@ extractCellV cell sc =   flip S.extractRow (getCellIdx sc cell) . unMatObsRow . L.view matrix $ sc +-- | Get a single cell vector from a SingleCells type using hashes.+extractCellHashV ::+  HMap.HashMap T.Text Int -> SingleCells -> Cell -> Maybe (S.SpVector Double)+extractCellHashV idxMap sc (Cell cell) =+  S.extractRow (unMatObsRow . L.view matrix $ sc) <$> idx+  where+    idx = HMap.lookup cell idxMap+ -- | Remove single cells from a SingleCells type. removeCellsSc :: [Cell] -> SingleCells -> SingleCells removeCellsSc cells sc = L.set rowNames newRowNames@@ -348,16 +389,37 @@                      $ sc   where     blacklistCells = Set.fromList cells-    blacklistIdxs = Set.fromList . fmap (getCellIdx sc) $ cells-    whitelistIdxs = filter-                      (not . flip Set.member blacklistIdxs)-                      [0..(V.length $ L.view rowNames sc) - 1]-    newRowNames =-      V.filter (not . flip Set.member blacklistCells) . L.view rowNames $ sc+    isValidV = fmap (not . flip Set.member blacklistCells) . L.view rowNames $ sc+    newRowNames = fmap snd . V.filter fst . V.zip isValidV . L.view rowNames $ sc     newMatrix = MatObsRow-              . S.fromRowsL-              . fmap (\x -> flip S.extractRow x . unMatObsRow . L.view matrix $ sc)-              $ whitelistIdxs+              . S.fromRowsV+              . fmap snd+              . V.filter fst+              . V.zip isValidV+              . V.fromList+              . S.toRowsL+              . unMatObsRow+              . L.view matrix+              $ sc++-- | Get ordered subset of single cells from a SingleCells type (if they exist).+subsetCellsSc :: [Cell] -> SingleCells -> SingleCells+subsetCellsSc cells sc = L.set rowNames newRowNames+                       . L.set matrix newMatrix+                       $ sc+  where+    idxMap = HMap.fromList+           . flip zip [0..]+           . fmap unCell+           . V.toList+           . L.view rowNames+           $ sc+    newRowNames = V.fromList+                . catMaybes+                . zipWith (liftM2 const) (fmap Just cells)+                $ validRows+    newMatrix = MatObsRow . S.fromRowsL . catMaybes $ validRows+    validRows = fmap (extractCellHashV idxMap sc) cells  -- | Get the index of a cell in a SingleCells type. getCellIdx :: SingleCells -> Cell -> Int

@@ -180,7 +180,7 @@                     (TU.format TU.fp tmpFasta)                     (TU.format TU.fp . unOutputPath $ outPath) -  TU.stdout . TU.inshell cmd $ mempty+  TU.stderr . TU.inshell cmd $ mempty  getMotifGenome :: DiffFile                -> Maybe BackgroundDiffFile@@ -218,4 +218,4 @@                     (unGenomeFile gf)                     (TU.format TU.fp . unOutputPath $ outPath) -  TU.stdout . TU.inshell cmd $ mempty+  TU.stderr . TU.inshell cmd $ mempty

@@ -39,7 +39,7 @@ import qualified Control.Foldl as Fold import qualified Control.Lens as L import qualified Data.ByteString.Lazy.Char8 as B-import qualified Data.ByteString.Streaming.Char8 as BS+import qualified Streaming.ByteString.Char8 as BS import qualified Data.HashMap.Strict as HMap import qualified Data.IntMap.Strict as IMap import qualified Data.IntSet as ISet@@ -123,7 +123,7 @@    _ <- MaybeT . fmap Just $ externalCompressedFragmentsSort file   _ <- MaybeT . fmap Just . callCommand $ TP.printf gcc file gf scale output-  _ <- MaybeT . fmap Just $ externalFragmentsSort output+  -- _ <- MaybeT . fmap Just $ externalFragmentsSort output  -- Removing as sorting ruins track line   _ <- MaybeT . fmap Just . callCommand $ TP.printf "bedGraphToBigWig %s %s %s" output gf wigOut   return () @@ -244,7 +244,9 @@           $ cr       cellClusterMap :: HMap.HashMap T.Text [Int]       cellClusterMap = HMap.fromList clusterAssocList-      processedStream = S.map (T.splitOn "\t" . T.decodeUtf8)+      processedStream = S.map (T.splitOn "\t")+                      . S.filter ((/=) (Just '#') . fmap fst . T.uncons)  -- No commented+                      . S.map T.decodeUtf8                       . S.mapped BS.toStrict                       . BS.lines                       . decompressStreamAll (WindowBits 31) -- For gunzip

@@ -6,14 +6,15 @@  {-# LANGUAGE BangPatterns #-} {-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE DuplicateRecordFields #-}  module TooManyCells.Program.Classify where  -- Remote import BirchBeer.Types (Label (..)) import Control.Monad (unless, when)+import Data.Bool (bool) import Data.Maybe (fromMaybe)-import Options.Generic import Safe (headMay) import Text.Read (readMaybe) import TextShow (showt)@@ -31,29 +32,29 @@ import TooManyCells.Program.Options  -- | Classify path.-classifyMain :: Options -> IO ()-classifyMain opts = do+classifyMain :: Subcommand -> IO ()+classifyMain sub@(ClassifyCommand opts) = do   let readOrErr err = fromMaybe (error err) . readMaybe-      refFiles' = unHelpful-                . referenceFile-                $ opts-      singleRefMatFlag' = unHelpful . singleReferenceMatrix $ opts+      refFiles' = referenceFile opts+      singleRefMatFlag' = singleReferenceMatrix $ opts+      skipAggregation' = skipAggregation $ opts       getRef x = AggReferenceMat-               . aggSc+               . bool aggSc AggSingleCells skipAggregation'                . fst                . fromMaybe (error $ "Could not load file in required field --matrix-path:" <> show x)-             <$> loadAllSSM (opts { matrixPath = Helpful [x] })+             <$> loadAllSSM sub ((\a -> a { matrixPath = [x] }) . (loadMatrixOptions :: Classify -> LoadMatrixOptions) $ opts)       getRefs xs = case (xs, singleRefMatFlag') of                     ([], _) -> error "No --matrix-path specified"                     ([x], True) -> fmap (AggReferenceMat . AggSingleCells)                                  . extractCellsSc                                  . fst                                  . fromMaybe (error $ "Could not load file in required field --matrix-path:" <> show x)-                               <$> loadAllSSM (opts { matrixPath = Helpful [x]})+                               <$> loadAllSSM sub ((\a -> a { matrixPath = [x] }) . (loadMatrixOptions :: Classify -> LoadMatrixOptions) $ opts)                     (xs, True) -> error "Cannot use --single-reference-matrix with more than one --matrix-path."                     (xs, False) -> mapM getRef refFiles' -  sc <- maybe (error "Requires --matrix-path") fst <$> loadAllSSM opts+  sc <- maybe (error "Requires --matrix-path") fst+    <$> loadAllSSM sub ((loadMatrixOptions :: Classify -> LoadMatrixOptions) opts)   refs <- getRefs refFiles'    let classified = classifyCells sc refs@@ -62,3 +63,4 @@    putStrLn "item,label,score"   mapM_ outputRow classified+classifyMain _ = error "Wrong path in classify, contact Gregory Schwartz for this error."

@@ -12,6 +12,7 @@ {-# LANGUAGE PackageImports    #-} {-# LANGUAGE TypeOperators     #-} {-# LANGUAGE TupleSections     #-}+{-# LANGUAGE DuplicateRecordFields #-}  module TooManyCells.Program.Differential where @@ -26,7 +27,6 @@ import Data.Monoid ((<>)) import Language.R as R import Math.Clustering.Hierarchical.Spectral.Types (getClusterItemsDend, EigenGroup (..))-import Options.Generic import System.IO (hPutStrLn, stderr) import Text.Read (readMaybe) import qualified Control.Lens as L@@ -44,36 +44,33 @@ import TooManyCells.MakeTree.Utility import TooManyCells.File.Types import TooManyCells.Matrix.Utility-import TooManyCells.Program.Options import TooManyCells.Program.LoadMatrix import TooManyCells.MakeTree.Load  -- | Differential path.-differentialMain :: Options -> IO ()-differentialMain opts = do+differentialMain :: Subcommand -> IO ()+differentialMain sub@(DifferentialCommand opts) = do     let readOrErr err = fromMaybe (error err) . readMaybe-        delimiter'     =-            Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts+        delimiter'        = Delimiter+                          . (delimiter :: LoadMatrixOptions -> Char)+                          . (loadMatrixOptions :: Differential -> LoadMatrixOptions)+                          $ opts         labelsFile' =-            fmap LabelFile . unHelpful . labelsFile $ opts+            fmap LabelFile . (labelsFile :: Differential -> Maybe String) $ opts         nodes'    =-          DiffNodes . readOrErr "Cannot read --nodes." . unHelpful . nodes $ opts-        prior'    = PriorPath-                  . fromMaybe (error "\nRequires a previous run to get the graph.")-                  . unHelpful-                  . prior-                  $ opts-        topN'     = TopN . fromMaybe 100 . unHelpful . topN $ opts-        features'    = fmap Feature . unHelpful . features $ opts-        aggregate' = Aggregate . unHelpful . aggregate $ opts-        separateNodes' = SeparateNodes . unHelpful . plotSeparateNodes $ opts-        separateLabels' = SeparateLabels . unHelpful . plotSeparateLabels $ opts-        violinFlag' = ViolinFlag . unHelpful . plotViolin $ opts-        noOutlierFlag' = NoOutlierFlag . unHelpful . plotNoOutlier $ opts-        updateTreeRows' = UpdateTreeRowsFlag . not . unHelpful . noUpdateTreeRows $ opts-        noEdger' = NoEdger . unHelpful . noEdger $ opts-        subsampleGroups' = fmap Subsample . unHelpful . subsampleGroups $ opts-        seed' = Seed . fromMaybe 0 . unHelpful . seed $ opts+          DiffNodes . readOrErr "Cannot read --nodes." . nodes $ opts+        prior'    = PriorPath . (prior :: Differential -> String) $ opts+        topN'     = TopN . (topN :: Differential -> Int) $ opts+        features'    = fmap Feature . features $ opts+        aggregate' = Aggregate . aggregate $ opts+        separateNodes' = SeparateNodes . plotSeparateNodes $ opts+        separateLabels' = SeparateLabels . plotSeparateLabels $ opts+        violinFlag' = ViolinFlag . plotViolin $ opts+        noOutlierFlag' = NoOutlierFlag . plotNoOutlier $ opts+        updateTreeRows' = UpdateTreeRowsFlag . not . (noUpdateTreeRows :: Differential -> Bool) $ opts+        edger' = Edger . edger $ opts+        subsampleGroups' = fmap Subsample . subsampleGroups $ opts+        seed' = Seed . seed $ opts         labels'   = fmap ( DiffLabels                          . L.over L.both ( (\x -> bool (Just x) Nothing . Set.null $ x)                                          . Set.fromList@@ -81,16 +78,16 @@                                          )                          . readOrErr "Cannot read --labels."                          )-                  . unHelpful                   . labels                   $ opts         (combined1, combined2) = combineNodesLabels nodes' labels'-        plotOutputR = fromMaybe "out.pdf" . unHelpful . plotOutput $ opts+        plotOutputR = plotOutput opts      when (isNothing labelsFile' && isJust labels') $       hPutStrLn stderr "Warning: labels requested with no label file, ignoring labels..." -    scRes <- loadAllSSM opts+    scRes <- loadAllSSM sub+           $ (loadMatrixOptions :: Differential -> LoadMatrixOptions) opts     let processedSc = fmap fst scRes         customLabelMap = join . fmap snd $ scRes @@ -128,7 +125,7 @@             (DiffNodes ([], _)) -> error "Need other nodes to compare with. If every node should be compared to all other nodes using Mann-Whitney U, use \"([], [])\"."             (DiffNodes (_, [])) -> error "Need other nodes to compare with. If every node should be compared to all other nodes using Mann-Whitney U, use \"([], [])\"."             _ -> do-              if unNoEdger noEdger'+              if not $ unEdger edger'                 then do                   res <- H.io                        $ getDEGraphKruskalWallis@@ -199,3 +196,4 @@           [H.r| suppressMessages(ggsave(diffNormPlot_hs[[1]], file = normOutputR_hs)) |]            return ()+differentialMain _ = error "Wrong path in differential, contact Gregory Schwartz for this error."

@@ -12,6 +12,7 @@ {-# LANGUAGE PackageImports    #-} {-# LANGUAGE TypeOperators     #-} {-# LANGUAGE TupleSections     #-}+{-# LANGUAGE DuplicateRecordFields #-}  module TooManyCells.Program.Diversity where @@ -24,7 +25,6 @@ import Data.Monoid ((<>)) import Language.R as R import Math.Clustering.Hierarchical.Spectral.Types (getClusterItemsDend, EigenGroup (..))-import Options.Generic import qualified Data.ByteString.Lazy.Char8 as B import qualified Data.Colour.Palette.BrewerSet as D import qualified Data.Colour.Palette.Harmony as D@@ -45,23 +45,23 @@ import TooManyCells.Diversity.Plot import TooManyCells.Diversity.Types import TooManyCells.File.Types-import TooManyCells.Program.Options+import qualified TooManyCells.Program.Options as Opt  -- | Diversity path.-diversityMain :: Options -> IO ()-diversityMain opts = do+diversityMain :: Subcommand -> IO ()+diversityMain (DiversityCommand opts) = do     let priors'         =-            fmap PriorPath . unHelpful . priors $ opts+            fmap PriorPath . priors $ opts         delimiter'     =-            Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts+            Delimiter . (delimiter :: Opt.Diversity -> Char) $ opts         labelsFile'       =-            fmap LabelFile . unHelpful . labelsFile $ opts+            fmap LabelFile . (labelsFile :: Opt.Diversity -> Maybe String) $ opts         output'         =-            OutputDirectory . fromMaybe "out" . unHelpful . output $ opts-        order'       = Order . fromMaybe 1 . unHelpful . order $ opts-        start'       = Start . fromMaybe 0 . unHelpful . start $ opts-        interval'    = Interval . fromMaybe 1 . unHelpful . interval $ opts-        endMay'      = fmap End . unHelpful . end $ opts+            OutputDirectory . (output :: Opt.Diversity -> String) $ opts+        order'       = Order . (order :: Opt.Diversity -> Double) $ opts+        start'       = Start . start $ opts+        interval'    = Interval . interval $ opts+        endMay'      = fmap End . end $ opts      -- Where to place output files.     FP.createDirectoryIfMissing True . unOutputDirectory $ output'@@ -127,3 +127,4 @@         return ()      return ()+diversityMain _ = error "Wrong path in diversity, contact Gregory Schwartz for this error."

@@ -12,6 +12,7 @@ {-# LANGUAGE PackageImports    #-} {-# LANGUAGE TypeOperators     #-} {-# LANGUAGE TupleSections     #-}+{-# LANGUAGE DuplicateRecordFields #-}  module TooManyCells.Program.Interactive where @@ -27,7 +28,6 @@ import Language.R as R import Math.Clustering.Hierarchical.Spectral.Types (getClusterItemsDend, EigenGroup (..)) import Math.Clustering.Spectral.Sparse (b1ToB2, B1 (..), B2 (..))-import Options.Generic import qualified Control.Lens as L import qualified Data.Aeson as A import qualified Data.ByteString.Lazy.Char8 as B@@ -46,20 +46,20 @@ import TooManyCells.Program.Utility  -- | Interactive tree interface.-interactiveMain :: Options -> IO ()-interactiveMain opts = H.withEmbeddedR defaultConfig $ do+interactiveMain :: Subcommand -> IO ()+interactiveMain sub@(InteractiveCommand opts) = H.withEmbeddedR defaultConfig $ do     let labelsFile'    =-            fmap LabelFile . unHelpful . labelsFile $ opts-        prior'         = maybe (error "\nRequires --prior") PriorPath-                       . unHelpful-                       . prior-                       $ opts+            fmap LabelFile . (labelsFile :: Interactive -> Maybe String) $ opts+        prior'         = PriorPath . (prior :: Interactive -> String) $ opts         updateTreeRows' =-          UpdateTreeRowsFlag . not . unHelpful . noUpdateTreeRows $ opts-        delimiter'     =-            Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts+          UpdateTreeRowsFlag . not . (noUpdateTreeRows :: Interactive -> Bool) $ opts+        delimiter'        = Delimiter+                          . (delimiter :: LoadMatrixOptions -> Char)+                          . (loadMatrixOptions :: Interactive -> LoadMatrixOptions)+                          $ opts -    scRes <- loadAllSSM opts+    scRes <- loadAllSSM sub+           $ (loadMatrixOptions :: Interactive -> LoadMatrixOptions) opts     let mat = fmap fst scRes         customLabelMap = join . fmap snd $ scRes @@ -86,3 +86,4 @@         $ mat      return ()+interactiveMain _ = error "Wrong path in interactive, contact Gregory Schwartz for this error."

@@ -12,6 +12,7 @@ {-# LANGUAGE PackageImports    #-} {-# LANGUAGE TypeOperators     #-} {-# LANGUAGE TupleSections     #-}+{-# LANGUAGE DuplicateRecordFields #-}  module TooManyCells.Program.LoadMatrix where @@ -28,7 +29,6 @@ import Data.Maybe (fromMaybe, isJust, isNothing) import GHC.Conc (getNumCapabilities) import Math.Clustering.Hierarchical.Spectral.Types (getClusterItemsDend, EigenGroup (..))-import Options.Generic import System.IO (hPutStrLn, stderr) import Text.Read (readMaybe) import qualified Control.Lens as L@@ -53,8 +53,8 @@ import Control.Lens  -- | Load the single cell matrix, post-whitelist-filtered cells.-loadSSM :: Options -> FilePath -> IO SingleCells-loadSSM opts matrixPath' = do+loadSSM :: Subcommand -> LoadMatrixOptions -> FilePath -> IO SingleCells+loadSSM sub opts matrixPath' = do   fileExist      <- FP.doesFileExist matrixPath'   directoryExist <- FP.doesDirectoryExist matrixPath'   compressedFileExist <- FP.doesFileExist $ matrixPath' FP.</> "matrix.mtx.gz"@@ -68,22 +68,19 @@                   $ matrixPath'              FP.</> (bool "barcodes.tsv" "barcodes.tsv.gz" compressedFileExist)       delimiter'      =-          Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts-      transpose'      = TransposeFlag . unHelpful . matrixTranspose $ opts-      featureColumn'  =-          FeatureColumn . fromMaybe 1 . unHelpful . featureColumn $ opts-      noBinarizeFlag' = NoBinarizeFlag . unHelpful . noBinarize $ opts-      binWidth' = fmap BinWidth . unHelpful . binwidth $ opts+          Delimiter . (delimiter :: LoadMatrixOptions -> Char) $ opts+      transpose'      = TransposeFlag . matrixTranspose $ opts+      featureColumn'  = FeatureColumn . featureColumn $ opts+      noBinarizeFlag' = NoBinarizeFlag . noBinarize $ opts+      binWidth' = fmap BinWidth . binwidth $ opts       excludeFragments' = fmap (ExcludeFragments . T.pack)-                        . unHelpful                         . excludeMatchFragments                         $ opts       blacklistRegionsFile' = fmap BlacklistRegions-                            . unHelpful                             . blacklistRegionsFile                             $ opts       cellWhitelistFile' =-            fmap CellWhitelistFile . unHelpful . cellWhitelistFile $ opts+            fmap CellWhitelistFile . cellWhitelistFile $ opts    cellWhitelist <- liftIO . sequence $ fmap getCellWhitelist cellWhitelistFile' @@ -147,37 +144,45 @@    fmap (windowSc binWidth' . whiteListFilter cellWhitelist . transposeFunc) unFilteredSc +-- | Ensure proper orientation of matrix based on the number of transposes for+-- normalizations.+ensureProperAxis :: [NormType] -> [NormType]+ensureProperAxis xs = bool (xs <> [TransposeNorm]) xs+                    . (== 0)+                    . flip mod 2+                    . length+                    . filter (== TransposeNorm)+                    $ xs+ -- | Load all single cell matrices.-loadAllSSM :: Options -> IO (Maybe (SingleCells, Maybe LabelMap))-loadAllSSM opts = runMaybeT $ do-  let matrixPaths'       = unHelpful . matrixPath $ opts-      normalizations'    = getNormalization opts-      pca'               = fmap PCADim . unHelpful . pca $ opts-      lsa'               = fmap LSADim . unHelpful . lsa $ opts-      svd'               = fmap SVDDim . unHelpful . svd $ opts-      dropDimensionFlag' = DropDimensionFlag . unHelpful . dropDimension $ opts-      binarizeFlag'      = BinarizeFlag . unHelpful . binarize $ opts-      binWidth'          = fmap BinWidth . unHelpful . binwidth $ opts+loadAllSSM :: Subcommand -> LoadMatrixOptions -> IO (Maybe (SingleCells, Maybe LabelMap))+loadAllSSM sub opts = runMaybeT $ do+  let matrixPaths'       = matrixPath $ opts+      normalizations'    =+        ensureProperAxis $ getNormalization sub (normalization opts)+      pca'               = fmap PCADim . pca $ opts+      lsa'               = fmap LSADim . lsa $ opts+      svd'               = fmap SVDDim . svd $ opts+      dropDimensionFlag' = DropDimensionFlag . dropDimension $ opts+      binarizeFlag'      = BinarizeFlag . binarize $ opts+      binWidth'          = fmap BinWidth . binwidth $ opts       shiftPositiveFlag' =-        ShiftPositiveFlag . unHelpful . shiftPositive $ opts+        ShiftPositiveFlag . shiftPositive $ opts       customLabel' = (\ xs -> bool                                 (fmap (Just . CustomLabel) xs)                                 (repeat Nothing)                             . null                             $ xs                       )-                    . unHelpful                     . customLabel                     $ opts       readOrErr err = fromMaybe (error err) . readMaybe       filterThresholds'  = fmap FilterThresholds                          . fmap (readOrErr "Cannot read --filter-thresholds")-                         . unHelpful                          . filterThresholds                          $ opts       customRegions' = CustomRegions                      . fmap (either (\x -> error $ "Cannot parse region format `chrN:START-END` in: " <> x) id . parseChrRegion)-                     . unHelpful                      . customRegion                      $ opts @@ -202,7 +207,7 @@                   . fmap wait                   . mapReduce workers                   . zipWith (\l -> fmap (labelRows l)) customLabel'  -- Depending on which axis to label from transpose.-                  . fmap (fmap (fastBinJoinCols binWidth' True) . loadSSM opts)  -- Load matrices, possible preparing for fast bin joining+                  . fmap (fmap (fastBinJoinCols binWidth' True) . loadSSM sub opts)  -- Load matrices, possible preparing for fast bin joining                   $ matrixPaths'    let sc = maybe@@ -217,6 +222,8 @@       normMat TotalNorm      = totalScaleSparseMat       normMat (LogCPMNorm b) = logCPMSparseMat b       normMat QuantileNorm   = quantileScaleSparseMat+      normMat MinMaxNorm     = minMaxNormSparseMat+      normMat TransposeNorm  = transposeSparseMat       normMat NoneNorm       = id       applyNorms mat = foldl' (\acc norm -> L.over matrix (normMat norm) acc) mat normalizations'       processSc = L.over matrix (MatObsRow . S.sparsifySM . unMatObsRow)

@@ -12,6 +12,7 @@ {-# LANGUAGE PackageImports    #-} {-# LANGUAGE TypeOperators     #-} {-# LANGUAGE TupleSections     #-}+{-# LANGUAGE DuplicateRecordFields #-}  module TooManyCells.Program.MakeTree where @@ -35,7 +36,6 @@ import Math.Clustering.Hierarchical.Spectral.Types (getClusterItemsDend, EigenGroup (..)) import Math.Clustering.Spectral.Sparse (b1ToB2, B1 (..), B2 (..)) import Math.Modularity.Types (Q (..))-import Options.Generic import System.IO (hPutStrLn, stderr) import Text.Read (readMaybe, readEither) import TextShow (showt)@@ -74,101 +74,77 @@ import TooManyCells.Program.Options import TooManyCells.Program.Utility -makeTreeMain :: Options -> IO ()-makeTreeMain opts = H.withEmbeddedR defaultConfig $ do+makeTreeMain :: Subcommand -> IO ()+makeTreeMain sub@(MakeTreeCommand opts) = H.withEmbeddedR defaultConfig $ do     let readOrErr err = fromMaybe (error err) . readMaybe-        matrixPaths'      = unHelpful . matrixPath $ opts+        matrixPaths'      = matrixPath . (loadMatrixOptions :: MakeTree -> LoadMatrixOptions) $ opts         labelsFile'       =-            fmap LabelFile . unHelpful . labelsFile $ opts+            fmap LabelFile . (labelsFile :: MakeTree -> Maybe String) $ opts         prior'            =-            fmap PriorPath . unHelpful . prior $ opts+            fmap PriorPath . (prior :: MakeTree -> Maybe String) $ opts         updateTreeRows'   =-          UpdateTreeRowsFlag . not . unHelpful . noUpdateTreeRows $ opts-        delimiter'        =-            Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts-        eigenGroup'       =-            maybe SignGroup (readOrErr "Cannot read --eigen-group.")-              . unHelpful-              . eigenGroup-              $ opts-        dense'            = DenseFlag . unHelpful . dense $ opts-        normalizations'   = getNormalization opts-        numEigen'         = fmap NumEigen . unHelpful . numEigen $ opts-        numRuns'          = fmap NumRuns . unHelpful . numRuns $ opts-        minSize'          = fmap MinClusterSize . unHelpful . minSize $ opts-        maxStep'          = fmap MaxStep . unHelpful . maxStep $ opts+          UpdateTreeRowsFlag . not . (noUpdateTreeRows :: MakeTree -> Bool) $ opts+        delimiter'        = Delimiter+                          . (delimiter :: LoadMatrixOptions -> Char)+                          . (loadMatrixOptions :: MakeTree -> LoadMatrixOptions)+                          $ opts+        eigenGroup'       = eigenGroup opts+        dense'            = DenseFlag . dense $ opts+        normalizations'   = getNormalization sub+                          . normalization+                          $ (loadMatrixOptions :: MakeTree -> LoadMatrixOptions) opts+        numEigen'         = fmap NumEigen . numEigen $ opts+        numRuns'          = fmap NumRuns . numRuns $ opts+        minSize'          = fmap MinClusterSize . minSize $ opts+        maxStep'          = fmap MaxStep . maxStep $ opts         maxProportion'    =-            fmap MaxProportion . unHelpful . maxProportion $ opts-        minDistance'       = fmap MinDistance . unHelpful . minDistance $ opts-        minModularity'     = fmap Q . unHelpful . minModularity $ opts-        minDistanceSearch' = fmap MinDistanceSearch . unHelpful . minDistanceSearch $ opts-        smartCutoff'      = fmap SmartCutoff . unHelpful . smartCutoff $ opts+            fmap MaxProportion . maxProportion $ opts+        minDistance'       = fmap MinDistance . minDistance $ opts+        minModularity'     = fmap Q . minModularity $ opts+        minDistanceSearch' = fmap MinDistanceSearch . minDistanceSearch $ opts+        smartCutoff'      = fmap SmartCutoff . smartCutoff $ opts         elbowCutoff'      =           fmap ( ElbowCutoff                . readOrErr "Cannot read --elbow-cutoff."                )-            . unHelpful             . elbowCutoff             $ opts-        customCut'        = CustomCut . Set.fromList . unHelpful . customCut $ opts-        rootCut'          = fmap RootCut . unHelpful . rootCut $ opts+        customCut'        = CustomCut . Set.fromList . customCut $ opts+        rootCut'          = fmap RootCut . rootCut $ opts         dendrogramOutput' = DendrogramFile                           . fromMaybe "dendrogram.svg"-                          . unHelpful                           . dendrogramOutput                           $ opts         matrixOutput'     = fmap (getMatrixOutputType . (unOutputDirectory output' FP.</>))-                          . unHelpful                           . matrixOutput                           $ opts         labelMapOutputFlag' =-            LabelMapOutputFlag . unHelpful . labelsOutput $ opts+            LabelMapOutputFlag . labelsOutput $ opts         fragmentsOutputFlag' =-            FragmentsOutputFlag . unHelpful . fragmentsOutput $ opts+            FragmentsOutputFlag . fragmentsOutput $ opts         drawLeaf'         =             maybe               (maybe DrawText (const (DrawItem DrawLabel)) labelsFile')               (readOrErr "Cannot read --draw-leaf. If using DrawContinuous, remember to put features in a list: DrawItem (DrawContinuous [\\\"FEATURE\\\"])")-                . unHelpful                 . drawLeaf                 $ opts-        drawCollection'   =-            maybe PieChart (readOrErr "Cannot read --draw-collection.")-              . unHelpful-              . drawCollection-              $ opts-        drawMark'         = maybe MarkNone (readOrErr "Cannot read --draw-mark.")-                          . unHelpful-                          . drawMark-                          $ opts-        drawNodeNumber'   = DrawNodeNumber . unHelpful . drawNodeNumber $ opts-        drawMaxNodeSize'  =-            DrawMaxNodeSize . fromMaybe 72 . unHelpful . drawMaxNodeSize $ opts+        drawCollection'   = drawCollection opts+        drawMark'         = drawMark opts+        drawNodeNumber'   = DrawNodeNumber . drawNodeNumber $ opts+        drawMaxNodeSize'  = DrawMaxNodeSize . drawMaxNodeSize $ opts         drawMaxLeafNodeSize' = DrawMaxLeafNodeSize                              . fromMaybe (unDrawMaxNodeSize drawMaxNodeSize')-                             . unHelpful                              . drawMaxLeafNodeSize                              $ opts         drawNoScaleNodes' =-            DrawNoScaleNodesFlag . unHelpful . drawNoScaleNodes $ opts-        drawLegendSep'    = DrawLegendSep-                          . fromMaybe 1-                          . unHelpful-                          . drawLegendSep-                          $ opts-        drawLegendAllLabels' =-            DrawLegendAllLabels . unHelpful . drawLegendAllLabels $ opts-        drawPalette' = maybe-                        Set1-                        (fromMaybe (error "Cannot read --palette.") . readMaybe)-                     . unHelpful-                     . drawPalette-                     $ opts+            DrawNoScaleNodesFlag . drawNoScaleNodes $ opts+        drawLegendSep'    = DrawLegendSep . drawLegendSep $ opts+        drawLegendAllLabels' = DrawLegendAllLabels . drawLegendAllLabels $ opts+        drawPalette' = drawPalette opts         drawColors'       = fmap ( CustomColors                                  . fmap sRGB24read                                  . (\x -> readOrErr "Cannot read --draw-colors." x :: [String])                                  )-                          . unHelpful                           . drawColors                           $ opts         drawDiscretize' = (=<<) (\x -> either error Just@@ -181,36 +157,31 @@                                     (Right . CustomColorMap . fmap sRGB24read)                                 $ (readEither x :: Either String [String])                                 )-                        . unHelpful                         . drawDiscretize                         $ opts           where             finalError err x = "Error in draw-discretize: " <> err <> " " <> x         drawScaleSaturation' =-            fmap DrawScaleSaturation . unHelpful . drawScaleSaturation $ opts+            fmap DrawScaleSaturation . drawScaleSaturation $ opts         drawItemLineWeight' = fmap DrawItemLineWeight-                            . unHelpful                             . drawItemLineWeight                             $ opts-        drawFont' = fmap DrawFont . unHelpful . drawFont $ opts-        drawBarBounds' = DrawBarBounds . unHelpful . drawBarBounds $ opts-        order'            = Order . fromMaybe 1 . unHelpful . order $ opts-        clumpinessMethod' =-            maybe Clump.Majority (readOrErr "Cannot read --clumpiness-method.")-              . unHelpful-              . clumpinessMethod-              $ opts-        projectionFile' =-            fmap ProjectionFile . unHelpful . projectionFile $ opts-        output'           =-            OutputDirectory . fromMaybe "out" . unHelpful . output $ opts+        drawFont' = fmap DrawFont . drawFont $ opts+        drawBarBounds' = DrawBarBounds . drawBarBounds $ opts+        order'            = Order . (order :: MakeTree -> Double) $ opts+        clumpinessMethod' = clumpinessMethod opts+        projectionFile' = fmap ProjectionFile+                        . (projectionFile :: MakeTree -> Maybe String)+                        $ opts+        output' = OutputDirectory . (output :: MakeTree -> String) $ opts      pb <- Progress.newProgressBar Progress.defStyle 10 (Progress.Progress 0 11 ())     -- Increment  progress bar.     Progress.incProgress pb 1      -- Load matrix once.-    scRes <- loadAllSSM opts+    scRes <- loadAllSSM sub+           $ (loadMatrixOptions :: MakeTree -> LoadMatrixOptions) opts     let processedSc = fmap fst scRes         customLabelMap = join . fmap snd $ scRes @@ -220,19 +191,6 @@     -- Where to place output files.     FP.createDirectoryIfMissing True . unOutputDirectory $ output' -    -- Get the label map from either a file or from expression thresholds.-    labelMap <- case drawLeaf' of-                    (DrawItem (DrawThresholdContinuous gs)) ->-                        return-                            . Just-                            . getLabelMapThresholdContinuous-                                (fmap (L.over L._1 Feature) gs)-                            . extractSc-                            $ processedSc-                    _ -> if isJust labelsFile'-                          then mapM (loadLabelData delimiter') labelsFile'-                          else return customLabelMap-     -- Increment  progress bar.     Progress.incProgress pb 1 @@ -262,12 +220,40 @@             fmap (L.over clusterDend (updateTreeRowBool updateTreeRows' processedSc))               $ loadClusterResultsFiles clInput treeInput +    projectionMap <- mapM loadProjectionMap projectionFile'++    -- Only used for proximity labels for now.+    let oldGr = treeToGraph+              . updateTreeRowBool updateTreeRows' processedSc+              . _clusterDend+              $ originalClusterResults++    -- Get the label map from either a file or from expression thresholds.+    labelMap <- case drawLeaf' of+                    (DrawItem (DrawThresholdContinuous gs)) ->+                        return+                            . Just+                            . getLabelMapThresholdContinuous+                                (fmap (L.over L._1 Feature) gs)+                            . extractSc+                            $ processedSc+                    (DrawItem (DrawProximity ps)) -> return $ do+                      coordMap <- fmap projectionMapToCoordinateMap projectionMap+                      return $ getLabelMapProximity+                                oldGr+                                coordMap+                                ps+                    _ -> if isJust labelsFile'+                          then mapM (loadLabelData delimiter') labelsFile'+                          else return customLabelMap++     -- Increment  progress bar.     Progress.incProgress pb 1      let birchMat = processedSc         birchSimMat =-            case (not . null . unHelpful . matrixPath $ opts, drawCollection') of+            case (not . null . matrixPath . (loadMatrixOptions :: MakeTree -> LoadMatrixOptions) $ opts, drawCollection') of                 (True, CollectionGraph{} )  ->                     Just                         . B2Matrix@@ -472,15 +458,12 @@         H.io $ Progress.incProgress pb 1          -- Plot clustering of the projections.-        case projectionFile' of+        case projectionMap of           Nothing  -> return ()-          (Just f) -> do-            -- Load projection map if provided.-            projectionMap <- H.io $ loadProjectionMap f-+          (Just pm) -> do             plotClustersR               (unOutputDirectory output' FP.</> "projection.pdf")-              projectionMap+              pm               . _clusterList               $ clusterResults @@ -489,7 +472,7 @@                 (Just lm, Just icm) ->                     plotLabelClustersR                         (unOutputDirectory output' FP.</> "label_projection.pdf")-                        projectionMap+                        pm                         lm                         icm                         (_clusterList clusterResults)@@ -499,3 +482,4 @@         H.io $ Progress.incProgress pb 1          return ()+makeTreeMain _ = error "Wrong path in make-tree, contact Gregory Schwartz for this error."

@@ -5,11 +5,11 @@ -}  {-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE DuplicateRecordFields #-}  module TooManyCells.Program.MatrixOutput where  -- Remote-import Options.Generic import qualified System.Directory as FP  -- Local@@ -19,13 +19,14 @@ import TooManyCells.Matrix.Utility import TooManyCells.Matrix.Types -matrixOutputMain :: Options -> IO ()-matrixOutputMain opts = do-    let matOutput' = getMatrixOutputType . unHelpful . matOutput $ opts+matrixOutputMain :: Subcommand -> IO ()+matrixOutputMain sub@(MatrixOutputCommand opts) = do+    let matOutput' = getMatrixOutputType . matOutput $ opts      -- Load matrix once.-    scRes <- loadAllSSM opts+    scRes <- loadAllSSM sub $ (loadMatrixOptions :: MatrixOutput -> LoadMatrixOptions) opts     let processedSc = fmap fst scRes      -- Write matrix     writeMatrixLike (MatrixTranspose False) matOutput' . extractSc $ processedSc+matrixOutputMain _ = error "Wrong path in matrix-output, contact Gregory Schwartz for this error."

@@ -5,13 +5,13 @@ -}  {-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE DuplicateRecordFields #-}  module TooManyCells.Program.Motifs where  -- Remote import Control.Monad (mfilter) import Data.Maybe (fromMaybe)-import Options.Generic import TextShow (showt) import qualified Data.Text as T import qualified Turtle as TU@@ -23,28 +23,20 @@ import TooManyCells.Program.Options  -- | Motif path.-motifsMain :: Options -> IO ()-motifsMain opts = do-  let genome' = GenomeFile . unHelpful . motifGenome $ opts-      topN' = TopN . fromMaybe 100 . unHelpful . topN $ opts-      motifCommand' = MotifCommand-                    . fromMaybe "meme %s -nmotifs 50 -oc %s"-                    . unHelpful-                    . motifCommand-                    $ opts+motifsMain :: Subcommand -> IO ()+motifsMain (MotifsCommand opts) = do+  let genome' = GenomeFile . motifGenome $ opts+      topN' = TopN . (topN :: Motifs -> Int) $ opts+      motifCommand' = MotifCommand . motifCommand $ opts       motifGenomeCommand' = fmap MotifGenomeCommand-                          . unHelpful                           . motifGenomeCommand                           $ opts-      diffFile' = DiffFile . TU.fromText . unHelpful . diffFile $ opts+      diffFile' = DiffFile . TU.fromText . diffFile $ opts       backgroundDiffFile' = fmap (BackgroundDiffFile . TU.fromText)-                          . unHelpful                           . backgroundDiffFile                           $ opts       outDir = OutputDirectory-             . fromMaybe "out"-             . unHelpful-             . output+             . (output :: Motifs -> String)              $ opts    TU.mktree . TU.fromText . T.pack . unOutputDirectory $ outDir

@@ -4,303 +4,434 @@ Options for the command line program. -} -{-# LANGUAGE DataKinds         #-}-{-# LANGUAGE DeriveGeneric     #-}-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE TypeOperators     #-}+-- {-# LANGUAGE DataKinds         #-}+-- {-# LANGUAGE DeriveGeneric     #-}+-- {-# LANGUAGE OverloadedStrings #-}+-- {-# LANGUAGE TypeOperators     #-}+{-# LANGUAGE ApplicativeDo #-}+{-# LANGUAGE DuplicateRecordFields #-}+{-# LANGUAGE PackageImports    #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE StandaloneDeriving #-}  module TooManyCells.Program.Options where  -- Remote-import Options.Generic+import BirchBeer.Types (DrawNodeMark (MarkNone), DrawCollection (PieChart), DrawLeaf (DrawText), Palette (..))+import TooManyCells.Paths.Types (PathDistance (PathStep))+import Math.Clustering.Hierarchical.Spectral.Types (EigenGroup (SignGroup))+import Options.Applicative import qualified Data.Text as T+import qualified "find-clumpiness" Types as Clump  -- Local+import TooManyCells.Matrix.Types (NormType)+import TooManyCells.Spatial.Types (TopDistances) +deriving instance Show Clump.Exclusivity+ -- | Command line arguments-data Options-    = MakeTree { matrixPath :: [String] <?> "(PATH) The path to the input directory containing the matrix output of cellranger (cellranger < 3 (matrix.mtx, genes.tsv, and barcodes.tsv) or cellranger >= 3 (matrix.mtx.gz, features.tsv.gz, and barcodes.tsv.gz) or an input csv file containing feature row names and cell column names. scATAC-seq is supported if input file contains \"fragments\", ends with \".tsv.gz\" (such as \"fragments.tsv.gz\" or \"sample1_fragments.tsv.gz\"), and is in the SORTED (sort -k1,1 -k2,2n) 10x fragments format (see also --binwidth, --no-binarize). If given as a list (--matrix-path input1 --matrix-path input2 etc.) then will join all matrices together. Assumes the same number and order of features in each matrix, so only cells are added."-               , binwidth :: Maybe Int <?> "(Nothing | BINSIZE) If input data has region features in the format of `chrN:START-END`, BINSIZE input is required to convert ranges to fixed width bins."-               , noBinarize :: Bool <?> "If a fragments.tsv.gz file, do not binarize data."-               , binarize :: Bool <?> "Binarize data. Default for fragments.tsv.gz."-               , excludeMatchFragments :: Maybe String <?> "([Nothing] | STRING) Exclude fragments from fragments.tsv.gz file if STRING is infix of fragment row. For instance, exclude all chrY fragments with \"--exclude-match-fragments chrY\"."-               , blacklistRegionsFile :: Maybe T.Text <?> "([Nothing] | FILE) Bed file containing regions to ignore. Any fragments overlapping these regions are ignored if the input is not a fragments file."-               , customRegion :: [T.Text] <?> "([Nothing] | chrN:START-END) Only look at these regions in the matrix for chromosome region features. A list in the format `chrN:START-END`, will assign values for a cell intersecting this region to this region. If a cell region overlaps two or more regions, will be counting that many times in the new features. Example input: `--custom-regions \"chr1:1003-10064\" --custom-regions \"chr2:3021-5034\"` etc."-               , projectionFile :: Maybe String <?> "([Nothing] | FILE) The input file containing positions of each cell for plotting. Format with header is \"barcode,x,y\". Useful for 10x where a TNSE projection is generated in \"projection.csv\". Cells without projections will not be plotted. If not supplied, no plot will be made."-               , cellWhitelistFile :: Maybe String <?> "([Nothing] | FILE) The input file containing the cells to include. No header, line separated list of barcodes."-               , labelsFile :: Maybe String <?> "([Nothing] | FILE) The input file containing the label for each cell barcode, with \"item,label\" header."-               , customLabel :: [T.Text] <?> "([] | [LABEL]) List of labels to assign each matrix if all cells from each matrix are given the same label per matrix. This argument intends to simplify the process of labeling by bypassing --labels-file if the user just wants each matrix to have its own label (i.e. sample). Must be the same length and order as --matrix-path: for instance, --matrix-path input1 --custom-label sample1 --matrix-path input2 --custom-label sample2 etc. will label all cells from input1 with sample1, input2 with sample2, etc. If there are multiple labels per matrix, you must use --labels-file."-               , delimiter :: Maybe Char <?> "([,] | CHAR) The delimiter for the most csv files in the program. For instance, if using a normal csv rather than cellranger output and for --labels-file."-               , featureColumn :: Maybe Int <?> "([1] | COLUMN) The column (1-indexed) in the features.tsv.gz file to use for feature names. If using matrix market format, cellranger stores multiple columns in the features file, usually the first column for the Ensembl identifier and the second column for the gene symbol. If the Ensembl identifier is not quickly accessible, use --feature-column 2 for the second column, which is usually more ubiquitous. Useful for overlaying gene expression so you can say --draw-leaf \"DrawItem (DrawContinuous \\\"CD4\\\")\") instead of --draw-leaf \"DrawItem (DrawContinuous \\\"ENSG00000010610\\\")\"). Does not affect CSV format (the column names will be the feature names)."-               , normalization :: [String] <?> "([TfIdfNorm] | UQNorm | MedNorm | TotalMedNorm | TotalNorm | LogCPMNorm DOUBLE | QuantileNorm | NoneNorm) Type of normalization before clustering. Can be used as a list to perform one after the other: --normalization QuantileNorm --normalization TfIdfNorm will first apply quantile then tf-idf normalization. TfIdfNorm normalizes based on the prevalence of each feature. UQNorm normalizes each observation by the upper quartile non-zero counts of that observation. MedNorm normalizes each observation by the median non-zero counts of that observation. TotalNorm normalizes each cell by the total count. TotalMedNorm normalized first each observation by total count then by median of non-zero counts across features. (LogCPMNorm DOUBLE) normalizes by logB(CPM + 1) where B is DOUBLE. QuantileNorm normalizes by quantile normalization, ignores zeros, may be slow. NoneNorm does not normalize. Default is TfIdfNorm for clustering and NoneNorm for differential (which instead uses the recommended edgeR single cell preprocessing including normalization and filtering, any normalization provided here will result in edgeR preprocessing on top). Older versions had BothNorm which has been replaced with --normalization TotalMedNorm --normalization TfIdfNorm. This change will also affect other analyses, as TfIdfNorm will now be in non-cluster-related entry points if specified, so --normalization UQNorm from < v2.0.0.0 is now --normalization UQNorm --normalization TfIdfNorm."-               , eigenGroup :: Maybe String <?> "([SignGroup] | KMeansGroup) Whether to group the eigenvector using the sign or kmeans while clustering. While the default is sign, kmeans may be more accurate (but starting points are arbitrary)."-               , numEigen :: Maybe Int <?> "([1] | INT) Number of eigenvectors to use while clustering with kmeans. Takes from the second to last eigenvector. Recommended to start at 1 and work up from there if needed. May help offset the possible instability and inaccuracy of SVDLIBC."-               , numRuns :: Maybe Int <?> "([Nothing] | INT) Number of runs for permutation test at each split for modularity. Defaults to no test."-               , minSize :: Maybe Int <?> "([1] | INT) The minimum size of a cluster. Defaults to 1."-               , maxStep :: Maybe Int <?> "([Nothing] | INT) Only keep clusters that are INT steps from the root. Defaults to all steps."-               , maxProportion :: Maybe Double <?> "([Nothing] | DOUBLE) Stopping criteria to stop at the node immediate after a node with DOUBLE proportion split. So a node N with L and R children will stop with this criteria at 0.5 if |L| / |R| < 0.5 or > 2 (absolute log2 transformed), that is, if one child has over twice as many items as the other child. Includes L and R in the final result."-               , minModularity :: Maybe Double <?> "([Nothing] | DOUBLE) Nearly the same as --min-distance, but for clustering instead of drawing (so the output json tree can be larger). Stopping criteria to stop at the node with DOUBLE modularity. So a node N with L and R children will stop with this criteria the distance at N to L and R is < DOUBLE. Does not include L and R in the final result."-               , minDistance :: Maybe Double <?> "([Nothing] | DOUBLE) Stopping criteria to stop at the node immediate after a node with DOUBLE distance. So a node N with L and R children will stop with this criteria the distance at N to L and R is < DOUBLE. Includes L and R in the final result."-               , minDistanceSearch :: Maybe Double <?> "([Nothing] | DOUBLE) Similar to --min-distance, but searches from the leaves to the root -- if a path from a subtree contains a distance of at least DOUBLE, keep that path, otherwise prune it. This argument assists in finding distant nodes."-               , smartCutoff :: Maybe Double <?> "([Nothing] | DOUBLE) Whether to set the cutoffs for --min-size, --max-proportion, --min-distance, and --min-distance-search based off of the distributions (median + (DOUBLE * MAD)) of all nodes. To use smart cutoffs, use this argument and then set one of the four arguments to an arbitrary number, whichever cutoff type you want to use. --max-proportion and --min-size distributions are log2 transformed."-               , elbowCutoff :: Maybe String <?> "(Max | Min) Whether to set the cutoffs for --min-size, --max-proportion, --min-distance, and --min-distance-search based off of the elbow point of distributions of all nodes. For a distribution in positive x and y on a graph, the top left hump would be Max and the bottom right dip would be Min. To use elbow cutoffs, use this argument and then set one of the three arguments to an arbitrary number, whichever cutoff type you want to use. --max-proportion and --min-size distributions are log2 transformed. Conflicts with --smart-cutoff, so this argument takes precedent."-               , customCut :: [Int] <?> "([Nothing] | NODE) List of nodes to prune (make these nodes leaves). Invoked by --custom-cut 34 --custom-cut 65 etc."-               , rootCut :: Maybe Int <?> "([Nothing] | NODE) Assign a new root to the tree, removing all nodes outside of the subtree."-               , dendrogramOutput :: Maybe String <?> "([dendrogram.svg] | FILE) The filename for the dendrogram. Supported formats are PNG, PS, PDF, and SVG."-               , matrixOutput :: Maybe String <?> "([Nothing] | FOLDER | FILE.csv) Output the filtered and normalized (not including TfIdfNorm) matrix in this folder under the --output directory in matrix market format or, if a csv file is specified, a dense csv format. Like input, features are rows."-               , labelsOutput :: Bool <?> "Whether to write the labels used for each observation as a labels.csv file in the output folder."-               , fragmentsOutput :: Bool <?> "Whether to output fragments_tsv.gz with barcodes altered by --custom-label in the output folder (excludes filtered-out cells). Useful for downstream analysis by the peaks entry point where the cluster barcodes differ from the original fragments.tsv.gz file when using --custom-label. Matches barcodes based on BARCODE-LABEL."-               , matrixTranspose :: Bool <?> "Whether to transpose the matrix before all processing (observations become features and vice-versa). Will be affected by other options (check your filtering thresholds, normalizations, etc!)"-               , drawLeaf :: Maybe String <?> "([DrawText] | DrawItem DrawItemType) How to draw leaves in the dendrogram. DrawText is the number of items in that leaf. DrawItem is the collection of items represented by circles, consisting of: DrawItem DrawLabel, where each item is colored by its label, DrawItem (DrawContinuous [FEATURE]), where each item is colored by the expression of FEATURE (corresponding to a feature name in the input matrix, [FEATURE] is a list, so if more than one FEATURE is listed, uses the average of the feature values), DrawItem (DrawThresholdContinuous [(FEATURE, THRESHOLD)]), where each item is colored by the binary high / low expression of FEATURE based on THRESHOLD (either `Exact DOUBLE` or `MadMedian DOUBLE`, where Exact just uses the DOUBLE as a cutoff value while MadMedian uses the DOUBLE as the number of MADs away from the median value of the feature) and multiple FEATUREs can be used to combinatorically label items (FEATURE1 high / FEATURE2 low, etc.), DrawItem DrawSumContinuous, where each item is colored by the sum of the post-normalized columns (use --normalization NoneNorm for UMI counts, default), and DrawItem DrawDiversity, where each node is colored by the diversity based on the labels of each item and the color is normalized separately for the leaves and the inner nodes. The default is DrawText, unless --labels-file is provided, in which DrawItem DrawLabel is the default. If the label or feature cannot be found, the default color will be black (check your spelling!)."-               , drawCollection :: Maybe String <?> "([PieChart] | PieRing | IndividualItems | Histogram | NoLeaf | CollectionGraph MAXWEIGHT THRESHOLD [NODE]) How to draw item leaves in the dendrogram. PieRing draws a pie chart ring around the items. PieChart only draws a pie chart instead of items. IndividualItems only draws items, no pie rings or charts. Histogram plots a histogram of the features requested. NoLeaf has no leaf, useful if there are so many items the tree takes very long to render. (CollectionGraph MAXWEIGHT THRESHOLD [NODE]) draws the nodes and edges within leaves that are descendents of NODE (empty list [] indicates draw all leaf networks) based on the input matrix, normalizes edges based on the MAXWEIGHT, and removes edges for display less than THRESHOLD (after normalization, so for CollectionGraph 2 0.5 [26], draw the leaf graphs for all leaves under node 26, then a edge of 0.7 would be removed because (0.7 / 2) < 0.5). For CollectionGraph with no colors, use --draw-leaf \"DrawItem DrawLabel\" and all nodes will be black. If you don't specify this option, DrawText from --draw-leaf overrides this argument and only the number of cells will be plotted."-               , drawMark :: Maybe String <?> "([MarkNone] | MarkModularity | MarkSignificance ) How to draw annotations around each inner node in the tree. MarkNone draws nothing and MarkModularity draws a black circle representing the modularity at that node, darker black means higher modularity for that next split. MarkSignificance is for significance, i.e. p-value, darker means higher value."-               , drawNodeNumber :: Bool <?> "Draw the node numbers on top of each node in the graph."-               , drawMaxNodeSize :: Maybe Double <?> "([72] | DOUBLE) The max node size when drawing the graph. 36 is the theoretical default, but here 72 makes for thicker branches."-               , drawMaxLeafNodeSize :: Maybe Double <?> "([--draw-max-node-size] | DOUBLE) The max leaf node size when drawing the graph. Defaults to the value of --draw-max-node-size."-               , drawNoScaleNodes :: Bool <?> "Do not scale inner node size when drawing the graph. Instead, uses draw-max-node-size as the size of each node and is highly recommended to change as the default may be too large for this option."-               , drawLegendSep :: Maybe Double <?> "([1] | DOUBLE) The amount of space between the legend and the tree."-               , drawLegendAllLabels :: Bool <?> "Whether to show all the labels in the label file instead of only showing labels within the current tree. The program generates colors from all labels in the label file first in order to keep consistent colors. By default, this value is false, meaning that only the labels present in the tree are shown (even though the colors are the same). The subset process occurs after --draw-colors, so when using that argument make sure to account for all labels."-               , drawPalette :: Maybe String <?> "([Set1] | Hsv | Ryb | Blues) Palette to use for legend colors. With high saturation in --draw-scale-saturation, consider using Hsv to better differentiate colors."-               , drawColors :: Maybe String <?> "([Nothing] | COLORS) Custom colors for the labels or continuous features. Will repeat if more labels than provided colors. For continuous feature plots, uses first two colors [high, low], defaults to [red, gray]. For instance: --draw-colors \"[\\\"#e41a1c\\\", \\\"#377eb8\\\"]\""-               , drawDiscretize :: Maybe String <?> "([Nothing] | COLORS | INT) Discretize colors by finding the nearest color for each item and node. For instance, --draw-discretize \"[\\\"#e41a1c\\\", \\\"#377eb8\\\"]\" will change all node and item colors to one of those two colors, based on Euclidean distance. If using \"--draw-discretize INT\", will instead take the default map and segment (or interpolate) it into INT colors, rather than a more continuous color scheme. May have unintended results when used with --draw-scale-saturation."-               , drawScaleSaturation :: Maybe Double <?> "([Nothing] | DOUBLE) Multiply the saturation value all nodes by this number in the HSV model. Useful for seeing more visibly the continuous colors by making the colors deeper against a gray scale."-               , drawItemLineWeight :: Maybe Double <?> "([0.1] | DOUBLE) The line weight for items in the leaves if shown. Supplied as if there are too many items, the collection may look like a black box. Set to 0 to disable outlines of items to avoid this."-               , drawFont :: Maybe String <?> "([Arial] | FONT) Specify the font to use for the labels when plotting."-               , drawBarBounds :: Bool <?> "Whether to plot only the minimum and maximum ticks for the color bars."-               , pca :: Maybe Int <?> "([Nothing] | INT) Not recommended, as it makes cosine similarity less meaningful (therefore less accurate -- instead, consider making your own similarity matrix and using cluster-tree, our sister algorithm, to cluster the matrix and plot with birch-beer). The number of dimensions to keep for PCA dimensionality reduction before clustering. Default is no PCA at all in order to keep all information. Should use with --shift-positive to ensure no negative values (as --pca will center and scale). Consider changing the modularity cutoff to a lower value (such as --min-modularity -0.5)."-               , lsa :: Maybe Int <?> "([Nothing] | INT) The number of dimensions to keep for LSA dimensionality reduction. Uses TD-IDF followed by SVD before clustering, same warnings as --pca apply, including the use of --shift-positive with possible --min-modularity -0.5."-               , svd :: Maybe Int <?> "([Nothing] | INT) The number of dimensions to keep for SVD dimensionality reduction. Will center and scale, same warnings as --pca apply, including the use of --shift-positive with possible --min-modularity -0.5."-               , dropDimension :: Bool <?> "Instead of keeping dimensions for --pca, --lsa, or --svd, drop the number of dimensions (i.e. --lsa 1 drops the first projection which may have a strong batch effect)."-               , shiftPositive :: Bool <?> "Shift features to positive values. Positive values are shifted to allow modularity to work correctly."-               , filterThresholds :: Maybe String <?> "([Nothing] | (DOUBLE, DOUBLE)) The minimum filter thresholds for (MINCELL, MINFEATURE) when filtering cells and features by low read counts. Default changed to Nothing due to additional supported assays. To use the original filter thresholds, use --filter-thresholds \"(250, 1)\"."-               , prior :: Maybe String <?> "([Nothing] | STRING) The input folder containing the output from a previous run. If specified, skips clustering by using the previous clustering files."-               , noUpdateTreeRows :: Bool <?> "Don't update the row indices of a tree if using an identical matrix to the one which generated the tree. Should not be used unless the matrix to make the tree is identical, then can result in speedup."-               , order :: Maybe Double <?> "([1] | DOUBLE) The order of diversity."-               , clumpinessMethod :: Maybe String <?> "([Majority] | Exclusive | AllExclusive) The method used when calculating clumpiness: Majority labels leaves according to the most abundant label, Exclusive only looks at leaves consisting of cells solely from one label, and AllExclusive treats the leaf as containing both labels."-               , dense :: Bool <?> "Whether to use dense matrix algorithms for clustering. Should be faster for dense matrices, so if batch correction, PCA, or other algorithms are applied upstream to the input matrix, consider using this option to speed up the tree generation."-               , output :: Maybe String <?> "([out] | STRING) The folder containing output."}-    | Interactive { matrixPath :: [String] <?> "(PATH) The path to the input directory containing the matrix output of cellranger (cellranger < 3 (matrix.mtx, genes.tsv, and barcodes.tsv) or cellranger >= 3 (matrix.mtx.gz, features.tsv.gz, and barcodes.tsv.gz) or an input csv file containing feature row names and cell column names. scATAC-seq is supported if input file contains \"fragments\", ends with \".tsv.gz\" (such as \"fragments.tsv.gz\" or \"sample1_fragments.tsv.gz\"), and is in the SORTED (sort -k1,1 -k2,2n) 10x fragments format (see also --binwidth, --no-binarize). If given as a list (--matrix-path input1 --matrix-path input2 etc.) then will join all matrices together. Assumes the same number and order of features in each matrix, so only cells are added."-                  , matrixTranspose :: Bool <?> "Whether to transpose the matrix before all processing (observations become features and vice-versa). Will be affected by other options (check your filtering thresholds, normalizations, etc!)"-                  , binwidth :: Maybe Int <?> "(Nothing | BINSIZE) If input data has region features in the format of `chrN:START-END`, BINSIZE input is required to convert ranges to fixed width bins."-                  , noBinarize :: Bool <?> "If a fragments.tsv.gz file, do not binarize data."-                  , binarize :: Bool <?> "Binarize data. Default for fragments.tsv.gz."-                  , excludeMatchFragments :: Maybe String <?> "([Nothing] | STRING) Exclude fragments from fragments.tsv.gz file if STRING is infix of fragment row. For instance, exclude all chrY fragments with \"--exclude-match-fragments chrY\"."-                  , blacklistRegionsFile :: Maybe T.Text <?> "([Nothing] | FILE) Bed file containing regions to ignore. Any fragments overlapping these regions are ignored if the input is not a fragments file."-                  , customRegion :: [T.Text] <?> "([Nothing] | chrN:START-END) Only look at these regions in the matrix for chromosome region features. A list in the format `chrN:START-END`, will assign values for a cell intersecting this region to this region. If a cell region overlaps two or more regions, will be counting that many times in the new features. Example input: `--custom-regions \"chr1:1003-10064\" --custom-regions \"chr2:3021-5034\"` etc."-                  , cellWhitelistFile :: Maybe String <?> "([Nothing] | FILE) The input file containing the cells to include. No header, line separated list of barcodes."-                  , labelsFile :: Maybe String <?> "([Nothing] | FILE) The input file containing the label for each cell barcode, with \"item,label\" header."-                  , customLabel :: [T.Text] <?> "([] | [LABEL]) List of labels to assign each matrix if all cells from each matrix are given the same label per matrix. This argument intends to simplify the process of labeling by bypassing --labels-file if the user just wants each matrix to have its own label (i.e. sample). Must be the same length and order as --matrix-path: for instance, --matrix-path input1 --custom-label sample1 --matrix-path input2 --custom-label sample2 etc. will label all cells from input1 with sample1, input2 with sample2, etc. If there are multiple labels per matrix, you must use --labels-file."-                  , delimiter :: Maybe Char <?> "([,] | CHAR) The delimiter for the most csv files in the program. For instance, if using a normal csv rather than cellranger output and for --labels-file."-                  , featureColumn :: Maybe Int <?> "([1] | COLUMN) The column (1-indexed) in the features.tsv.gz file to use for feature names. If using matrix market format, cellranger stores multiple columns in the features file, usually the first column for the Ensembl identifier and the second column for the gene symbol. If the Ensembl identifier is not quickly accessible, use --feature-column 2 for the second column, which is usually more ubiquitous. Useful for overlaying gene expression so you can say --draw-leaf \"DrawItem (DrawContinuous \\\"CD4\\\")\") instead of --draw-leaf \"DrawItem (DrawContinuous \\\"ENSG00000010610\\\")\"). Does not affect CSV format (the column names will be the feature names)."-                  , normalization :: [String] <?> "([TfIdfNorm] | UQNorm | MedNorm | TotalMedNorm | TotalNorm | LogCPMNorm DOUBLE | QuantileNorm | NoneNorm) Type of normalization before clustering. Can be used as a list to perform one after the other: --normalization QuantileNorm --normalization TfIdfNorm will first apply quantile then tf-idf normalization. TfIdfNorm normalizes based on the prevalence of each feature. UQNorm normalizes each observation by the upper quartile non-zero counts of that observation. MedNorm normalizes each observation by the median non-zero counts of that observation. TotalNorm normalizes each cell by the total count. TotalMedNorm normalized first each observation by total count then by median of non-zero counts across features. (LogCPMNorm DOUBLE) normalizes by logB(CPM + 1) where B is DOUBLE. QuantileNorm normalizes by quantile normalization, ignores zeros, may be slow. NoneNorm does not normalize. Default is TfIdfNorm for clustering and NoneNorm for differential (which instead uses the recommended edgeR single cell preprocessing including normalization and filtering, any normalization provided here will result in edgeR preprocessing on top). Older versions had BothNorm which has been replaced with --normalization TotalMedNorm --normalization TfIdfNorm. This change will also affect other analyses, as TfIdfNorm will now be in non-cluster-related entry points if specified, so --normalization UQNorm from < v2.0.0.0 is now --normalization UQNorm --normalization TfIdfNorm."-                  , pca :: Maybe Int <?> "([Nothing] | INT) Not recommended, as it makes cosine similarity less meaningful (therefore less accurate -- instead, consider making your own similarity matrix and using cluster-tree, our sister algorithm, to cluster the matrix and plot with birch-beer). The number of dimensions to keep for PCA dimensionality reduction before clustering. Default is no PCA at all in order to keep all information. Should use with --shift-positive to ensure no negative values (as --pca will center and scale). Consider changing the modularity cutoff to a lower value (such as --min-modularity -0.5)."-                  , lsa :: Maybe Int <?> "([Nothing] | INT) The number of dimensions to keep for LSA dimensionality reduction. Uses TD-IDF followed by SVD before clustering, same warnings as --pca apply, including the use of --shift-positive with possible --min-modularity -0.5."-                  , svd :: Maybe Int <?> "([Nothing] | INT) The number of dimensions to keep for SVD dimensionality reduction. Will center and scale, same warnings as --pca apply, including the use of --shift-positive with possible --min-modularity -0.5."-                  , dropDimension :: Bool <?> "Instead of keeping dimensions for --pca, --lsa, or --svd, drop the number of dimensions (i.e. --lsa 1 drops the first projection which may have a strong batch effect)."-                  , shiftPositive :: Bool <?> "Shift features to positive values. Positive values are shifted to allow modularity to work correctly."-                  , filterThresholds :: Maybe String <?> "([Nothing] | (DOUBLE, DOUBLE)) The minimum filter thresholds for (MINCELL, MINFEATURE) when filtering cells and features by low read counts. Default changed to Nothing due to additional supported assays. To use the original filter thresholds, use --filter-thresholds \"(250, 1)\"."-                  , prior :: Maybe String <?> "([Nothing] | STRING) The input folder containing the output from a previous run. If specified, skips clustering by using the previous clustering files."-                  , noUpdateTreeRows :: Bool <?> "Don't update the row indices of a tree if using an identical matrix to the one which generated the tree. Should not be used unless the matrix to make the tree is identical, then can result in speedup."      }-    | Differential { matrixPath :: [String] <?> "(PATH) The path to the input directory containing the matrix output of cellranger (cellranger < 3 (matrix.mtx, genes.tsv, and barcodes.tsv) or cellranger >= 3 (matrix.mtx.gz, features.tsv.gz, and barcodes.tsv.gz) or an input csv file containing feature row names and cell column names. scATAC-seq is supported if input file contains \"fragments\", ends with \".tsv.gz\" (such as \"fragments.tsv.gz\" or \"sample1_fragments.tsv.gz\"), and is in the SORTED (sort -k1,1 -k2,2n) 10x fragments format (see also --binwidth, --no-binarize). If given as a list (--matrix-path input1 --matrix-path input2 etc.) then will join all matrices together. Assumes the same number and order of features in each matrix, so only cells are added."-                   , matrixTranspose :: Bool <?> "Whether to transpose the matrix before all processing (observations become features and vice-versa). Will be affected by other options (check your filtering thresholds, normalizations, etc!)"-                   , binwidth :: Maybe Int <?> "(Nothing | BINSIZE) If input data has region features in the format of `chrN:START-END`, BINSIZE input is required to convert ranges to fixed width bins."-                   , noBinarize :: Bool <?> "If a fragments.tsv.gz file, do not binarize data."-                   , binarize :: Bool <?> "Binarize data. Default for fragments.tsv.gz."-                   , excludeMatchFragments :: Maybe String <?> "([Nothing] | STRING) Exclude fragments from fragments.tsv.gz file if STRING is infix of fragment row. For instance, exclude all chrY fragments with \"--exclude-match-fragments chrY\"."-                   , blacklistRegionsFile :: Maybe T.Text <?> "([Nothing] | FILE) Bed file containing regions to ignore. Any fragments overlapping these regions are ignored if the input is not a fragments file."-                   , customRegion :: [T.Text] <?> "([Nothing] | chrN:START-END) Only look at these regions in the matrix for chromosome region features. A list in the format `chrN:START-END`, will assign values for a cell intersecting this region to this region. If a cell region overlaps two or more regions, will be counting that many times in the new features. Example input: `--custom-regions \"chr1:1003-10064\" --custom-regions \"chr2:3021-5034\"` etc."-                   , cellWhitelistFile :: Maybe String <?> "([Nothing] | FILE) The input file containing the cells to include. No header, line separated list of barcodes."-                   , labelsFile :: Maybe String <?> "([Nothing] | FILE) The input file containing the label for each cell barcode, with \"item,label\" header."-                   , customLabel :: [T.Text] <?> "([] | [LABEL]) List of labels to assign each matrix if all cells from each matrix are given the same label per matrix. This argument intends to simplify the process of labeling by bypassing --labels-file if the user just wants each matrix to have its own label (i.e. sample). Must be the same length and order as --matrix-path: for instance, --matrix-path input1 --custom-label sample1 --matrix-path input2 --custom-label sample2 etc. will label all cells from input1 with sample1, input2 with sample2, etc. If there are multiple labels per matrix, you must use --labels-file."-                   , featureColumn :: Maybe Int <?> "([1] | COLUMN) The column (1-indexed) in the features.tsv.gz file to use for feature names. If using matrix market format, cellranger stores multiple columns in the features file, usually the first column for the Ensembl identifier and the second column for the gene symbol. If the Ensembl identifier is not quickly accessible, use --feature-column 2 for the second column, which is usually more ubiquitous. Useful for overlaying gene expression so you can say --draw-leaf \"DrawItem (DrawContinuous \\\"CD4\\\")\") instead of --draw-leaf \"DrawItem (DrawContinuous \\\"ENSG00000010610\\\")\"). Does not affect CSV format (the column names will be the feature names)."-                   , pca :: Maybe Int <?> "([Nothing] | INT) Not recommended, as it makes cosine similarity less meaningful (therefore less accurate -- instead, consider making your own similarity matrix and using cluster-tree, our sister algorithm, to cluster the matrix and plot with birch-beer). The number of dimensions to keep for PCA dimensionality reduction before clustering. Default is no PCA at all in order to keep all information. Should use with --shift-positive to ensure no negative values (as --pca will center and scale). Consider changing the modularity cutoff to a lower value (such as --min-modularity -0.5)."-                   , lsa :: Maybe Int <?> "([Nothing] | INT) The number of dimensions to keep for LSA dimensionality reduction. Uses TD-IDF followed by SVD before clustering, same warnings as --pca apply, including the use of --shift-positive with possible --min-modularity -0.5."-                   , svd :: Maybe Int <?> "([Nothing] | INT) The number of dimensions to keep for SVD dimensionality reduction. Will center and scale, same warnings as --pca apply, including the use of --shift-positive with possible --min-modularity -0.5."-                   , dropDimension :: Bool <?> "Instead of keeping dimensions for --pca, --lsa, or --svd, drop the number of dimensions (i.e. --lsa 1 drops the first projection which may have a strong batch effect)."-                   , shiftPositive :: Bool <?> "Shift features to positive values. Positive values are shifted to allow modularity to work correctly."-                   , filterThresholds :: Maybe String <?> "([Nothing] | (DOUBLE, DOUBLE)) The minimum filter thresholds for (MINCELL, MINFEATURE) when filtering cells and features by low read counts. Default changed to Nothing due to additional supported assays. To use the original filter thresholds, use --filter-thresholds \"(250, 1)\"."-                   , delimiter :: Maybe Char <?> "([,] | CHAR) The delimiter for the most csv files in the program. For instance, if using a normal csv rather than cellranger output and for --labels-file."-                   , normalization :: [String] <?> "([TfIdfNorm] | UQNorm | MedNorm | TotalMedNorm | TotalNorm | LogCPMNorm DOUBLE | QuantileNorm | NoneNorm) Type of normalization before clustering. Can be used as a list to perform one after the other: --normalization QuantileNorm --normalization TfIdfNorm will first apply quantile then tf-idf normalization. TfIdfNorm normalizes based on the prevalence of each feature. UQNorm normalizes each observation by the upper quartile non-zero counts of that observation. MedNorm normalizes each observation by the median non-zero counts of that observation. TotalNorm normalizes each cell by the total count. TotalMedNorm normalized first each observation by total count then by median of non-zero counts across features. (LogCPMNorm DOUBLE) normalizes by logB(CPM + 1) where B is DOUBLE. QuantileNorm normalizes by quantile normalization, ignores zeros, may be slow. NoneNorm does not normalize. Default is TfIdfNorm for clustering and NoneNorm for differential (which instead uses the recommended edgeR single cell preprocessing including normalization and filtering, any normalization provided here will result in edgeR preprocessing on top). Older versions had BothNorm which has been replaced with --normalization TotalMedNorm --normalization TfIdfNorm. This change will also affect other analyses, as TfIdfNorm will now be in non-cluster-related entry points if specified, so --normalization UQNorm from < v2.0.0.0 is now --normalization UQNorm --normalization TfIdfNorm."-                   , prior :: Maybe String <?> "([Nothing] | STRING) The input folder containing the output from a previous run. If specified, skips clustering by using the previous clustering files."-                   , noUpdateTreeRows :: Bool <?> "Don't update the row indices of a tree if using an identical matrix to the one which generated the tree. Should not be used unless the matrix to make the tree is identical, then can result in speedup."-                   , noEdger :: Bool <?> "Use Kruskall-Wallis instead of edgeR for differential expression between two sets of nodes (automatically on and required for all to all comparisons)."-                   , nodes :: String <?> "([NODE], [NODE]) Find the differential expression between cells belonging downstream of a list of nodes versus another list of nodes. Directionality is \"([1], [2])\" -> 2 / 1. \"([], [])\" switches the process to instead find the log2 average division between all nodes with all other cells in the provided data set matrix regardless of whether they are present within the tree (node / other cells) using the Kruskal-Wallis Test (--features does not work for this, --labels works, and UQNorm for the normalization is recommended. Only returns nodes where the comparison had both groups containing at least five cells.). If not using --no-update-tree-rows, remember to filter the matrix for cells outside of the tree if you only want to compare against other nodes within the tree."-                   , labels :: Maybe String <?> "([Nothing] | ([LABEL], [LABEL])) Use --labels-file to restrict the differential analysis to cells with these labels. Same format as --nodes, so the first list in --nodes and --labels gets the cells within that list of nodes with this list of labels. The same for the second list. For instance, --nodes \"([1], [2])\" --labels \"([\\\"A\\\"], [\\\"B\\\"])\" will compare cells from node 1 of label \"A\" only with cells from node 2 of label \"B\" only. To use all cells for that set of nodes, use an empty list, i.e. --labels \"([], [\\\"A\\\"])\". When comparing all nodes with all other cells, remember that the notation would be ([Other Cells], [Node]), so to compare cells of label X in Node with cells of label Y in Other Cells, use --labels \"([\\\"Y\\\", \\\"X\\\"])\". Requires both --labels and --labels-file, otherwise will include all labels."-                   , topN :: Maybe Int <?> "([100] | INT ) The top INT differentially expressed features."-                   , seed :: Maybe Int <?> "([0] | INT) The seed to use for subsampling. See --subsample-groups."-                   , subsampleGroups :: Maybe Int <?> "([Nothing] | INT) Whether to subsample each group in the differential comparison. Subsets the specified number of cells. When set to 0, subsamples the larger group to equal the size of the smaller group. When using with --nodes \"([], [])\" to compare all nodes against each other, note that the compared nodes may be resampled. Highly experimental at this stage, use with caution. See --seed."-                   , features :: [T.Text] <?> "([Nothing] | FEATURE) List of features (e.g. genes) to plot for all cells within selected nodes. Invoked by --features CD4 --features CD8 etc. When this argument is supplied, only the plot is outputted and edgeR differential expression is ignored. Outputs to --output."-                   , aggregate :: Bool <?> "([False] | True) Whether to plot the aggregate (mean here) of features for each cell from \"--features\" instead of plotting different distributions for each feature."-                   , plotSeparateNodes :: Bool <?> "([False] | True) Whether to plot each node separately. This will plot each node provided in --nodes from both entries in the tuple (as they may be different from --labels)."-                   , plotSeparateLabels :: Bool <?> "([False] | True) Whether to plot each label separately. This will plot each label provided in --labels from both entries in the tuple (as they may be different from --nodes)."-                   , plotViolin :: Bool <?> "([False] | True) Whether to plot features as a violin plots instead of boxplots."-                   , plotNoOutlier :: Bool <?> "([False] | True) Whether to avoid plotting outliers as there can be too many, making the plot too large."-                   , plotOutput :: Maybe String <?> "([out.pdf] | STRING) The file containing the output plot."}-    | Diversity { priors :: [String] <?> "(PATH) Either input folders containing the output from a run of too-many-cells or a csv files containing the clusters for each cell in the format \"cell,cluster\". Advanced features not available in the latter case. If --labels-file is specified, those labels designate entity type, otherwise the assigned cluster is the entity type."-                , delimiter :: Maybe Char <?> "([,] | CHAR) The delimiter for the most csv files in the program. For instance, if using a normal csv rather than cellranger output and for --labels-file."-                , labelsFile :: Maybe String <?> "([Nothing] | FILE) The input file containing the label for each cell barcode, with \"item,label\" header."-                , start :: Maybe Integer <?> "([0] | INT) For the rarefaction curve, start the curve at this subsampling."-                , interval :: Maybe Integer <?> "([1] | INT) For the rarefaction curve, the amount to increase each subsampling. For instance, starting at 0 with an interval of 4, we would sampling 0, 4, 8, 12, ..."-                , end :: Maybe Integer <?> "([N] | INT) For the rarefaction curve, which subsample to stop at. By default, the curve stops at the observed number of species for each population."-                , order :: Maybe Double <?> "([1] | DOUBLE) The order of diversity."-                , output :: Maybe String <?> "([out] | STRING) The folder containing output."}-    | Paths { prior :: Maybe String <?> "([Nothing] | STRING) The input folder containing the output from a previous run. If specified, skips clustering by using the previous clustering files."-            , labelsFile :: Maybe String <?> "([Nothing] | FILE) The input file containing the label for each cell barcode, with \"item,label\" header."-            , flipDirection :: Bool <?> "Flip the starting node when calculating the distances."-            , shallowStart :: Bool <?> "Choose the shallowest leaf as the starting node (rather than the deepest)."-            , pathDistance :: Maybe String <?> "([PathStep] | PathModularity) How to measure the distance from the starting leaf. PathModularity weighs the steps by the modularity, while PathStep counts the number of steps."-            , bandwidth :: Maybe Double <?> "([1] | DOUBLE) Bandwidth of the density plot."-            , delimiter :: Maybe Char <?> "([,] | CHAR) The delimiter for the most csv files in the program. For instance, if using a normal csv rather than cellranger output and for --labels-file."-            , pathsPalette :: Maybe String <?> "([Set1] | Hsv | Ryb | Blues) Palette to use for legend colors."-            , output :: Maybe String <?> "([out] | STRING) The folder containing output."}-    | Classify { matrixPath :: [String] <?> "(PATH) The path to the input directory containing the matrix output of cellranger (cellranger < 3 (matrix.mtx, genes.tsv, and barcodes.tsv) or cellranger >= 3 (matrix.mtx.gz, features.tsv.gz, and barcodes.tsv.gz) or an input csv file containing feature row names and cell column names. scATAC-seq is supported if input file contains \"fragments\", ends with \".tsv.gz\" (such as \"fragments.tsv.gz\" or \"sample1_fragments.tsv.gz\"), and is in the SORTED (sort -k1,1 -k2,2n) 10x fragments format (see also --binwidth, --no-binarize). If given as a list (--matrix-path input1 --matrix-path input2 etc.) then will join all matrices together. Assumes the same number and order of features in each matrix, so only cells are added."-               , referenceFile :: [String] <?> "(PATH) The path to the reference file to compare each cell to. Every transformation (e.g. filters and normalizations) applied to --matrix-path apply here as well."-               , singleReferenceMatrix :: Bool <?> "Treat the reference file as a single matrix such that each observation (barcode) is an aggregated reference population."-               , matrixTranspose :: Bool <?> "Whether to transpose the matrix before all processing (observations become features and vice-versa). Will be affected by other options (check your filtering thresholds, normalizations, etc!)"-               , binwidth :: Maybe Int <?> "(Nothing | BINSIZE) If input data has region features in the format of `chrN:START-END`, BINSIZE input is required to convert ranges to fixed width bins."-               , noBinarize :: Bool <?> "If a fragments.tsv.gz file, do not binarize data."-               , binarize :: Bool <?> "Binarize data. Default for fragments.tsv.gz."-               , excludeMatchFragments :: Maybe String <?> "([Nothing] | STRING) Exclude fragments from fragments.tsv.gz file if STRING is infix of fragment row. For instance, exclude all chrY fragments with \"--exclude-match-fragments chrY\"."-               , blacklistRegionsFile :: Maybe T.Text <?> "([Nothing] | FILE) Bed file containing regions to ignore. Any fragments overlapping these regions are ignored if the input is not a fragments file."-               , customRegion :: [T.Text] <?> "([Nothing] | chrN:START-END) Only look at these regions in the matrix for chromosome region features. A list in the format `chrN:START-END`, will assign values for a cell intersecting this region to this region. If a cell region overlaps two or more regions, will be counting that many times in the new features. Example input: `--custom-regions \"chr1:1003-10064\" --custom-regions \"chr2:3021-5034\"` etc."-               , cellWhitelistFile :: Maybe String <?> "([Nothing] | FILE) The input file containing the cells to include. No header, line separated list of barcodes."-               , labelsFile :: Maybe String <?> "([Nothing] | FILE) The input file containing the label for each cell barcode, with \"item,label\" header."-               , customLabel :: [T.Text] <?> "([] | [LABEL]) List of labels to assign each matrix if all cells from each matrix are given the same label per matrix. This argument intends to simplify the process of labeling by bypassing --labels-file if the user just wants each matrix to have its own label (i.e. sample). Must be the same length and order as --matrix-path: for instance, --matrix-path input1 --custom-label sample1 --matrix-path input2 --custom-label sample2 etc. will label all cells from input1 with sample1, input2 with sample2, etc. If there are multiple labels per matrix, you must use --labels-file."-               , featureColumn :: Maybe Int <?> "([1] | COLUMN) The column (1-indexed) in the features.tsv.gz file to use for feature names. If using matrix market format, cellranger stores multiple columns in the features file, usually the first column for the Ensembl identifier and the second column for the gene symbol. If the Ensembl identifier is not quickly accessible, use --feature-column 2 for the second column, which is usually more ubiquitous. Useful for overlaying gene expression so you can say --draw-leaf \"DrawItem (DrawContinuous \\\"CD4\\\")\") instead of --draw-leaf \"DrawItem (DrawContinuous \\\"ENSG00000010610\\\")\"). Does not affect CSV format (the column names will be the feature names)."-               , pca :: Maybe Int <?> "([Nothing] | INT) Not recommended, as it makes cosine similarity less meaningful (therefore less accurate -- instead, consider making your own similarity matrix and using cluster-tree, our sister algorithm, to cluster the matrix and plot with birch-beer). The number of dimensions to keep for PCA dimensionality reduction before clustering. Default is no PCA at all in order to keep all information. Should use with --shift-positive to ensure no negative values (as --pca will center and scale). Consider changing the modularity cutoff to a lower value (such as --min-modularity -0.5)."-               , lsa :: Maybe Int <?> "([Nothing] | INT) The number of dimensions to keep for LSA dimensionality reduction. Uses TD-IDF followed by SVD before clustering, same warnings as --pca apply, including the use of --shift-positive with possible --min-modularity -0.5."-               , svd :: Maybe Int <?> "([Nothing] | INT) The number of dimensions to keep for SVD dimensionality reduction. Will center and scale, same warnings as --pca apply, including the use of --shift-positive with possible --min-modularity -0.5."-               , dropDimension :: Bool <?> "Instead of keeping dimensions for --pca, --lsa, or --svd, drop the number of dimensions (i.e. --lsa 1 drops the first projection which may have a strong batch effect)."-               , shiftPositive :: Bool <?> "Shift features to positive values. Positive values are shifted to allow modularity to work correctly."-               , filterThresholds :: Maybe String <?> "([Nothing] | (DOUBLE, DOUBLE)) The minimum filter thresholds for (MINCELL, MINFEATURE) when filtering cells and features by low read counts. Default changed to Nothing due to additional supported assays. To use the original filter thresholds, use --filter-thresholds \"(250, 1)\"."-               , delimiter :: Maybe Char <?> "([,] | CHAR) The delimiter for the most csv files in the program. For instance, if using a normal csv rather than cellranger output and for --labels-file."-               , normalization :: [String] <?> "([TfIdfNorm] | UQNorm | MedNorm | TotalMedNorm | TotalNorm | LogCPMNorm DOUBLE | QuantileNorm | NoneNorm) Type of normalization before clustering. Can be used as a list to perform one after the other: --normalization QuantileNorm --normalization TfIdfNorm will first apply quantile then tf-idf normalization. TfIdfNorm normalizes based on the prevalence of each feature. UQNorm normalizes each observation by the upper quartile non-zero counts of that observation. MedNorm normalizes each observation by the median non-zero counts of that observation. TotalNorm normalizes each cell by the total count. TotalMedNorm normalized first each observation by total count then by median of non-zero counts across features. (LogCPMNorm DOUBLE) normalizes by logB(CPM + 1) where B is DOUBLE. QuantileNorm normalizes by quantile normalization, ignores zeros, may be slow. NoneNorm does not normalize. Default is TfIdfNorm for clustering and NoneNorm for differential (which instead uses the recommended edgeR single cell preprocessing including normalization and filtering, any normalization provided here will result in edgeR preprocessing on top). Older versions had BothNorm which has been replaced with --normalization TotalMedNorm --normalization TfIdfNorm. This change will also affect other analyses, as TfIdfNorm will now be in non-cluster-related entry points if specified, so --normalization UQNorm from < v2.0.0.0 is now --normalization UQNorm --normalization TfIdfNorm." }-    | Peaks { fragmentsPath :: [String] <?> "(PATH) The path to the input fragments.tsv.gz file. The input file must contain \"fragments\" and end with \".tsv.gz\" (such as \"fragments.tsv.gz\" or \"sample1_fragments.tsv.gz\"), and is in the SORTED (sort -k1,1 -k2,2n) 10x fragments format (CHR\tSTART\tEND\tBARCODE\tVALUE). See --fragments-output from the make-tree and matrix-output entry points for assistance in merging files."-            , prior :: Maybe String <?> "([Nothing] | STRING) The input folder containing the output from a previous run. If specified, skips clustering by using the previous clustering files."-            , delimiter :: Maybe Char <?> "([,] | CHAR) The delimiter for the most csv files in the program. For instance, if using a normal csv rather than cellranger output and for --labels-file."-            , labelsFile :: Maybe String <?> "([Nothing] | FILE) The input file containing the label for each cell barcode, with \"item,label\" header."-            , excludeMatchFragments :: Maybe String <?> "([Nothing] | STRING) Exclude fragments from fragments.tsv.gz file if STRING is infix of fragment row. For instance, exclude all chrY fragments with \"--exclude-match-fragments chrY\"."-            , blacklistRegionsFile :: Maybe T.Text <?> "([Nothing] | FILE) Bed file containing regions to ignore. Any fragments overlapping these regions are ignored if the input is not a fragments file."-            , customRegion :: [T.Text] <?> "([Nothing] | chrN:START-END) Only look at these regions in the matrix for chromosome region features. A list in the format `chrN:START-END`, will assign values for a cell intersecting this region to this region. If a cell region overlaps two or more regions, will be counting that many times in the new features. Example input: `--custom-regions \"chr1:1003-10064\" --custom-regions \"chr2:3021-5034\"` etc."-            , peakCallCommand :: Maybe String <?> "([macs2 callpeak --nomodel --nolambda -p 0.001 -B -t %s -n %s --outdir %s] | STRING) The command to call peaks with. Can be any command that will be run on each generated fragment file per cluster, but the first \"%s\" must be the input argument, second \"%s\" is the name of the sample, and the third \"%s\" should be the output directory. Uses macs2 by default. Must return a .narrowPeak file with each row being \"CHR\tSTART\tEND\t*\tVALUE\n\" at least (* can be anything, after VALUE there can be anything as well. Check macs2 output for guidance)."-            , genomecovCommand :: Maybe String <?> "([bedtools genomecov -i %s -g %s -scale %f -bg -trackline > %s] | STRING) The command to convert to coverage bedgraph output. Can be any command that will be run on each bed per cluster, but the first \"%s\" must be the input argument, the second \"%s\" is the genome file (see https://github.com/arq5x/bedtools2/tree/master/genomes), followed by the \"%f\" scaling argument, with the last \"%s\" as the output argument, in order. Uses bedtools genomecov by default."-            , genome :: Maybe String <?> "([./human.hg38.genome] | PATH) The location of the genome file for the --genomecov-command, see https://github.com/arq5x/bedtools2/tree/master/genomes"-            , skipFragments :: Bool <?> "Whether to skip the generation of the fragments (e.g. if changing only --peak-call-command and fragment separation by cluster already exists)."-            , peakNode :: [Int] <?> "([ALLLEAFNODES] | NODE) List of nodes to peak call, i.e. \"--peak-node 3 --peak-node 5 --peak-node 7\". If the node is not a leaf node, make sure to use --all-nodes in addition."-            , peakNodeLabels :: [String] <?> "([ALLLABELS] | NODE) List of labels to keep in each node when outputting fragments and peaks, i.e. --peak-node-labels \"(3, [\\\"Red\\\"])\" --peak-node-labels \"(5, [\\\"Red\\\", \\\"Blue\\\"]. Nodes not listed will include all labels."-            , allNodes :: Bool <?> "Whether to get fragments and peaks for all nodes, not just the leaves."-            , bedgraph :: Bool <?> "Whether to output cluster normalized per million bedgraph output."-            , output :: Maybe String <?> "([out] | STRING) The folder containing output." }-    | Motifs { diffFile :: T.Text <?> "(FILE) The input file containing the differential features between nodes. Must be in the format `node,feature,log2FC,pVal,FDR`. The node column is optional (if wanting to separate per node)."-             , backgroundDiffFile :: Maybe T.Text <?> "(FILE) The input file containing the differential features between nodes for use as a background in motif finding. Must be in the format `node,feature,log2FC,pVal,FDR`. The node column is optional (if wanting to separate per node). If using this argument, be sure to update the --motif-command appropriately (background file comes last, e.g. with homer use `/path/to/findMotifs.pl %s fasta %s -bgFasta %s`)."-             , motifGenome :: T.Text <?> "(FILE) The location of the genome file in fasta format to convert bed to fasta."-             , motifCommand :: Maybe String <?> "([meme %s -nmotifs 50 -oc %s] | STRING) The command to find motifs in a fasta file. Can be any command that will be run on each fasta file converted from the bed optionally per node, but the first \"%s\" must be the input file, the second \"%s\" is the output. An example of homer: `/path/to/findMotifs.pl %s fasta %s`. Uses meme by default."-             , motifGenomeCommand :: Maybe String <?> "([Nothing] | STRING) The command to find motifs from a bed file instead of --motif-command (replaces that argument), as in homer's findMotifsGenome.pl. Can be any command that will be run on each bed file optionally per node, but the first \"%s\" must be the input file, the second \"%s\" is the genome file, and the last is the output. An example of homer: `/path/to/findMotifsGenome.pl %s %s %s`."-             , topN :: Maybe Int <?> "([100] | INT ) The top INT differentially expressed features."-             , output :: Maybe String <?> "([out] | STRING) The folder containing output." }-    | MatrixOutput { matrixPath :: [String] <?> "(PATH) The path to the input directory containing the matrix output of cellranger (cellranger < 3 (matrix.mtx, genes.tsv, and barcodes.tsv) or cellranger >= 3 (matrix.mtx.gz, features.tsv.gz, and barcodes.tsv.gz) or an input csv file containing feature row names and cell column names. scATAC-seq is supported if input file contains \"fragments\", ends with \".tsv.gz\" (such as \"fragments.tsv.gz\" or \"sample1_fragments.tsv.gz\"), and is in the SORTED (sort -k1,1 -k2,2n) 10x fragments format (see also --binwidth, --no-binarize). If given as a list (--matrix-path input1 --matrix-path input2 etc.) then will join all matrices together. Assumes the same number and order of features in each matrix, so only cells are added."-                   , binwidth :: Maybe Int <?> "(Nothing | BINSIZE) If input data has region features in the format of `chrN:START-END`, BINSIZE input is required to convert ranges to fixed width bins."-                   , noBinarize :: Bool <?> "If a fragments.tsv.gz file, do not binarize data."-                   , binarize :: Bool <?> "Binarize data. Default for fragments.tsv.gz."-                   , excludeMatchFragments :: Maybe String <?> "([Nothing] | STRING) Exclude fragments from fragments.tsv.gz file if STRING is infix of fragment row. For instance, exclude all chrY fragments with \"--exclude-match-fragments chrY\"."-                   , blacklistRegionsFile :: Maybe T.Text <?> "([Nothing] | FILE) Bed file containing regions to ignore. Any fragments overlapping these regions are ignored if the input is not a fragments file."-                   , customRegion :: [T.Text] <?> "([Nothing] | chrN:START-END) Only look at these regions in the matrix for chromosome region features. A list in the format `chrN:START-END`, will assign values for a cell intersecting this region to this region. If a cell region overlaps two or more regions, will be counting that many times in the new features. Example input: `--custom-regions \"chr1:1003-10064\" --custom-regions \"chr2:3021-5034\"` etc."-                   , cellWhitelistFile :: Maybe String <?> "([Nothing] | FILE) The input file containing the cells to include. No header, line separated list of barcodes."-                   , customLabel :: [T.Text] <?> "([] | [LABEL]) List of labels to assign each matrix if all cells from each matrix are given the same label per matrix. This argument intends to simplify the process of labeling by bypassing --labels-file if the user just wants each matrix to have its own label (i.e. sample). Must be the same length and order as --matrix-path: for instance, --matrix-path input1 --custom-label sample1 --matrix-path input2 --custom-label sample2 etc. will label all cells from input1 with sample1, input2 with sample2, etc. If there are multiple labels per matrix, you must use --labels-file."-                   , delimiter :: Maybe Char <?> "([,] | CHAR) The delimiter for the most csv files in the program. For instance, if using a normal csv rather than cellranger output and for --labels-file."-                   , featureColumn :: Maybe Int <?> "([1] | COLUMN) The column (1-indexed) in the features.tsv.gz file to use for feature names. If using matrix market format, cellranger stores multiple columns in the features file, usually the first column for the Ensembl identifier and the second column for the gene symbol. If the Ensembl identifier is not quickly accessible, use --feature-column 2 for the second column, which is usually more ubiquitous. Useful for overlaying gene expression so you can say --draw-leaf \"DrawItem (DrawContinuous \\\"CD4\\\")\") instead of --draw-leaf \"DrawItem (DrawContinuous \\\"ENSG00000010610\\\")\"). Does not affect CSV format (the column names will be the feature names)."-                   , normalization :: [String] <?> "([TfIdfNorm] | UQNorm | MedNorm | TotalMedNorm | TotalNorm | LogCPMNorm DOUBLE | QuantileNorm | NoneNorm) Type of normalization before clustering. Can be used as a list to perform one after the other: --normalization QuantileNorm --normalization TfIdfNorm will first apply quantile then tf-idf normalization. TfIdfNorm normalizes based on the prevalence of each feature. UQNorm normalizes each observation by the upper quartile non-zero counts of that observation. MedNorm normalizes each observation by the median non-zero counts of that observation. TotalNorm normalizes each cell by the total count. TotalMedNorm normalized first each observation by total count then by median of non-zero counts across features. (LogCPMNorm DOUBLE) normalizes by logB(CPM + 1) where B is DOUBLE. QuantileNorm normalizes by quantile normalization, ignores zeros, may be slow. NoneNorm does not normalize. Default is TfIdfNorm for clustering and NoneNorm for differential (which instead uses the recommended edgeR single cell preprocessing including normalization and filtering, any normalization provided here will result in edgeR preprocessing on top). Older versions had BothNorm which has been replaced with --normalization TotalMedNorm --normalization TfIdfNorm. This change will also affect other analyses, as TfIdfNorm will now be in non-cluster-related entry points if specified, so --normalization UQNorm from < v2.0.0.0 is now --normalization UQNorm --normalization TfIdfNorm."-                   , matrixTranspose :: Bool <?> "Whether to transpose the matrix before all processing (observations become features and vice-versa). Will be affected by other options (check your filtering thresholds, normalizations, etc!)"-                   , pca :: Maybe Int <?> "([Nothing] | INT) Not recommended, as it makes cosine similarity less meaningful (therefore less accurate -- instead, consider making your own similarity matrix and using cluster-tree, our sister algorithm, to cluster the matrix and plot with birch-beer). The number of dimensions to keep for PCA dimensionality reduction before clustering. Default is no PCA at all in order to keep all information. Should use with --shift-positive to ensure no negative values (as --pca will center and scale). Consider changing the modularity cutoff to a lower value (such as --min-modularity -0.5)."-                   , lsa :: Maybe Int <?> "([Nothing] | INT) The number of dimensions to keep for LSA dimensionality reduction. Uses TD-IDF followed by SVD before clustering, same warnings as --pca apply, including the use of --shift-positive with possible --min-modularity -0.5."-                   , svd :: Maybe Int <?> "([Nothing] | INT) The number of dimensions to keep for SVD dimensionality reduction. Will center and scale, same warnings as --pca apply, including the use of --shift-positive with possible --min-modularity -0.5."-                   , dropDimension :: Bool <?> "Instead of keeping dimensions for --pca, --lsa, or --svd, drop the number of dimensions (i.e. --lsa 1 drops the first projection which may have a strong batch effect)."-                   , filterThresholds :: Maybe String <?> "([Nothing] | (DOUBLE, DOUBLE)) The minimum filter thresholds for (MINCELL, MINFEATURE) when filtering cells and features by low read counts. Default changed to Nothing due to additional supported assays. To use the original filter thresholds, use --filter-thresholds \"(250, 1)\"."-                   , shiftPositive :: Bool <?> "Shift features to positive values. Positive values are shifted to allow modularity to work correctly."-                   , matOutput :: String <?> "([out_matrix] | FOLDER | FILE.csv) Output the filtered and normalized (not including TfIdfNorm) matrix in this folder in matrix market format or, if a csv file is specified, a dense csv format. Like input, features are rows." }-    deriving (Generic)+data Subcommand+    = MakeTreeCommand MakeTree+    | InteractiveCommand Interactive+    | DifferentialCommand Differential+    | DiversityCommand Diversity+    | PathsCommand Paths+    | ClassifyCommand Classify+    | PeaksCommand Peaks+    | MotifsCommand Motifs+    | MatrixOutputCommand MatrixOutput+    | SpatialCommand Spatial+    deriving (Read, Show) -modifiers :: Modifiers-modifiers = lispCaseModifiers { shortNameModifier = short }-  where-    short "allNodes"              = Nothing-    short "aggregate"             = Nothing-    short "atac"                  = Nothing-    short "binarize"              = Nothing-    short "bedgraph"              = Nothing-    short "blacklistRegionsFile"  = Nothing-    short "customRegion"          = Nothing-    short "clumpinessMethod"      = Just 'u'-    short "clusterNormalization"  = Just 'C'-    short "customCut"             = Nothing-    short "customLabel"           = Just 'Z'-    short "dendrogramOutput"      = Just 'U'-    short "diffFile"              = Nothing-    short "drawCollection"        = Just 'E'-    short "drawColors"            = Just 'R'-    short "drawDendrogram"        = Just 'D'-    short "drawDiscretize"        = Nothing-    short "drawFont"              = Nothing-    short "drawLeaf"              = Just 'L'-    short "drawLegendAllLabels"   = Just 'J'-    short "drawLegendSep"         = Just 'Q'-    short "drawMark"              = Just 'K'-    short "drawMaxLeafNodeSize"   = Nothing-    short "drawMaxNodeSize"       = Just 'A'-    short "drawNoScaleNodes"      = Just 'W'-    short "drawNodeNumber"        = Just 'N'-    short "drawPalette"           = Just 'Y'-    short "drawScaleSaturation"   = Just 'V'-    short "drawBarBounds"         = Nothing-    short "dropDimension"         = Nothing-    short "eigenGroup"            = Just 'B'-    short "elbowCutoff"           = Nothing-    short "featureColumn"         = Nothing-    short "filterThresholds"      = Just 'H'-    short "fragmentsOutput"       = Nothing-    short "fragmentsPath"         = Just 'f'-    short "genome"                = Nothing-    short "genomecovCommand"      = Nothing-    short "labels"                = Nothing-    short "labelsOutput"          = Nothing-    short "lsa"                   = Nothing-    short "matOutput"             = Nothing-    short "matrixOutput"          = Nothing-    short "matrixTranspose"       = Just 'T'-    short "maxDistance"           = Just 't'-    short "maxProportion"         = Just 'X'-    short "maxStep"               = Just 'S'-    short "minDistance"           = Nothing-    short "minDistanceSearch"     = Nothing-    short "minModularity"         = Nothing-    short "minSize"               = Just 'M'-    short "motifCommand"          = Nothing-    short "motifGenome"           = Nothing-    short "noBinarize"            = Nothing-    short "noEdger"               = Nothing-    short "normalization"         = Just 'z'-    short "numEigen"              = Just 'G'-    short "numRuns"               = Nothing-    short "order"                 = Just 'O'-    short "pathsPalette"          = Nothing-    short "pca"                   = Nothing-    short "peakNode"              = Nothing-    short "peakNodeLabels"        = Nothing-    short "plotOutput"            = Nothing-    short "plotSeparateNodes"     = Nothing-    short "plotSeparateLabels"    = Nothing-    short "plotViolin"            = Nothing-    short "priors"                = Just 'P'-    short "projectionFile"        = Just 'j'-    short "rootCut"               = Nothing-    short "seed"                  = Nothing-    short "shallowStart"          = Nothing-    short "shiftPositive"         = Nothing-    short "singleReferenceMatrix" = Nothing-    short "subsampleGroups"       = Nothing-    short "svd"                   = Nothing-    short "noUpdateTreeRows"      = Nothing-    short x                       = firstLetter x+sub :: Parser Subcommand+sub = hsubparser+  ( commandGroup "Analyses using the single-cell matrix"+ <> command "make-tree" (info makeTreeParser ( progDesc "Generate and plot the too-many-cells tree" ))+ <> command "interactive" (info interactiveParser ( progDesc "Interactive tree plotting (legacy, slow)" ))+ <> command "differential" (info differentialParser ( progDesc "Find differential features between groups of nodes" ))+ <> command "classify" (info classifyParser ( progDesc "Classify single-cells based on reference profiles" ))+ <> command "spatial" (info spatialParser ( progDesc "Spatially analyze single-cells" ))+ <> command "matrix-output" (info matrixOutputParser ( progDesc "Transform the input matrix only" ))+ ) <|> hsubparser ( commandGroup "No single-cell matrix analyses"+ <> command "diversity" (info diversityParser ( progDesc "Quantify the diversity and rarefaction curves of the tree" ))+ <> command "paths" (info pathsParser ( progDesc "Infer pseudo-time information from the tree" ))+  ) <|> hsubparser (commandGroup "too-many-peaks analyses for scATAC-seq"+ <> command "peaks" (info peaksParser ( progDesc "Find peaks in nodes for scATAC-seq" ))+ <> command "motifs" (info motifsParser ( progDesc "Find motifs from peaks for scATAC-seq" ))+  ) -instance ParseRecord Options where-    parseRecord = parseRecordWithModifiers modifiers+-- Global parsers+delimiterParser :: Parser Char+delimiterParser =+  option auto ( long "delimiter" <> metavar "CHAR" <> value ',' <> showDefault <> help "The delimiter for the most csv files in the program. For instance, if using a normal csv rather than cellranger output and for --labels-file." )++outputParser :: Parser String+outputParser = option str ( long "output" <> short 'o' <> metavar "FOLDER" <> value "out" <> showDefault <> help "The folder containing output." )++projectionParser :: Parser (Maybe String)+projectionParser = optional $ option str ( long "projection-file" <> short 'j' <> metavar "FILE" <> help "The input file containing positions of each cell for plotting. Format is \"item,x,y,sample\" as the header, with the sample column being optional (to separate, for instance, slides in a spatial analysis so there is no proximity between slides). Useful for 10x where a t-NSE projection is generated in \"projection.csv\". Cells without projections will not be plotted. If not supplied, no plot will be made." )++labelsFileParser :: Parser (Maybe String)+labelsFileParser = optional $ option str ( long "labels-file" <> short 'l' <> metavar "FILE" <> help "The input file containing the label for each cell barcode, with \"item,label\" header." )++data LoadMatrixOptions = LoadMatrixOptions {matrixPath :: [String]+                               , binwidth :: Maybe Int+                               , noBinarize :: Bool+                               , binarize :: Bool+                               , excludeMatchFragments :: Maybe String+                               , blacklistRegionsFile :: Maybe T.Text+                               , customRegion :: [T.Text]+                               , cellWhitelistFile :: Maybe String+                               , customLabel :: [T.Text]+                               , delimiter :: Char+                               , featureColumn :: Int+                               , normalization :: [NormType]+                               , matrixTranspose :: Bool+                               , pca :: Maybe Int+                               , lsa :: Maybe Int+                               , svd :: Maybe Int+                               , dropDimension :: Bool+                               , filterThresholds :: Maybe String+                               , shiftPositive :: Bool+                               } deriving (Read, Show)++loadMatrixParser :: Parser LoadMatrixOptions+loadMatrixParser = do+  matrixPath <- many $ option str ( long "matrix-path" <> short 'm' <> metavar "PATH" <> help "The path to the input directory containing the matrix output of cellranger (cellranger < 3 (matrix.mtx, genes.tsv, and barcodes.tsv) or cellranger >= 3 (matrix.mtx.gz, features.tsv.gz, and barcodes.tsv.gz) or an input csv file containing feature row names and cell column names. scATAC-seq is supported if input file contains \"fragments\", ends with \".tsv.gz\" (such as \"fragments.tsv.gz\" or \"sample1_fragments.tsv.gz\"), and is in the SORTED (sort -k1,1 -k2,2n) 10x fragments format (see also --binwidth, --no-binarize). If given as a list (--matrix-path input1 --matrix-path input2 etc.) then will join all matrices together. Assumes the same number and order of features in each matrix, so only cells are added." )+  binwidth <- optional $ option auto ( long "binwidth" <> short 'b' <> metavar "BINSIZE" <> help "If input data has region features in the format of `chrN:START-END`, BINSIZE input is required to convert ranges to fixed width bins." )+  noBinarize <- switch ( long "no-binarize" <> help "If a fragments.tsv.gz file, do not binarize data." )+  binarize <- switch ( long "binarize" <> help "Binarize data. Default for fragments.tsv.gz." )+  excludeMatchFragments <- optional $ option str ( long "exclude-match-fragments" <> short 'e' <> metavar "STRING" <> help "Exclude fragments from fragments.tsv.gz file if STRING is infix of fragment row. For instance, exclude all chrY fragments with \"--exclude-match-fragments chrY\"." )+  blacklistRegionsFile <- optional $ option str ( long "blacklist-regions-file" <> metavar "FILE" <> help "Bed file containing regions to ignore. Any fragments overlapping these regions are ignored if the input is not a fragments file." )+  customRegion <- many $ option str ( long "custom-region" <> metavar "chrN:START-END" <> help "Only look at these regions in the matrix for chromosome region features. A list in the format `chrN:START-END`, will assign values for a cell intersecting this region to this region. If a cell region overlaps two or more regions, will be counting that many times in the new features. Example input: `--custom-regions \"chr1:1003-10064\" --custom-regions \"chr2:3021-5034\"` etc." )+  cellWhitelistFile <- optional $ option str ( long "cell-whitelist-file" <> short 'c' <> metavar "FILE" <> help "The input file containing the cells to include. No header, line separated list of barcodes." )+  customLabel <- many $ option str ( long "customLabel" <> short 'Z' <> metavar "[LABEL]" <> help "List of labels to assign each matrix if all cells from each matrix are given the same label per matrix. This argument intends to simplify the process of labeling by bypassing --labels-file if the user just wants each matrix to have its own label (i.e. sample). Must be the same length and order as --matrix-path: for instance, --matrix-path input1 --custom-label sample1 --matrix-path input2 --custom-label sample2 etc. will label all cells from input1 with sample1, input2 with sample2, etc. If there are multiple labels per matrix, you must use --labels-file." )+  delimiter <- delimiterParser+  featureColumn <- option auto ( long "feature-column" <> metavar "COLUMN" <> value 1 <> showDefault <> help "The column (1-indexed) in the features.tsv.gz file to use for feature names. If using matrix market format, cellranger stores multiple columns in the features file, usually the first column for the Ensembl identifier and the second column for the gene symbol. If the Ensembl identifier is not quickly accessible, use --feature-column 2 for the second column, which is usually more ubiquitous. Useful for overlaying gene expression so you can say --draw-leaf \"DrawItem (DrawContinuous \\\"CD4\\\")\") instead of --draw-leaf \"DrawItem (DrawContinuous \\\"ENSG00000010610\\\")\"). Does not affect CSV format (the column names will be the feature names)." )+  normalization <- many $ option auto ( long "normalization" <> short 'z' <> metavar "[TfIdfNorm] | UQNorm | MedNorm | TotalMedNorm | TotalNorm | LogCPMNorm DOUBLE | QuantileNorm | NoneNorm" <> help "Type of normalization before clustering. Can be used as a list to perform one after the other: --normalization QuantileNorm --normalization TfIdfNorm will first apply quantile then tf-idf normalization. TfIdfNorm normalizes based on the prevalence of each feature. UQNorm normalizes each observation by the upper quartile non-zero counts of that observation. MedNorm normalizes each observation by the median non-zero counts of that observation. TotalNorm normalizes each cell by the total count. TotalMedNorm normalized first each observation by total count then by median of non-zero counts across features. (LogCPMNorm DOUBLE) normalizes by logB(CPM + 1) where B is DOUBLE. QuantileNorm normalizes by quantile normalization, ignores zeros, may be slow. MinMaxNorm does min-max normalization of each cell. TransposeNorm just transposes the matrix, so any row-based normalization here can instead be applied to columns: for instance, --normalization QuantileNorm --normalization TransposeNorm --normalization MinMaxNorm --normalization TransposeNorm will first apply quantile normalization to each cell, then min-max normalization to each column (before returning the cells to the proper axis with another transpose. If the number of TransposeNorms are odd, then one is automatically added to the end to ensure the same axis as before normalizations (see --matrix-transpose to permanently re-orientate the matrix)). NoneNorm does not normalize. Default is TfIdfNorm for clustering and NoneNorm for everything else (note, using --edger in differential analysis will result in edgeR preprocessing on top of any other normalization chosen). Older versions had BothNorm which has been replaced with --normalization TotalMedNorm --normalization TfIdfNorm. This change will also affect other analyses, as TfIdfNorm will now be in non-cluster-related entry points if specified, so --normalization UQNorm from < v2.0.0.0 is now --normalization UQNorm --normalization TfIdfNorm." )+  matrixTranspose <- switch ( long "matrix-transpose" <> short 'T' <> help "Whether to transpose the matrix before all processing (observations become features and vice-versa). Will be affected by other options (check your filtering thresholds, normalizations, etc!)" )+  pca <- optional $ option auto ( long "pca" <> metavar "INT" <> help "Not recommended, as it makes cosine similarity less meaningful (therefore less accurate -- instead, consider making your own similarity matrix and using cluster-tree, our sister algorithm, to cluster the matrix and plot with birch-beer). The number of dimensions to keep for PCA dimensionality reduction before clustering. Default is no PCA at all in order to keep all information. Should use with --shift-positive to ensure no negative values (as --pca will center and scale). Consider changing the modularity cutoff to a lower value (such as --min-modularity -0.5)." )+  lsa <- optional $ option auto ( long "lsa" <> metavar "INT" <> help "The number of dimensions to keep for LSA dimensionality reduction. Uses TD-IDF followed by SVD before clustering, same warnings as --pca apply, including the use of --shift-positive with possible --min-modularity -0.5." )+  svd <- optional $ option auto ( long "svd" <> metavar "INT" <> help "The number of dimensions to keep for SVD dimensionality reduction. Will center and scale, same warnings as --pca apply, including the use of --shift-positive with possible --min-modularity -0.5." )+  dropDimension <- switch ( long "drop-dimension" <> help "Instead of keeping dimensions for --pca, --lsa, or --svd, drop the number of dimensions (i.e. --lsa 1 drops the first projection which may have a strong batch effect)." )+  filterThresholds <- optional $ option str ( long "filter-thresholds" <> short 'H' <> metavar "(DOUBLE, DOUBLE)" <> help "The minimum filter thresholds for (MINCELL, MINFEATURE) when filtering cells and features by low read counts. Default changed to Nothing due to additional supported assays. To use the original filter thresholds, use --filter-thresholds \"(250, 1)\"." )+  shiftPositive <- switch ( long "shift-positive" <> help "Shift features to positive values. Positive values are shifted to allow modularity to work correctly." )++  pure $ LoadMatrixOptions {..}++data MakeTree = MakeTree { loadMatrixOptions :: LoadMatrixOptions+               , projectionFile :: Maybe String+               , labelsFile :: Maybe String+               , eigenGroup :: EigenGroup+               , numEigen :: Maybe Int+               , numRuns :: Maybe Int+               , minSize :: Maybe Int+               , maxStep :: Maybe Int+               , maxProportion :: Maybe Double+               , minModularity :: Maybe Double+               , minDistance :: Maybe Double+               , minDistanceSearch :: Maybe Double+               , smartCutoff :: Maybe Double+               , elbowCutoff :: Maybe String+               , customCut :: [Int]+               , rootCut :: Maybe Int+               , dendrogramOutput :: Maybe String+               , matrixOutput :: Maybe String+               , labelsOutput :: Bool+               , fragmentsOutput :: Bool+               , drawLeaf :: Maybe String+               , drawCollection :: DrawCollection+               , drawMark :: DrawNodeMark+               , drawNodeNumber :: Bool+               , drawMaxNodeSize :: Double+               , drawMaxLeafNodeSize :: Maybe Double+               , drawNoScaleNodes :: Bool+               , drawLegendSep :: Double+               , drawLegendAllLabels :: Bool+               , drawPalette :: Palette+               , drawColors :: Maybe String+               , drawDiscretize :: Maybe String+               , drawScaleSaturation :: Maybe Double+               , drawItemLineWeight :: Maybe Double+               , drawFont :: Maybe String+               , drawBarBounds :: Bool+               , prior :: Maybe String+               , noUpdateTreeRows :: Bool+               , order :: Double+               , clumpinessMethod :: Clump.Exclusivity+               , dense :: Bool+               , output :: String+               } deriving (Read, Show)++makeTreeParser :: Parser Subcommand+makeTreeParser = do+  loadMatrixOptions <- loadMatrixParser+  projectionFile <- projectionParser+  labelsFile <- labelsFileParser+  eigenGroup <- option auto ( long "eigen-group" <> short 'B' <> metavar "SignGroup | KMeansGroup" <> value SignGroup <> showDefault <> help "Whether to group the eigenvector using the sign or kmeans while clustering. While the default is sign, kmeans may be more accurate (but starting points are arbitrary)." )+  numEigen <- optional $ option auto ( long "num-eigen" <> short 'G' <> metavar "INT" <> help "Number of eigenvectors to use while clustering with kmeans. Takes from the second to last eigenvector. Recommended to start at 1 and work up from there if needed. May help offset the possible instability and inaccuracy of SVDLIBC." )+  numRuns <- optional $ option auto ( long "numRuns" <> metavar "([Nothing] | INT)" <> help "Number of runs for permutation test at each split for modularity. Defaults to no test." )+  minSize <- optional $ option auto ( long "min-size" <> short 'M' <> metavar "INT" <> help "The minimum size of a cluster. Defaults to 1." )+  maxStep <- optional $ option auto ( long "max-step" <> short 'S' <> metavar "INT" <> help "Only keep clusters that are INT steps from the root. Defaults to all steps." )+  maxProportion <- optional $ option auto ( long "max-proportion" <> short 'X' <> metavar "DOUBLE" <> help "Stopping criteria to stop at the node immediate after a node with DOUBLE proportion split. So a node N with L and R children will stop with this criteria at 0.5 if |L| / |R| < 0.5 or > 2 (absolute log2 transformed), that is, if one child has over twice as many items as the other child. Includes L and R in the final result." )+  minModularity <- optional $ option auto ( long "min-modularity" <> metavar "DOUBLE" <> help "Nearly the same as --min-distance, but for clustering instead of drawing (so the output json tree can be larger). Stopping criteria to stop at the node with DOUBLE modularity. So a node N with L and R children will stop with this criteria the distance at N to L and R is < DOUBLE. Does not include L and R in the final result." )+  minDistance <- optional $ option auto ( long "min-distance" <> short 't' <> metavar "DOUBLE" <> help "Stopping criteria to stop at the node immediate after a node with DOUBLE distance. So a node N with L and R children will stop with this criteria the distance at N to L and R is < DOUBLE. Includes L and R in the final result." )+  minDistanceSearch <- optional $ option auto ( long "min-distance-search" <> metavar "DOUBLE" <> help "Similar to --min-distance, but searches from the leaves to the root -- if a path from a subtree contains a distance of at least DOUBLE, keep that path, otherwise prune it. This argument assists in finding distant nodes." )+  smartCutoff <- optional $ option auto ( long "smart-cutoff" <> short 's' <> metavar "DOUBLE" <> help "Whether to set the cutoffs for --min-size, --max-proportion, --min-distance, and --min-distance-search based off of the distributions (median + (DOUBLE * MAD)) of all nodes. To use smart cutoffs, use this argument and then set one of the four arguments to an arbitrary number, whichever cutoff type you want to use. --max-proportion and --min-size distributions are log2 transformed." )+  elbowCutoff <- optional $ option str ( long "elbow-cutoff" <> metavar "Max | Min" <> help "Whether to set the cutoffs for --min-size, --max-proportion, --min-distance, and --min-distance-search based off of the elbow point of distributions of all nodes. For a distribution in positive x and y on a graph, the top left hump would be Max and the bottom right dip would be Min. To use elbow cutoffs, use this argument and then set one of the three arguments to an arbitrary number, whichever cutoff type you want to use. --max-proportion and --min-size distributions are log2 transformed. Conflicts with --smart-cutoff, so this argument takes precedent." )+  customCut <- many $ option auto ( long "custom-cut" <> metavar "NODE" <> help "List of nodes to prune (make these nodes leaves). Invoked by --custom-cut 34 --custom-cut 65 etc." )+  rootCut <- optional $ option auto ( long "root-cut" <> metavar "NODE" <> help "Assign a new root to the tree, removing all nodes outside of the subtree." )+  dendrogramOutput <- optional $ option str ( long "dendrogram-output" <> short 'U' <> metavar "FILE" <> value "dendrogram.svg" <> showDefault <> help "The filename for the dendrogram. Supported formats are PNG, PS, PDF, and SVG." )+  matrixOutput <- optional $ option str ( long "matrix-output" <> metavar "FOLDER | FILE.csv" <> help "Output the filtered and normalized (not including TfIdfNorm) matrix in this folder under the --output directory in matrix market format or, if a csv file is specified, a dense csv format. Like input, features are rows." )+  labelsOutput <- switch ( long "labels-output" <> help "Whether to write the labels used for each observation as a labels.csv file in the output folder." )+  fragmentsOutput <- switch ( long "fragments-output" <> help "Whether to output fragments_tsv.gz with barcodes altered by --custom-label in the output folder (excludes filtered-out cells). Useful for downstream analysis by the peaks entry point where the cluster barcodes differ from the original fragments.tsv.gz file when using --custom-label. Matches barcodes based on BARCODE-LABEL." )+  drawLeaf <- optional $ option str ( long "draw-leaf" <> short 'L' <> metavar "DrawText | DrawItem DrawItemType" <> help "How to draw leaves in the dendrogram. DrawText is the number of items in that leaf. DrawItem is the collection of items represented by circles, consisting of: DrawItem DrawLabel, where each item is colored by its label, DrawItem (DrawContinuous [FEATURE]), where each item is colored by the expression of FEATURE (corresponding to a feature name in the input matrix, [FEATURE] is a list, so if more than one FEATURE is listed, uses the average of the feature values), DrawItem (DrawThresholdContinuous [(FEATURE, THRESHOLD)]), where each item is colored by the binary high / low expression of FEATURE based on THRESHOLD (either `Exact DOUBLE` or `MadMedian DOUBLE`, where Exact just uses the DOUBLE as a cutoff value while MadMedian uses the DOUBLE as the number of MADs away from the median value of the feature) and multiple FEATUREs can be used to combinatorically label items (FEATURE1 high / FEATURE2 low, etc.), DrawItem DrawSumContinuous, where each item is colored by the sum of the post-normalized columns (use --normalization NoneNorm for UMI counts, default), DrawItem (DrawProximity ([NODE], DISTANCE)), where each item is colored by its proximity in Euclidean distance (a neighbor is defined as less than DISTANCE to the selected nodes) to a set of items drawn from a list of NODE, and DrawItem DrawDiversity, where each node is colored by the diversity based on the labels of each item and the color is normalized separately for the leaves and the inner nodes. The default is DrawText, unless --labels-file is provided, in which DrawItem DrawLabel is the default. If the label or feature cannot be found, the default color will be black (check your spelling!)." )+  drawCollection <- option auto ( long "draw-collection" <> short 'E' <> metavar "PieChart | PieRing | IndividualItems | Histogram | NoLeaf | CollectionGraph MAXWEIGHT THRESHOLD [NODE]" <> value PieChart <> showDefault <> help "How to draw item leaves in the dendrogram. PieRing draws a pie chart ring around the items. PieChart only draws a pie chart instead of items. IndividualItems only draws items, no pie rings or charts. Histogram plots a histogram of the features requested. NoLeaf has no leaf, useful if there are so many items the tree takes very long to render. (CollectionGraph MAXWEIGHT THRESHOLD [NODE]) draws the nodes and edges within leaves that are descendents of NODE (empty list [] indicates draw all leaf networks) based on the input matrix, normalizes edges based on the MAXWEIGHT, and removes edges for display less than THRESHOLD (after normalization, so for CollectionGraph 2 0.5 [26], draw the leaf graphs for all leaves under node 26, then a edge of 0.7 would be removed because (0.7 / 2) < 0.5). For CollectionGraph with no colors, use --draw-leaf \"DrawItem DrawLabel\" and all nodes will be black. If you don't specify this option, DrawText from --draw-leaf overrides this argument and only the number of cells will be plotted." )+  drawMark <- option auto ( long "draw-mark" <> short 'K' <> metavar "MarkNone | MarkModularity | MarkSignificance" <> value MarkNone <> showDefault <> help "How to draw annotations around each inner node in the tree. MarkNone draws nothing and MarkModularity draws a black circle representing the modularity at that node, darker black means higher modularity for that next split. MarkSignificance is for significance, i.e. p-value, darker means higher value." )+  drawNodeNumber <- switch ( long "draw-node-number" <> short 'N' <> help "Draw the node numbers on top of each node in the graph." )+  drawMaxNodeSize <- option auto ( long "draw-max-node-size" <> short 'A' <> metavar "DOUBLE" <> value 72 <> showDefault <> help "The max node size when drawing the graph. 36 is the theoretical default, but here 72 makes for thicker branches." )+  drawMaxLeafNodeSize <- optional $ option auto ( long "draw-max-leaf-node-size" <> metavar "DOUBLE" <> help "The max leaf node size when drawing the graph. Defaults to the value of --draw-max-node-size." )+  drawNoScaleNodes <- switch ( long "draw-no-scale-nodes" <> short 'W' <> help "Do not scale inner node size when drawing the graph. Instead, uses draw-max-node-size as the size of each node and is highly recommended to change as the default may be too large for this option." )+  drawLegendSep <- option auto ( long "draw-legend-sep" <> short 'Q' <> metavar "DOUBLE" <> value 1 <> showDefault <> help "The amount of space between the legend and the tree." )+  drawLegendAllLabels <- switch ( long "draw-legend-all-labels" <> short 'J' <> help "Whether to show all the labels in the label file instead of only showing labels within the current tree. The program generates colors from all labels in the label file first in order to keep consistent colors. By default, this value is false, meaning that only the labels present in the tree are shown (even though the colors are the same). The subset process occurs after --draw-colors, so when using that argument make sure to account for all labels." )+  drawPalette <- option auto ( long "draw-palette" <> short 'Y' <> metavar "Set1 | Hsv | Ryb | Blues" <> value Set1 <> showDefault <> help "Palette to use for legend colors. With high saturation in --draw-scale-saturation, consider using Hsv to better differentiate colors." )+  drawColors <- optional $ option str ( long "draw-colors" <> short 'R' <> metavar "COLORS" <> help "Custom colors for the labels or continuous features. Will repeat if more labels than provided colors. For continuous feature plots, uses first two colors [high, low], defaults to [red, gray]. For instance: --draw-colors \"[\\\"#e41a1c\\\", \\\"#377eb8\\\"]\"" )+  drawDiscretize <- optional $ option str ( long "draw-discretize" <> metavar "COLORS | INT" <> help "Discretize colors by finding the nearest color for each item and node. For instance, --draw-discretize \"[\\\"#e41a1c\\\", \\\"#377eb8\\\"]\" will change all node and item colors to one of those two colors, based on Euclidean distance. If using \"--draw-discretize INT\", will instead take the default map and segment (or interpolate) it into INT colors, rather than a more continuous color scheme. May have unintended results when used with --draw-scale-saturation." )+  drawScaleSaturation <- optional $ option auto ( long "draw-scale-saturation" <> short 'V' <> metavar "DOUBLE" <> help "Multiply the saturation value all nodes by this number in the HSV model. Useful for seeing more visibly the continuous colors by making the colors deeper against a gray scale." )+  drawItemLineWeight <- optional $ option auto ( long "draw-item-line-weight" <> metavar "DOUBLE" <> help "The line weight for items in the leaves if shown. Supplied as if there are too many items, the collection may look like a black box. Set to 0 to disable outlines of items to avoid this. Default: 0.1" )+  drawFont <- optional $ option str ( long "draw-font" <> metavar "FONT" <> help "Specify the font to use for the labels when plotting. Default: Arial" )+  drawBarBounds <- switch ( long "draw-bar-bounds" <> help "Whether to plot only the minimum and maximum ticks for the color bars." )+  prior <- optional $ option str ( long "prior" <> short 'p' <> metavar "STRING" <> help "The input folder containing the output from a previous run. If specified, skips clustering by using the previous clustering files." )+  noUpdateTreeRows <- switch ( long "no-update-tree-rows" <> help "Don't update the row indices of a tree if using an identical matrix to the one which generated the tree. Should not be used unless the matrix to make the tree is identical, then can result in speedup." )+  order <- option auto ( long "order" <> short 'O' <> metavar "DOUBLE" <> value 1 <> showDefault <> help "The order of diversity." )+  clumpinessMethod <- option auto ( long "clumpiness-method" <> short 'u' <> metavar "Majority | Exclusive | AllExclusive" <> value Clump.Majority <> showDefault <> help "The method used when calculating clumpiness: Majority labels leaves according to the most abundant label, Exclusive only looks at leaves consisting of cells solely from one label, and AllExclusive treats the leaf as containing both labels." )+  dense <- switch ( long "dense" <> help "Whether to use dense matrix algorithms for clustering. Should be faster for dense matrices, so if batch correction, PCA, or other algorithms are applied upstream to the input matrix, consider using this option to speed up the tree generation." )+  output <- outputParser++  pure $ MakeTreeCommand (MakeTree {..})++data Interactive = Interactive { loadMatrixOptions :: LoadMatrixOptions+                               , labelsFile :: Maybe String+                               , prior :: String+                               , noUpdateTreeRows :: Bool+                               } deriving (Read, Show)++interactiveParser :: Parser Subcommand+interactiveParser = do+  loadMatrixOptions <-  loadMatrixParser+  labelsFile <- labelsFileParser+  prior <- option str ( long "prior" <> metavar "STRING" <> help "The input folder containing the output from a previous run. If specified, skips clustering by using the previous clustering files." )+  noUpdateTreeRows <- switch ( long "no-update-tree-rows" <> help "Don't update the row indices of a tree if using an identical matrix to the one which generated the tree. Should not be used unless the matrix to make the tree is identical, then can result in speedup." )++  pure $ InteractiveCommand (Interactive {..})++data Differential = Differential { loadMatrixOptions :: LoadMatrixOptions+                   , labelsFile :: Maybe String+                   , prior :: String+                   , noUpdateTreeRows :: Bool+                   , edger :: Bool+                   , nodes :: String+                   , labels :: Maybe String+                   , topN :: Int+                   , seed :: Int+                   , subsampleGroups :: Maybe Int+                   , features :: [T.Text]+                   , aggregate :: Bool+                   , plotSeparateNodes :: Bool+                   , plotSeparateLabels :: Bool+                   , plotViolin :: Bool+                   , plotNoOutlier :: Bool+                   , plotOutput :: String+                   } deriving (Read, Show)++differentialParser :: Parser Subcommand+differentialParser = do+  loadMatrixOptions <-  loadMatrixParser+  labelsFile <- labelsFileParser+  prior <- option str ( long "prior" <> metavar "STRING" <> help "The input folder containing the output from a previous run. If specified, skips clustering by using the previous clustering files." )+  noUpdateTreeRows <- switch ( long "no-update-tree-rows" <> help "Don't update the row indices of a tree if using an identical matrix to the one which generated the tree. Should not be used unless the matrix to make the tree is identical, then can result in speedup." )+  edger <- switch ( long "edger" <> help "Use edgeR instead of Kruskall-Wallis for differential expression between two sets of nodes (automatically off and required to be off for all to all comparisons)." )+  nodes <- option str ( long "nodes" <> short 'n' <> metavar "([NODE], [NODE])" <> help "Find the differential expression between cells belonging downstream of a list of nodes versus another list of nodes. Directionality is \"([1], [2])\" -> 2 / 1. \"([], [])\" switches the process to instead find the log2 average division between all nodes with all other cells in the provided data set matrix regardless of whether they are present within the tree (node / other cells) using the Kruskal-Wallis Test (--features does not work for this, --labels works, and UQNorm for the normalization is recommended. Only returns nodes where the comparison had both groups containing at least five cells.). If not using --no-update-tree-rows, remember to filter the matrix for cells outside of the tree if you only want to compare against other nodes within the tree." )+  labels <- optional $ option str ( long "labels" <> metavar "([LABEL], [LABEL])" <> help "Use --labels-file to restrict the differential analysis to cells with these labels. Same format as --nodes, so the first list in --nodes and --labels gets the cells within that list of nodes with this list of labels. The same for the second list. For instance, --nodes \"([1], [2])\" --labels \"([\\\"A\\\"], [\\\"B\\\"])\" will compare cells from node 1 of label \"A\" only with cells from node 2 of label \"B\" only. To use all cells for that set of nodes, use an empty list, i.e. --labels \"([], [\\\"A\\\"])\". When comparing all nodes with all other cells, remember that the notation would be ([Other Cells], [Node]), so to compare cells of label X in Node with cells of label Y in Other Cells, use --labels \"([\\\"Y\\\", \\\"X\\\"])\". Requires both --labels and --labels-file, otherwise will include all labels." )+  topN <- option auto ( long "top-n" <> metavar "INT" <> value 100 <> showDefault <> help "The top INT differentially expressed features." )+  seed <- option auto ( long "seed" <> metavar "INT" <> value 0 <> showDefault <> help "The seed to use for subsampling. See --subsample-groups." )+  subsampleGroups <- optional $ option auto ( long "subsample-groups" <> metavar "INT" <> help "Whether to subsample each group in the differential comparison. Subsets the specified number of cells. When set to 0, subsamples the larger group to equal the size of the smaller group. When using with --nodes \"([], [])\" to compare all nodes against each other, note that the compared nodes may be resampled. Highly experimental at this stage, use with caution. See --seed." )+  features <- many $ option str ( long "features" <> metavar "FEATURE" <> help "List of features (e.g. genes) to plot for all cells within selected nodes. Invoked by --features CD4 --features CD8 etc. When this argument is supplied, only the plot is outputted and differential expression is ignored. Outputs to --output." )+  aggregate <- switch ( long "aggregate" <> help "Whether to plot the aggregate (mean here) of features for each cell from \"--features\" instead of plotting different distributions for each feature." )+  plotSeparateNodes <- switch ( long "plot-separate-nodes" <> help "Whether to plot each node separately. This will plot each node provided in --nodes from both entries in the tuple (as they may be different from --labels)." )+  plotSeparateLabels <- switch ( long "plot-separate-labels" <> help "Whether to plot each label separately. This will plot each label provided in --labels from both entries in the tuple (as they may be different from --nodes)." )+  plotViolin <- switch ( long "plot-violin" <> help "Whether to plot features as a violin plots instead of boxplots." )+  plotNoOutlier <- switch ( long "plot-no-outlier" <> help "Whether to avoid plotting outliers as there can be too many, making the plot too large." )+  plotOutput <- option str ( long "plot-output" <> metavar "FILE" <> value "out.pdf" <> showDefault <> help "The file containing the output plot." )++  pure $ DifferentialCommand (Differential {..})++data Diversity = Diversity { priors :: [String]+                , delimiter :: Char+                , labelsFile :: Maybe String+                , start :: Integer+                , interval :: Integer+                , end :: Maybe Integer+                , order :: Double+                , output :: String+                } deriving (Read, Show)++diversityParser :: Parser Subcommand+diversityParser = do+  priors <- many $ option str ( long "priors" <> short 'P' <> metavar "PATH" <> help "Either input folders containing the output from a run of too-many-cells or a csv files containing the clusters for each cell in the format \"cell,cluster\". Advanced features not available in the latter case. If --labels-file is specified, those labels designate entity type, otherwise the assigned cluster is the entity type." )+  delimiter <- delimiterParser+  labelsFile <- labelsFileParser+  start <- option auto ( long "start" <> short 's' <> metavar "INT" <> value 0 <> showDefault <> help "For the rarefaction curve, start the curve at this subsampling." )+  interval <- option auto ( long "interval" <> short 'i' <> metavar "INT" <> value 1 <> showDefault <> help "For the rarefaction curve, the amount to increase each subsampling. For instance, starting at 0 with an interval of 4, we would sampling 0, 4, 8, 12, ..." )+  end <- optional $ option auto ( long "end" <> metavar "INT" <> help "For the rarefaction curve, which subsample to stop at. By default, the curve stops at the observed number of species for each population." )+  order <- option auto ( long "order" <> metavar "DOUBLE" <> value 1 <> showDefault <> help "The order of diversity." )+  output <- outputParser++  pure $ DiversityCommand (Diversity {..})++data Paths = Paths { prior :: Maybe String+            , labelsFile :: Maybe String+            , flipDirection :: Bool+            , shallowStart :: Bool+            , pathDistance :: PathDistance+            , bandwidth :: Double+            , delimiter :: Char+            , pathsPalette :: Palette+            , output :: String+            } deriving (Read, Show)++pathsParser :: Parser Subcommand+pathsParser = do+  prior <- optional $ option str ( long "prior" <> metavar "FOLDER" <> help "The input folder containing the output from a previous run. If specified, skips clustering by using the previous clustering files." )+  labelsFile <- labelsFileParser+  flipDirection <- switch ( long "flip-direction" <> help "Flip the starting node when calculating the distances." )+  shallowStart <- switch ( long "shallow-start" <> help "Choose the shallowest leaf as the starting node (rather than the deepest)." )+  pathDistance <- option auto ( long "path-distance" <> metavar "PathStep | PathModularity" <> value PathStep <> showDefault <> help "How to measure the distance from the starting leaf. PathModularity weighs the steps by the modularity, while PathStep counts the number of steps." )+  bandwidth <- option auto ( long "bandwidth" <> metavar "DOUBLE" <> value 1 <> showDefault <> help "Bandwidth of the density plot." )+  delimiter <- delimiterParser+  pathsPalette <- option auto ( long "paths-palette" <> metavar "Set1 | Hsv | Ryb | Blues" <> value Set1 <> showDefault <> help "Palette to use for legend colors." )+  output <- outputParser++  pure $ PathsCommand (Paths {..})++data Classify = Classify { loadMatrixOptions :: LoadMatrixOptions+               , referenceFile :: [String]+               , singleReferenceMatrix :: Bool+               , skipAggregation :: Bool+               , labelsFile :: Maybe String+               } deriving (Read, Show)++classifyParser :: Parser Subcommand+classifyParser = do+  loadMatrixOptions <- loadMatrixParser+  referenceFile <- many $ option str (long "reference-file" <> short 'r' <> metavar "PATH" <> help "The path to the reference file to compare each cell to. Every transformation (e.g. filters and normalizations) applied to --matrix-path apply here as well." )+  singleReferenceMatrix <- switch ( long "single-reference-matrix" <> help "Treat the reference file as a single matrix such that each observation (barcode) is an aggregated reference population." )+  skipAggregation <- switch ( long "skip-aggregation" <> help "If there is more than one reference file, treat each reference file as a bulk population. Useful when comparing to bulk populations for BigWig or BedGraph reference files where each file is a separate population." )+  labelsFile <- labelsFileParser++  pure $ ClassifyCommand (Classify {..})++data Peaks = Peaks { fragmentsPath :: [String]+            , prior :: Maybe String+            , delimiter :: Char+            , labelsFile :: Maybe String+            , peaksExcludeMatchFragments :: Maybe String+            , peaksBlacklistRegionsFile :: Maybe T.Text+            , peaksCustomRegion :: [T.Text]+            , peakCallCommand :: String+            , genomecovCommand :: String+            , genome :: String+            , skipFragments :: Bool+            , peakNode :: [Int]+            , peakNodeLabels :: [String]+            , allNodes :: Bool+            , bedgraph :: Bool+            , output :: String+            } deriving (Read, Show)++peaksParser :: Parser Subcommand+peaksParser = do+  fragmentsPath <- many $ option str ( long "fragments-path" <> short 'f' <> metavar "PATH" <> help "The path to the input fragments.tsv.gz file. The input file must contain \"fragments\" and end with \".tsv.gz\" (such as \"fragments.tsv.gz\" or \"sample1_fragments.tsv.gz\"), and is in the SORTED (sort -k1,1 -k2,2n) 10x fragments format (CHR\tSTART\tEND\tBARCODE\tVALUE). See --fragments-output from the make-tree and matrix-output entry points for assistance in merging files." )+  prior <- optional $ option str ( long "prior" <> metavar "STRING" <> help "The input folder containing the output from a previous run. If specified, skips clustering by using the previous clustering files." )+  delimiter <- delimiterParser+  labelsFile <- labelsFileParser+  peaksExcludeMatchFragments <- optional $ option str ( long "peaks-exclude-match-fragments" <> metavar "STRING" <> help "Exclude fragments from fragments.tsv.gz file if STRING is infix of fragment row. For instance, exclude all chrY fragments with \"--exclude-match-fragments chrY\"." )+  peaksBlacklistRegionsFile <- optional $ option str ( long "peaks-blacklist-regions-file" <> metavar "FILE" <> help "Bed file containing regions to ignore. Any fragments overlapping these regions are ignored if the input is not a fragments file." )+  peaksCustomRegion <- many $ option str ( long "peaks-custom-region" <> metavar "chrN:START-END" <> help "Only look at these regions in the matrix for chromosome region features. A list in the format `chrN:START-END`, will assign values for a cell intersecting this region to this region. If a cell region overlaps two or more regions, will be counting that many times in the new features. Example input: `--custom-regions \"chr1:1003-10064\" --custom-regions \"chr2:3021-5034\"` etc." )+  peakCallCommand <- option str ( long "peak-call-command" <> metavar "STRING" <> value "macs2 callpeak --nomodel --nolambda -p 0.001 -B -t %s -n %s --outdir %s" <> showDefault <> help "The command to call peaks with. Can be any command that will be run on each generated fragment file per cluster, but the first \"%s\" must be the input argument, second \"%s\" is the name of the sample, and the third \"%s\" should be the output directory. Uses macs2 by default. Must return a .narrowPeak file with each row being \"CHR\tSTART\tEND\t*\tVALUE\n\" at least (* can be anything, after VALUE there can be anything as well. Check macs2 output for guidance)." )+  genomecovCommand <- option str ( long "genomecov-command" <> metavar "STRING" <> value "bedtools genomecov -i %s -g %s -scale %f -bg -trackline > %s" <> showDefault <> help "The command to convert to coverage bedgraph output. Can be any command that will be run on each bed per cluster, but the first \"%s\" must be the input argument, the second \"%s\" is the genome file (see https://github.com/arq5x/bedtools2/tree/master/genomes), followed by the \"%f\" scaling argument, with the last \"%s\" as the output argument, in order. Uses bedtools genomecov by default." )+  genome <- option str ( long "genome" <> metavar "PATH" <> value "./human.hg38.genome" <> showDefault <> help "The location of the genome file for the --genomecov-command, see https://github.com/arq5x/bedtools2/tree/master/genomes" )+  skipFragments <- switch ( long "skip-fragments" <> help "Whether to skip the generation of the fragments (e.g. if changing only --peak-call-command and fragment separation by cluster already exists)." )+  peakNode <- many $ option auto ( long "peak-node" <> metavar "[NODE]" <> help "List of nodes to peak call, i.e. \"--peak-node 3 --peak-node 5 --peak-node 7\". If the node is not a leaf node, make sure to use --all-nodes in addition. Defaults to all leaf nodes." )+  peakNodeLabels <- many $ option str ( long "peak-node-labels" <> metavar "[LABEL])" <> help "List of labels to keep in each node when outputting fragments and peaks, i.e. --peak-node-labels \"(3, [\\\"Red\\\"])\" --peak-node-labels \"(5, [\\\"Red\\\", \\\"Blue\\\"]. Nodes not listed will include all labels. Defaults to all labels." )+  allNodes <- switch ( long "all-nodes" <> help "Whether to get fragments and peaks for all nodes, not just the leaves." )+  bedgraph <- switch ( long "bedgraph" <> help "Whether to output cluster normalized per million bedgraph output." )+  output <- outputParser++  pure $ PeaksCommand (Peaks {..})++data Motifs = Motifs { diffFile :: T.Text+             , backgroundDiffFile :: Maybe T.Text+             , motifGenome :: T.Text+             , motifCommand :: String+             , motifGenomeCommand :: Maybe String+             , topN :: Int+             , output :: String+             } deriving (Read, Show)++motifsParser :: Parser Subcommand+motifsParser = do+  diffFile <- option str ( long "diff-file" <> metavar "FILE" <> help "The input file containing the differential features between nodes. Must be in the format `node,feature,log2FC,pVal,qVal`. The node column is optional (if wanting to separate per node)." )+  backgroundDiffFile <- optional $ option str ( long "background-diff-file" <> metavar "FILE" <> help "The input file containing the differential features between nodes for use as a background in motif finding. Must be in the format `node,feature,log2FC,pVal,qVal`. The node column is optional (if wanting to separate per node). If using this argument, be sure to update the --motif-command appropriately (background file comes last, e.g. with homer use `/path/to/findMotifs.pl %s fasta %s -bgFasta %s`)." )+  motifGenome <- option str ( long "motif-genome" <> metavar "FILE" <> help "The location of the genome file in fasta format to convert bed to fasta." )+  motifCommand <- option str ( long "motif-command" <> metavar "STRING" <> value "meme %s -nmotifs 50 -oc %s" <> showDefault <> help "The command to find motifs in a fasta file. Can be any command that will be run on each fasta file converted from the bed optionally per node, but the first \"%s\" must be the input file, the second \"%s\" is the output. An example of homer: `/path/to/findMotifs.pl %s fasta %s`. Uses meme by default." )+  motifGenomeCommand <- optional $ option str ( long "motif-genome-command" <> metavar "STRING" <> help "The command to find motifs from a bed file instead of --motif-command (replaces that argument), as in homer's findMotifsGenome.pl. Can be any command that will be run on each bed file optionally per node, but the first \"%s\" must be the input file, the second \"%s\" is the genome file, and the last is the output. An example of homer: `/path/to/findMotifsGenome.pl %s %s %s`." )+  topN <- option auto ( long "top-n" <> metavar "INT " <> value 100 <> showDefault <> help "The top INT differentially expressed features." )+  output <- outputParser++  pure $ MotifsCommand (Motifs {..})++data MatrixOutput = MatrixOutput { loadMatrixOptions :: LoadMatrixOptions+                   , matOutput :: String+                   } deriving (Read, Show)++matrixOutputParser :: Parser Subcommand+matrixOutputParser = do+  loadMatrixOptions <- loadMatrixParser+  matOutput <- option str ( long "matrix-output" <> metavar "FOLDER | FILE.csv" <> value "out_matrix" <> showDefault <> help "Output the filtered and normalized (not including TfIdfNorm) matrix in this folder in matrix market format or, if a csv file is specified, a dense csv format. Like input, features are rows." )++  pure $ MatrixOutputCommand (MatrixOutput {..})++data Spatial = Spatial { loadMatrixOptions :: LoadMatrixOptions+              , projectionFile :: Maybe String+              , labelsFile :: Maybe String+              , stateLabelsFile :: Maybe String+              , annoSpatMarkerFile :: Maybe String+              , annoSpatCommand :: String+              , mark :: [T.Text]+              , topDistances :: Maybe TopDistances+              , pcfCrossFlag :: Bool+              , onlySummaryFlag :: Bool+              , skipFinishedFlag :: Bool+              , includeOthersFlag :: Bool+              , output :: String+              } deriving (Read, Show)++spatialParser :: Parser Subcommand+spatialParser = do+  loadMatrixOptions <- loadMatrixParser+  labelsFile <- labelsFileParser+  projectionFile <- projectionParser+  stateLabelsFile <- optional $ option str ( long "state-labels-file" <> metavar "FILE" <> help "The input file containing a metadata label for sample (i.e. disease or control) to group samples by, with \"item,label\" header." )+  mark <- many $ option str ( long "mark" <> metavar "FEATURE | LABEL" <> help "Marks for the spatial relationship analysis. A list (`--mark MARK --mark MARK`) where `MARK` is either a feature such as a gene or protein or a label from the `--labels-file`. If the cells have labels (from `--labels-file`, `--annospat-marker-file`, or `--custom-label`), the `MARK` will be interpreted as a label. If `--mark ALL` (capitalized ALL, only ALL and nothing else) and `--labels-file` is given, then a pairwise comparison between labels will be computed, or if the `--labels-file` is not given then a pairwise comparison between features will be computed. If no marks are given, no spatial relationship will be carried out and only the interactive projection plot will be given. Importantly, if --include-others is enabled, the summary output files will also contain Kruskal-Wallis calculations with the shown distributions, but also with the marksOther-marksOther relationships to ensure the observation is not random. This explains why the Kruskal-Wallis value may be different from a single label vs. label comparison value, which would be equal otherwise." )+  topDistances <- optional $ option auto ( long "top-distances" <> metavar "[TopQuantile 10] | TopQuantile DOUBLE | TopDistance DOUBLE" <> help "For \"top\" (close) distance statistics, this number defines the top quantile (TopQuantile) or exact top distance (TopDistance). For instance, TopQuantile 4 would be the closest 25% of all distances, while TopQuantile 10 would be the closest 10% of all distances, and TopDistance 15 would be all distances from the origin to 15 (whatever unit that would be)." )+  annoSpatMarkerFile <- optional $ option str ( long "annospat-marker-file" <> metavar "STRING" <> help "The location of the AnnoSpat marker file for cell classification. Triggers the use of AnnoSpat to generate the labels file rather than --labels-file. Overrides the --labels-file argument." )+  annoSpatCommand <- option str ( long "annospat-command" <> metavar "STRING" <> value "AnnoSpat generateLabels -i %s -m %s -o %s -f %s -l %s -r %s -s \"\"" <> showDefault <> help "The AnnoSpat command to label cells. To customize, use the default value then follow the with custom additional arguments, making sure not to alter the default arguments used here." )+  pcfCrossFlag <- switch ( long "pcf-cross" <> help "Whether to use the multitype point correlation function (pcfcross) instead of the mark correlation function (markcrosscorr) for spatial relationships using categorical marks." )+  onlySummaryFlag <- switch ( long "summary-only" <> help "Whether to only calculate the summarization of the spatial statistics, for use only if skipping the spatial relationship calculations." )+  skipFinishedFlag <- switch ( long "skip-finished" <> help "Whether to skip the finished comparisons (contains the stats.csv file in the output folder)." )+  includeOthersFlag <- switch ( long "include-others" <> help "Whether to include marksOther as background, useful if there are no state labels." )+  output <- outputParser++  pure $ SpatialCommand (Spatial {..})

@@ -12,6 +12,7 @@ {-# LANGUAGE PackageImports    #-} {-# LANGUAGE TypeOperators     #-} {-# LANGUAGE TupleSections     #-}+{-# LANGUAGE DuplicateRecordFields #-}  module TooManyCells.Program.Paths where @@ -24,7 +25,6 @@ import Data.Maybe (fromMaybe) import Language.R as R import Math.Clustering.Hierarchical.Spectral.Types (getClusterItemsDend, EigenGroup (..))-import Options.Generic import Text.Read (readMaybe) import qualified Data.Aeson as A import qualified Data.ByteString.Lazy.Char8 as B@@ -44,32 +44,25 @@ import TooManyCells.Matrix.Types  -- | Paths path.-pathsMain :: Options -> IO ()-pathsMain opts = do+pathsMain :: Subcommand -> IO ()+pathsMain (PathsCommand opts) = do     let readOrErr err = fromMaybe (error err) . readMaybe         labelsFile'   =             maybe (error "\nNeed a label file.") LabelFile-                . unHelpful-                . labelsFile+                . (labelsFile :: Paths -> Maybe String)                 $ opts         prior'        =             maybe (error "\nNeed a prior path containing tree.") PriorPath-                . unHelpful-                . prior+                . (prior :: Paths -> Maybe String)                 $ opts-        delimiter'    =-            Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts-        bandwidth'    = Bandwidth . fromMaybe 1 . unHelpful . bandwidth $ opts-        direction'    = FlipFlag . unHelpful . flipDirection $ opts-        shallow'      = ShallowFlag . unHelpful . shallowStart $ opts-        palette'      = maybe Set1 ( readOrErr "Cannot read --paths-palette")-                      . unHelpful-                      . pathsPalette-                      $ opts-        pathDistance' =-            maybe PathStep read . unHelpful . pathDistance $ opts+        delimiter'    = Delimiter . (delimiter :: Paths -> Char) $ opts+        bandwidth'    = Bandwidth . bandwidth $ opts+        direction'    = FlipFlag . flipDirection $ opts+        shallow'      = ShallowFlag . shallowStart $ opts+        palette'      = pathsPalette opts+        pathDistance' = pathDistance opts         output'       =-            OutputDirectory . fromMaybe "out" . unHelpful . output $ opts+            OutputDirectory . (output :: Paths -> String) $ opts      -- Where to place output files.     FP.createDirectoryIfMissing True . unOutputDirectory $ output'@@ -100,3 +93,4 @@         return ()      return ()+pathsMain _ = error "Wrong path in paths, contact Gregory Schwartz for this error."

@@ -4,12 +4,13 @@ Peaks entry point for command line program. -} +{-# LANGUAGE DuplicateRecordFields #-}+ module TooManyCells.Program.Peaks where  -- Remote import BirchBeer.Load (loadLabelData) import BirchBeer.Types (Delimiter (..), LabelFile (..), Label (..))-import Options.Generic import Control.Monad (unless, when) import Data.Maybe (fromMaybe, isNothing) import System.IO (hPutStrLn, stderr)@@ -34,50 +35,38 @@ import TooManyCells.Program.Utility  -- | Peaks path.-peaksMain :: Options -> IO ()-peaksMain opts = do+peaksMain :: Subcommand -> IO ()+peaksMain (PeaksCommand opts) = do   let readOrErr err = fromMaybe (error err) . readMaybe-      fragmentPaths' = unHelpful . fragmentsPath $ opts+      fragmentPaths' = fragmentsPath opts       prior'    = PriorPath                 . fromMaybe (error "\nRequires a previous run to get the clusters.")-                . unHelpful-                . prior+                . (prior :: Peaks -> Maybe String)                 $ opts-      genome' =-        GenomeFile . fromMaybe "./human.hg38.genome" . unHelpful . genome $ opts-      genomecovCommand' = GenomecovCommand-                   . fromMaybe "bedtools genomecov -i %s -g %s -scale %f -bg > %s"-                   . unHelpful-                   . genomecovCommand-                   $ opts-      peakCommand' = PeakCommand-                   . fromMaybe "macs2 callpeak --nomodel --nolambda -p 0.001 -B -t %s -n %s --outdir %s"-                   . unHelpful-                   . peakCallCommand-                   $ opts-      delimiter'    = Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts-      labelsFile'   = fmap LabelFile . unHelpful . labelsFile $ opts-      bedgraphFlag' = BedGraphFlag . unHelpful . bedgraph $ opts-      allNodesFlag' = AllNodesFlag . unHelpful . allNodes $ opts-      peakNodes'    = PeakNodes . ISet.fromList . unHelpful . peakNode $ opts+      genome' = GenomeFile . genome $ opts+      genomecovCommand' = GenomecovCommand . genomecovCommand $ opts+      peakCommand' = PeakCommand . peakCallCommand $ opts+      delimiter'    = Delimiter . (delimiter :: Peaks -> Char) $ opts+      labelsFile'   = fmap LabelFile . (labelsFile :: Peaks -> Maybe String) $ opts+      bedgraphFlag' = BedGraphFlag . bedgraph $ opts+      allNodesFlag' = AllNodesFlag . allNodes $ opts+      peakNodes'    = PeakNodes . ISet.fromList . peakNode $ opts       peakNodesLabels' = PeakNodesLabels                        . IMap.fromList                        . fmap (L.over L._2 (Set.fromList . fmap Label))                        . fmap (readOrErr "\nCannot read --peak-nodes-labels")-                       . unHelpful                        . peakNodeLabels                        $ opts-      output' = OutputDirectory . fromMaybe "out" . unHelpful . output $ opts+      output' = OutputDirectory . (output :: Peaks -> String) $ opts       clInput = (FP.</> "cluster_list.json") . unPriorPath $ prior'       treeInput = (FP.</> "cluster_tree.json") . unPriorPath $ prior' -  when (isNothing . unHelpful . genome $ opts) $-    hPutStrLn stderr "--genome file not specified, using ./human.hg38.genome ..."+  hPutStrLn stderr ("Using --genome " <> genome opts <> " ...")    cr <- loadClusterResultsFiles clInput treeInput :: IO ClusterResults   labelMap <- mapM (loadLabelData delimiter') labelsFile' -  unless (unHelpful $ skipFragments opts) $ do+  unless (skipFragments opts) $ do     hPutStrLn stderr "Splitting fragment file by cluster ..."     mapM_       ( (=<<) ( saveClusterFragments@@ -97,3 +86,4 @@    hPutStrLn stderr "Calling peaks ..."   peakCallFiles peakCommand' genome' output' =<< mapM getMatrixFileType fragmentPaths'+peaksMain _ = error "Wrong path in peaks, contact Gregory Schwartz for this error."

@@ -0,0 +1,245 @@+{- TooManyCells.Program.Spatial+Gregory W. Schwartz++Spatial entry point for command line program.+-}++{-# LANGUAGE BangPatterns      #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE DuplicateRecordFields #-}+{-# LANGUAGE ViewPatterns #-}++module TooManyCells.Program.Spatial where++-- Remote+import BirchBeer.Load (loadLabelData)+import BirchBeer.Types (LabelFile (..), Delimiter (..), LabelMap (..), Sample (..), Feature (..), Label (..))+import Control.Concurrent.Async.Pool (withTaskGroup, mapTasks)+import Control.Monad (join, forM_, guard, unless, when, void)+import Data.Bool (bool)+import Data.Maybe (isJust, fromMaybe)+import GHC.Conc (getNumCapabilities)+import Language.R.Instance as R+import Language.R.QQ+import System.Directory (getTemporaryDirectory)+import System.IO (hPutStrLn, stderr)+import Text.Read (readMaybe)+import qualified Control.Foldl as Fold+import qualified Control.Lens as L+import qualified Data.Map.Strict as Map+import qualified Data.Set as Set+import qualified Data.Text as T+import qualified Data.Vector as V+import qualified H.Prelude as H+import qualified System.FilePath as FP+import qualified Turtle as TU++-- Local+import TooManyCells.Matrix.Load (loadProjectionMap)+import TooManyCells.Program.LoadMatrix (loadAllSSM)+import TooManyCells.Program.Options+import TooManyCells.Spatial.AnnoSpat (scToAnnoSpatFile, runAnnoSpat)+import TooManyCells.Spatial.ProjectionPlot (plotSpatialProjection)+import TooManyCells.Spatial.SummaryPlot (plotSummary)+import TooManyCells.Spatial.Relationships (spatialRelationshipsR)+import TooManyCells.Spatial.Utility (subsampleProjectionMap, markToText)+import qualified TooManyCells.File.Types as Too+import qualified TooManyCells.Matrix.Types as Too+import qualified TooManyCells.Program.Utility as Too+import qualified TooManyCells.Spatial.Types as Too++-- | General call for spatialRelationshipsR.+spatialRelationshipsCall :: Too.OutputDirectory+                         -> Too.SkipFinishedFlag+                         -> Too.PCFCrossFlag+                         -> Too.ProjectionMap+                         -> Maybe LabelMap+                         -> Too.SingleCells+                         -> Maybe Sample+                         -> Maybe Too.TopDistances+                         -> [Too.Mark]+                         -> IO ()+spatialRelationshipsCall (Too.OutputDirectory outDir) skipFlag' pcfCrossFlag' pm lm sc sample td marks = do+  let markString = T.unpack . T.intercalate "_" . fmap markToText $ marks+      folderName (Just (Sample x)) = T.unpack x <> "_" <> markString+      folderName Nothing = markString+      outDir' = Too.OutputDirectory+              $ outDir+         FP.</> folderName sample+      finalFile = TU.fromText+                . T.pack+                $ Too.unOutputDirectory outDir' FP.</> "stats.csv"++  -- If stats.csv exists and is not empty, then run everything+  !skip <- (&&)+       <$> (pure $ Too.unSkipFinishedFlag skipFlag')+       <*> (TU.testfile finalFile)+  when (not skip)+    $ spatialRelationshipsR outDir' pcfCrossFlag' pm lm sc td marks++-- | Helper for spatial relationships function.+spatialRelationshipsHelper :: Too.OutputDirectory+                           -> Too.SkipFinishedFlag+                           -> Too.PCFCrossFlag+                           -> Too.ProjectionMap+                           -> Maybe LabelMap+                           -> Too.SingleCells+                           -> Maybe Sample+                           -> Maybe Too.TopDistances+                           -> [Too.Mark]+                           -> IO ()+spatialRelationshipsHelper _ _ _ _ _ _ _ _ [] = pure ()+spatialRelationshipsHelper outDir skipFlag' pcfCrossFlag' pm lm sc sample td (fmap markToText -> ["ALL"]) = do+  cores <- getNumCapabilities+  withTaskGroup cores $ \workers ->+    void+      . mapTasks workers  -- mapTasks_ not working+      . fmap (\(!x, !y) -> spatialRelationshipsCall outDir skipFlag' pcfCrossFlag' pm lm sc sample td [x, y])+      $ (,) <$> getMarks lm <*> getMarks lm+  where+    getMarks Nothing = fmap Too.MarkFeature . V.toList . L.view Too.colNames $ sc+    getMarks (Just lm) =+      fmap Too.MarkLabel . Set.toList . Set.fromList . Map.elems . unLabelMap $ lm+spatialRelationshipsHelper outDir skipFlag' pcfCrossFlag' pm lm sc sample td marks =+  spatialRelationshipsCall outDir skipFlag' pcfCrossFlag' pm lm sc sample td marks++-- | Helper for plotting spatial summary.+spatialSummary :: Too.OutputDirectory+               -> Delimiter+               -> Too.IncludeOthersFlag+               -> Maybe Too.StateLabelsFile+               -> IO ()+spatialSummary outputDir' delimiter' iof stateLabelsFile' = do+  stateLabelMap <-+    mapM+      ( fmap (Too.StateLabelMap . unLabelMap)+      . loadLabelData delimiter'+      . LabelFile+      . Too.unStateLabelsFile+      )+      stateLabelsFile'+  let featureList = fmap+                      Feature+                      [ "meanCorr"+                      , "maxCorr"+                      , "minCorr"+                      , "topMaxCorr"+                      , "topMinCorr"+                      , "topMeanCorr"+                      , "negSwap"+                      , "posSwap"+                      , "longestPosLength"+                      , "longestNegLength"+                      , "maxPosWithVal"+                      , "logMaxPosWithVal"+                      , "maxPos"+                      , "minPos"+                      ]+  mapM_ (plotSummary outputDir' iof stateLabelMap) featureList++-- | Spatial path.+spatialMain :: Subcommand -> IO ()+spatialMain sub@(SpatialCommand opts) = H.withEmbeddedR R.defaultConfig $ do+  let onlySummaryFlag' = onlySummaryFlag opts+      outputDir' = Too.OutputDirectory . (output :: Spatial -> String) $ opts+      stateLabelsFile' = fmap Too.StateLabelsFile . stateLabelsFile $ opts+      delimiter' = Delimiter+                 . (delimiter :: LoadMatrixOptions -> Char)+                 . (loadMatrixOptions :: Spatial -> LoadMatrixOptions)+                 $ opts+      skipFlag' = Too.SkipFinishedFlag $ skipFinishedFlag opts+      iof' = Too.IncludeOthersFlag $ includeOthersFlag opts+      td' = topDistances opts++  unless onlySummaryFlag' $ do++    let readOrErr err = fromMaybe (error err) . readMaybe+        delimiter'        = Delimiter+                          . (delimiter :: LoadMatrixOptions -> Char)+                          . (loadMatrixOptions :: Spatial -> LoadMatrixOptions)+                          $ opts+        projectionFile' =+            maybe (error "--projection-file required") Too.ProjectionFile+              . (projectionFile :: Spatial -> Maybe String)+              $ opts+        annoSpatMarkerFile' = fmap Too.AnnoSpatMarkerFile . annoSpatMarkerFile $ opts+        annoSpatCommand' = Too.AnnoSpatCommand . annoSpatCommand $ opts++    scRes <- fmap (fromMaybe (error "Requires --matrix-path"))+          . loadAllSSM sub+          $ (loadMatrixOptions :: Spatial -> LoadMatrixOptions) opts+    let processedSc = fst scRes+        customLabelMap = snd scRes++    projectionMap <- loadProjectionMap projectionFile'++    let labelsFileError = hPutStrLn stderr "Warning: Problem in AnnoSpat, skipping label generation ..."+    labelsFile' <-+      case annoSpatMarkerFile' of+        Nothing ->+          pure . fmap LabelFile . (labelsFile :: Spatial -> Maybe String) $ opts+        (Just mf) -> (=<<) (maybe (labelsFileError >> pure Nothing) pure) . TU.reduce Fold.head $ do+          let annoOutDir = Too.OutputDirectory+                        $ Too.unOutputDirectory outputDir' FP.</> "AnnoSpat_out"+          tmpDir <- TU.liftIO $ fmap (TU.fromText . T.pack) getTemporaryDirectory+          tmpOut <- fmap Too.TempPath $ TU.mktempfile tmpDir "AnnoSpat_input.csv"+          startEndCols <- TU.liftIO $ scToAnnoSpatFile projectionMap processedSc tmpOut++          case startEndCols of+            Nothing -> pure Nothing+            (Just (startCol, endCol)) ->+              TU.liftIO+                $ runAnnoSpat+                    annoSpatCommand'+                    tmpOut+                    mf+                    annoOutDir+                    startCol+                    endCol++    labelMap <- if isJust labelsFile'+                  then mapM (loadLabelData delimiter') $ labelsFile'+                  else return customLabelMap++    let marks' = fmap+                  ( bool (Too.MarkFeature . Feature) (Too.MarkLabel . Label)+                  . isJust+                  $ labelMap+                  )+              . mark+              $ opts++    pcfCrossFlag' <- case (isJust labelMap, pcfCrossFlag opts) of+                      (False, True) -> do+                        hPutStrLn stderr "Warning: Continuous feature marks detected but pcfcross requested, ignoring pcfcross request ..."+                        return $ Too.PCFCrossFlag False+                      _ ->+                        return . Too.PCFCrossFlag . pcfCrossFlag $ opts++    let samples = Set.toList+                . Set.fromList+                . fmap fst+                . Map.elems+                . Too.unProjectionMap+                $ projectionMap++    forM_ samples $ \s -> do+      let s' = fromMaybe (Sample "total") s+          sOutLabel = T.unpack $ unSample s'+          projectionOutput = Too.OutputDirectory+                            . (FP.</> sOutLabel FP.</> "projections")+                            . Too.unOutputDirectory+                            $ outputDir'+          relationshipsOutput = Too.OutputDirectory+                              . (FP.</> sOutLabel FP.</> "relationships")+                              . Too.unOutputDirectory+                              $ outputDir'+          subPm = subsampleProjectionMap s projectionMap+      plotSpatialProjection projectionOutput labelMap subPm processedSc s'++      case marks' of+        [] -> pure ()+        m -> spatialRelationshipsHelper relationshipsOutput skipFlag' pcfCrossFlag' subPm labelMap processedSc s td' m++  spatialSummary outputDir' delimiter' iof' stateLabelsFile'+spatialMain _ = error "Wrong path in spatial, contact Gregory Schwartz for this error."

@@ -4,16 +4,21 @@ Utility functions for the command line program. -} +{-# LANGUAGE BangPatterns      #-}+ module TooManyCells.Program.Utility where  -- Remote+import Control.Exception (SomeException (..), try, evaluate)+import Control.Monad (guard)+import Control.Monad.Trans.Maybe (runMaybeT, MaybeT (..)) import Data.Bool (bool) import Data.List (isInfixOf) import Data.Maybe (fromMaybe)-import Options.Generic import Text.Read (readMaybe) import qualified System.Directory as FP import qualified System.FilePath as FP+import qualified Turtle as TU  -- Local import TooManyCells.File.Types@@ -25,25 +30,12 @@ readOrErr err = fromMaybe (error err) . readMaybe  -- | Normalization defaults.-getNormalization :: Options -> [NormType]-getNormalization opts@(MakeTree{}) =-  (\x -> if null x then [TfIdfNorm] else x)-    . fmap (readOrErr "Cannot read --normalization.")-    . unHelpful-    . normalization-    $ opts-getNormalization opts@(Interactive{}) =-  (\x -> if null x then [TfIdfNorm] else x)-    . fmap (readOrErr "Cannot read --normalization.")-    . unHelpful-    . normalization-    $ opts-getNormalization opts =-  (\x -> if null x then [NoneNorm] else x)-    . fmap (readOrErr "Cannot read --normalization.")-    . unHelpful-    . normalization-    $ opts+getNormalization :: Subcommand -> [NormType] -> [NormType]+getNormalization (MakeTreeCommand opts) xs =+  if null xs then [TfIdfNorm] else xs+getNormalization (InteractiveCommand opts) xs =+  if null xs then [TfIdfNorm] else xs+getNormalization _ xs = if null xs then [NoneNorm] else xs  -- | Get the file type of an input matrix. Returns either a left file (i.e. CSV) -- or a right matrix market.@@ -68,3 +60,11 @@         | otherwise = error "\nMatrix path does not exist."    return matrixFile++-- | Check if file exists and is not empty. Currently may have a handle closing+-- issue, do not use for now.+nonEmptyFileExists :: TU.FilePath -> IO Bool+nonEmptyFileExists path = either (const' False) (/= TU.B 0)+                      <$> (try $ TU.du path :: IO (Either SomeException TU.Size))+  where+    const' !x !y = x

@@ -0,0 +1,108 @@+{- TooManyCells.Spatial.AnnoSpat+Gregory W. Schwartz++Collects functions in the program for using AnnoSpat to label cells with classification.+-}++{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE TupleSections #-}++module TooManyCells.Spatial.AnnoSpat+    ( scToAnnoSpatFile+    , runAnnoSpat+    ) where++-- Remote+import BirchBeer.Types+import Control.Monad (liftM2, liftM3, join)+import Data.Bool (bool)+import Data.List (zipWith3, foldl')+import Data.Maybe (catMaybes, fromMaybe, isJust)+import TextShow (showt)+import TooManyCells.Matrix.Types+import qualified Control.Foldl as Fold+import qualified Control.Lens as L+import qualified Data.ByteString.Lazy.Char8 as BL+import qualified Data.Csv as CSV+import qualified Data.Map.Strict as Map+import qualified Data.Sparse.Common as S+import qualified Data.Text as T+import qualified Data.Text.Encoding as T+import qualified Data.Vector as V+import qualified Text.Printf as TP+import qualified Turtle as TU++-- Local+import qualified TooManyCells.File.Types as Too+import qualified TooManyCells.Spatial.Types as Too++-- | SingleCells to input file for AnnoSpat.+scToAnnoSpatFile ::+  ProjectionMap -> SingleCells -> Too.TempPath -> IO (Maybe (Too.StartCol, Too.EndCol))+scToAnnoSpatFile pm sc (Too.TempPath tmpOut) = do+  let header = V.fromList ["item", "sample"]+            <> (fmap (T.encodeUtf8 . unFeature) . L.view colNames $ sc)+      startEndCols = L.over L._2 Too.EndCol+                   . L.over L._1 Too.StartCol+                   . L.over L.both unFeature+                 <$> ( L.sequenceOf L.both+                     . Fold.fold ((,) <$> Fold.head <*> Fold.last)+                     . L.view colNames+                     $ sc+                     )+      itemOrdered =+        zip (repeat "item") . fmap unCell . V.toList . L.view rowNames $ sc+      sampleOrdered =+        fmap (\ !x -> ("sample",)+                    . maybe "NA" unSample+                    . join+                    . fmap fst+                    . Map.lookup x+                    $ unProjectionMap pm+             )+          . V.toList+          . L.view rowNames+          $ sc++  BL.writeFile (T.unpack . TU.format TU.fp $ tmpOut)+    . CSV.encodeByName header+    . zipWith3 (\x y zs -> Map.fromList $ x:y:zs) itemOrdered sampleOrdered+    . fmap (zip (V.toList . getColNames $ sc) . fmap showt . S.toDenseListSV)+    . S.toRowsL+    . getMatrix+    $ sc++  pure startEndCols++-- | Run AnnoSpat to get cell labels.+runAnnoSpat :: Too.AnnoSpatCommand+            -> Too.TempPath+            -> Too.AnnoSpatMarkerFile+            -> Too.OutputDirectory+            -> Too.StartCol+            -> Too.EndCol+            -> IO (Maybe LabelFile)+runAnnoSpat defCmd tmpInput mf outDir startCol endCol = TU.reduce Fold.head $ do+  let cmd = T.pack+          $ TP.printf+              (Too.unAnnoSpatCommand defCmd)+              (T.unpack . TU.format TU.fp $ Too.unTempPath tmpInput)+              (Too.unAnnoSpatMarkerFile mf)+              (Too.unOutputDirectory outDir)+              (T.unpack $ Too.unStartCol startCol)+              (T.unpack $ Too.unEndCol endCol)+              ("sample" :: String)++  TU.mktree . TU.fromText . T.pack $ Too.unOutputDirectory outDir+  TU.stderr . TU.inshell cmd $ mempty++  labelsFile <- fmap (LabelFile . T.unpack . TU.format TU.fp)+              . TU.single+              . TU.findtree (TU.invert $ TU.has "numericLabels")+              . TU.find (TU.has "trte_labels_")+              . TU.fromText+              . T.pack+              $ Too.unOutputDirectory outDir++  pure labelsFile

@@ -0,0 +1,264 @@+{- TooManyCells.Spatial.ProjectionPlot+Gregory W. Schwartz++Collects functions pertaining to plotting interactive figures of point proximity+on the left with cumulative distribution functions of features on the right.+-}++{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE DeriveGeneric #-}++module TooManyCells.Spatial.ProjectionPlot+    ( plotSpatialProjection+    ) where++-- Remote+import BirchBeer.Types (Feature (..), LabelMap (..), Label (..), getColNames, Id (..), Sample (..))+import Data.Bool (bool)+import Data.Colour.Palette.BrewerSet (Kolor, brewerSet, ColorCat (..) )+import Data.Colour.Palette.Harmony (colorRamp)+import Data.Colour.SRGB (sRGB24show)+import Data.List (genericLength)+import Data.Maybe (fromMaybe, isJust, catMaybes)+import qualified Control.Foldl as Fold+import qualified Control.Lens as L+import qualified Data.ByteString.Char8 as BS+import qualified Data.ByteString.Lazy.Char8 as B+import qualified Data.Foldable as F+import qualified Data.Map.Strict as Map+import qualified Data.Set as Set+import qualified Data.Text as T+import qualified Data.Text.Read as T+import qualified Data.Vector as V+import qualified Graphics.Vega.VegaLite as VL+import qualified Graphics.Vega.VegaLite.Theme as VL+import qualified System.FilePath as FP+import qualified Turtle as TU++-- Local+import TooManyCells.File.Types (OutputDirectory (..))+import TooManyCells.Matrix.Types (SingleCells (..), ProjectionMap (..), X (..), Y (..))+import TooManyCells.Spatial.Types (ColorMap (..), Range (..))+import TooManyCells.Spatial.Utility++-- | Get the minimum and maximum values for a projection map.+getMinMax :: ProjectionMap -> ((Double, Double), (Double, Double))+getMinMax =+  fromMaybe ((0, 0), (0, 0))+    . (L.sequenceOf L.both :: (Maybe (Double, Double), Maybe (Double, Double)) -> Maybe ((Double, Double), (Double, Double)))  -- m ((a, b), (c, d))+    . (L.over L.both (L.sequenceOf L.both) :: ((Maybe Double, Maybe Double), (Maybe Double, Maybe Double)) -> (Maybe (Double, Double), Maybe (Double, Double)))  -- (m (a, b), m (b, c))+    . (L.over L.both (Fold.fold ((,) <$> Fold.minimum <*> Fold.maximum)) :: ([Double], [Double]) -> ((Maybe Double, Maybe Double), (Maybe Double, Maybe Double))) -- ((m a, m b), (m c, m d))+    . (L.over (L._2 . L.mapped) unY :: ([Double], [Y]) -> ([Double], [Double]) )+    . L.over (L._1 . L.mapped) unX+    . unzip+    . fmap snd+    . Map.elems+    . unProjectionMap++-- | Get the color mapping from feature to color.+getColorMap :: [Feature] -> ColorMap+getColorMap xs =+  ColorMap . Map.fromList . zip xs . colorRamp (length xs) . brewerSet Set1 $ 9++-- | Get the feature selection of a distribution.+featureSel :: Feature -> [VL.SelectSpec] -> VL.PropertySpec+featureSel (Feature f) = VL.selection+                       . VL.select+                           ("pick_" <> f)+                           VL.Interval [VL.Encodings [VL.ChX], VL.Empty]++-- | Color by features.+colorByFeatures+  :: Maybe LabelMap -> ColorMap -> [Feature] -> VL.BuildEncodingSpecs+colorByFeatures lm colorMap features =+        VL.color ( maybe  -- Choose whether to color by labels or features+                    [ VL.MDataCondition+                        featureSels+                        [ VL.MString "white" ]+                    ]+                    (\ x+                    -> [ VL.MName "label", VL.MmType VL.Nominal, labelColorScale x]+                    )+                    lm+                  )+  where+    featureSels = fmap+                    (\ f+                    -> ( VL.Selection ("pick_" <> unFeature f)+                       ,  [ VL.MString+                              ( maybe "white" (T.pack . sRGB24show)+                              . Map.lookup f+                              $ unColorMap colorMap+                              )+                          ]+                        )+                    )+                    features+++-- | Pseudo-legend for selection.+legendSpec :: LabelMap -> VL.VLSpec+legendSpec lm =+  VL.asSpec [ VL.mark VL.Circle []+            , selection []+            , encoding []+            ]+  where+    selection = VL.selection+              . VL.select+                  "pick_legend"+                  VL.Multi+                  [ VL.Fields [ "label" ], VL.Empty ]+    labels = Set.toAscList+           . Set.fromList+           . fmap unLabel+           . Map.elems+           . unLabelMap+           $ lm+    encoding =+      VL.encoding+        . VL.position VL.Y [ VL.PName "label", VL.PmType VL.Nominal ]+        . VL.color [ VL.MName "label"+                   , VL.MmType VL.Nominal+                   , labelColorScale lm+                   ,  VL.MLegend []+                   ]++-- | Get the circle spec for plotting.+getCircleSpec :: Maybe LabelMap -> Range -> ColorMap -> [Feature] -> VL.VLSpec+getCircleSpec lm range colorMap features =+  VL.asSpec [ VL.mark VL.Circle []+            , VL.width plotWidth+            , VL.height plotHeight+            , circleEnc []+            , circleFilterTrans []+            ]+  where+    aspectRatio = (maxY range - minY range) / (maxX range - minX range)+    plotWidth = 800+    plotHeight = aspectRatio * plotWidth+    circleEnc =+      VL.encoding+        . VL.position VL.X [ VL.PName "x"+                            , VL.PmType VL.Quantitative+                            , VL.PAxis [ VL.AxTitle "X Axis" ]+                            , VL.PScale+                                [ VL.SDomain (VL.DNumbers [minX range, maxX range])+                                ]+                            ]+        . VL.position VL.Y [ VL.PName "y"+                            , VL.PmType VL.Quantitative+                            , VL.PAxis [ VL.AxTitle "Y Axis" ]+                            , VL.PScale+                                [ VL.SDomain (VL.DNumbers [minY range, maxY range])+                                ]+                            ]+        . colorByFeatures lm colorMap features+    circleFilterTrans =+      VL.transform+        . VL.filter+            ( VL.FCompose+                ( bool (VL.Expr "false") (VL.Selection "pick_legend") (isJust lm)+          `VL.Or` F.foldl'+                  (\acc x -> VL.Or acc (VL.Selection $ "pick_" <> unFeature x))+                  (VL.Expr "false")+                  features+                )+            )++-- | Get a window spec of a feature for plotting.+getWindowSpec :: ColorMap -> Feature -> VL.VLSpec+getWindowSpec colorMap (Feature feature) =+  VL.asSpec [ VL.title feature []+            , VL.height 40+            , windowEnc []+            , windowTrans []+            , featureSel (Feature feature) []+            , VL.mark VL.Area []+            ]+  where+    windowTrans =+      VL.transform+        . VL.filter (VL.FExpr $ "datum." <> feature <> " > 0")+        . VL.window+            [([VL.WOp VL.CumeDist, VL.WField feature], "window_" <> feature)]+            [VL.WFrame Nothing (Just 0), VL.WSort [VL.WAscending feature]]+    windowEnc = VL.encoding+               . VL.position VL.X [ VL.PName $  feature+                                  , VL.PmType VL.Quantitative+                                  , VL.PAxis [ VL.AxTitle "Expression"]+                                  ]+               . VL.position VL.Y [ VL.PName $ "window_" <> feature+                                  , VL.PmType VL.Quantitative+                                  , VL.PAxis [ VL.AxTitle "Probability"]+                                  ]+               . VL.color+                  [ VL.MString ( maybe "white" (T.pack . sRGB24show)+                               . Map.lookup (Feature feature)+                               $ unColorMap colorMap+                               )+                  ]++-- | The color scheme for the label field.+labelColorScale :: LabelMap -> VL.MarkChannel+labelColorScale lm = VL.MScale [ VL.SDomain (VL.DStrings labels)+                               , VL.SRange (VL.RStrings colors)+                               ]+  where+    labels =+      Set.toAscList . Set.fromList . fmap unLabel . Map.elems . unLabelMap $ lm+    colors =+      fmap (T.pack . sRGB24show) . colorRamp (length labels) . brewerSet Set1 $ 9++plotSpatialProjection :: OutputDirectory+                      -> Maybe LabelMap+                      -> ProjectionMap+                      -> SingleCells+                      -> Sample+                      -> IO ()+plotSpatialProjection outputDir' labelMap' pm sc sample = do+  let labelMap = LabelMap . Map.insert (Id "") (Label "NA") . unLabelMap+             <$> labelMap'  -- Insert dummy cell to account for NA labels.+      dataSet = scToVLData labelMap pm sc+      features = fmap Feature . V.toList . getColNames $ sc+      ((miX, maX), (miY, maY)) = getMinMax pm+      range = Range miX maX miY maY+      numWindowCols = ceiling . sqrt . genericLength $ features++      colorMap = getColorMap features+      allSelections =+        ( VL.FilterOp+        $ VL.FCompose+            ( F.foldl'+              (\acc x -> VL.And acc (VL.Selection $ "pick_" <> unFeature x))+              (VL.Expr "true")+              features+            )+        , [VL.MString "white"]+        )++      circleSpec = getCircleSpec labelMap range colorMap features+      windowSpecs = fmap (getWindowSpec colorMap) features+      allSpec = VL.hConcat+              $ [ maybe+                    circleSpec+                    (\ lm+                    -> VL.asSpec+                        [VL.columns 2, VL.vlConcat [circleSpec, legendSpec lm]]+                    )+                    labelMap+                , VL.asSpec [VL.columns numWindowCols, VL.vlConcat windowSpecs]+                ]++      p = VL.toVegaLite+            $ [ dataSet+              , VL.theme VL.defaultConfig []+              , allSpec+              ]++  TU.mktree . TU.fromText . T.pack $ unOutputDirectory outputDir'+  let fileName (Sample x) = T.unpack x <> "_projection.html"+      outputPath = unOutputDirectory outputDir' FP.</> fileName sample+  VL.toHtmlFile outputPath p

@@ -0,0 +1,248 @@+{- TooManyCells.Spatial.Relationships+Gregory W. Schwartz++Collects the functions pertaining to quantifying relationships between points.+-}++{-# LANGUAGE QuasiQuotes #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE OverloadedStrings #-}++module TooManyCells.Spatial.Relationships+    ( spatialRelationshipsR+    ) where++-- Remote+import BirchBeer.Types (LabelMap (..))+import Data.Bool (bool)+import Data.Maybe (isJust)+import Language.R as R+import Language.R.QQ (r)++-- Local+import TooManyCells.File.Types (OutputDirectory (..))+import TooManyCells.Matrix.Types (SingleCells (..), ProjectionMap (..), X (..), Y (..))+import TooManyCells.Matrix.Utility (sparseMatToSparseRMat)+import TooManyCells.Spatial.Utility (scToRPat)+import qualified TooManyCells.Spatial.Types as Too++spatialRelationshipsR :: OutputDirectory+                      -> Too.PCFCrossFlag+                      -> ProjectionMap+                      -> Maybe LabelMap+                      -> SingleCells+                      -> Maybe Too.TopDistances+                      -> [Too.Mark]+                      -> IO ()+spatialRelationshipsR (OutputDirectory outDir) pcfCrossFlag' pm lm sc td marks = R.runRegion $ do+  let isLabelMarks = bool 0 1 . isJust $ lm :: Double+      pcfCrossFlagDouble = bool 0 1 $ Too.unPCFCrossFlag pcfCrossFlag' :: Double+      (tVal, tQFlag) = case td of+                        Nothing -> (10, 1 :: Double)+                        Just (Too.TopQuantile x) -> (x :: Double, 1)+                        Just (Too.TopDistance x) -> (x :: Double, 0)++  pat <- scToRPat lm pm sc marks++  [r|+    suppressMessages(library(spatstat))+    suppressMessages(library(reshape2))+    suppressMessages(library(plyr))++    # Find index of first switch from above to below or equal value.+    findNegSwap = function(x, xs) {+      if(all(xs$obs > x)) {+        return(xs$r[length(xs$r)])+      } else if(all(xs$obs <= x)) {+        return(xs$r[1])+      } else {+        for (i in c(1:length(xs$r))) {+          if (i == length(xs$r)) {+            return(xs$r[i])+          } else if ((xs$obs[i] > x) && (xs$obs[i + 1] <= x)) {+            return(xs$r[i + 1])+          }+        }+      }+    }++    # Find index of first switch from below or equal value to above.+    findPosSwap = function(x, xs) {+      if (all(xs$obs <= x)) {+        return(xs$r[length(xs$r)])+      } else if (all(xs$obs > x)) {+        return(xs$r[1])+      } else {+        for (i in c(1:length(xs$r))) {+          if (i == length(xs$r)) {+            return(xs$r[i])+          } else if (xs$obs[i] <= x && xs$obs[i + 1] > x) {+            return(xs$r[i + 1])+          }+        }+      }+    }++    # Find the longest length of a continuously positive (or negative / equal to) stretch.+    findLongestLength = function(posFlag, x, xs) {+      if(posFlag) {+        lengthInfo = rle(xs$obs > x)+      } else {+        lengthInfo = rle(xs$obs <= x)+      }+      stretches = lengthInfo$lengths[which(lengthInfo$values)]+      if(length(stretches) == 0) {+        return(0)+      } else {+        return(xs$r[max(stretches)] - xs$r[1])+      }+    }++    # Find the first r value of the maximum value.+    findMaxPos = function(xs) {+      return(xs$r[min(which(xs$obs == max(xs$obs)))])+    }++    # Find the first r value of the maximum value with maximum value.+    findMaxPosWithVal = function(as) {+      # Do not use 0 values+      xs = as[-1,]+      r = findMaxPos(xs)+      val = max(xs$obs)++      return(val / r)+    }++    # Find the first r value of the minimum value.+    findMinPos = function(xs) {+      return(xs$r[min(which(xs$obs == min(xs$obs)))])+    }++    topDistances = function(xs) {+      if(tQFlag_hs == 1) {+        return(xs$obs[1:round(length(xs$obs) / tVal_hs)])+      } else {+        return(xs$obs[1:max(which(xs$r <= tVal_hs))])+      }+    }++    # Find statistics for a single curve+    getCurveStats = function (outFolder, df) {+      statsDf = data.frame(meanCorr = mean(df$obs))+      statsDf$maxCorr = max(df$obs)+      statsDf$minCorr = min(df$obs)+      statsDf$meanCorr = mean(df$obs)+      statsDf$topMaxCorr = max(topDistances(df))+      statsDf$topMeanCorr = mean(topDistances(df))+      statsDf$topMinCorr = min(topDistances(df))+      statsDf$negSwap = findNegSwap(1, df)+      statsDf$posSwap = findPosSwap(1, df)+      statsDf$longestPosLength = findLongestLength(TRUE, 1, df)+      statsDf$longestNegLength = findLongestLength(FALSE, 1, df)+      statsDf$maxPosWithVal = findMaxPosWithVal(df)+      statsDf$logMaxPosWithVal = log(findMaxPosWithVal(df))+      statsDf$maxPos = findMaxPos(df)+      statsDf$minPos = findMinPos(df)+      statsDf$label = basename(outFolder)+      statsDf$topVal = tVal_hs+      statsDf$topQuantile = tQFlag_hs++      return(statsDf)+    }++    crossCorrCurveStats = function (outFolder) {++      # Cross correlation function (several columns of marks).+      if(pcfCrossFlagDouble_hs == 1) {+        crossFn = envelope(pat_hs, pcfcross, funargs=marks)+        crossFn$obs = crossFn$iso+      } else {+        crossFn = markcrosscorr(pat_hs)+        for(i in c(1:length(crossFn$fns))) {+          crossFn$fns[[i]]$obs = crossFn$fns[[i]]$iso+        }+      }+      plot(crossFn)++      # Object+      outObject = crossFn+      outObject$label = basename(outFolder)+      saveRDS(outObject, file = file.path(outFolder,"crosscorr.rds"))++      # Statistics+      if(pcfCrossFlagDouble_hs == 1) {+        statsDf = getCurveStats(outFolder, crossFn)+      } else {+        varDf = melt(crossFn$which)+        varDf = varDf[order(varDf$value),]+        statsDf = cbind(varDf, do.call(rbind, lapply(crossFn$fns, function (x) getCurveStats(outFolder, x))))+      }++      # To get the sample size, different for quantitative vs. nominal marks+      markVars = marks(pat_hs)+      if (is.data.frame(markVars)) {+        statsDf$n = nrow(markVars)+      } else {+        statsDf$n = length(markVars[markVars != "Other"])+      }++      # Output statistics+      write.csv(statsDf, file.path(outFolder, "stats.csv"), quote = FALSE, row.names = FALSE)++      # Curve+      if(pcfCrossFlagDouble_hs == 1) {+        curveDf = data.frame(x = crossFn$r, y = crossFn$obs, label = basename(outFolder))+      } else {+        curveDf = adply(melt(crossFn$which), 1, function(xs) {+          data.frame(x = crossFn$fns[[xs$value]]$r, y = crossFn$fns[[xs$value]]$obs, label = basename(outFolder))+        })+      }+      write.csv(curveDf, file.path(outFolder, "curve.csv"), quote = FALSE, row.names = FALSE)++      return ()+    }++    mainRelationships = function(isLabelMarks, outFolder) {+      message(paste("Processing", outFolder))+      dir.create(outFolder, showWarnings = FALSE, recursive = TRUE)++      message("Plotting point process")+      pdf(file = file.path(outFolder, "basic_plot.pdf"))+      tryCatch(plot(pat_hs), error = function (e) message(paste("Could not plot basic_plot.pdf:", e)))+      dev.off()++      message("Plotting mark correlation function")+      # Mark correlation function.+      pdf(file = file.path(outFolder, "mark_correlation_function.pdf"))+      tryCatch(plot(markcorr(pat_hs)), error = function (e) message(paste("Could not plot mark_correlation_function.pdf:", e)))+      dev.off()++      # Mark variogram (lower value, more similar mark values at a distance).+      if (!isLabelMarks) {+        message("Plotting variogram")+        pdf(file = file.path(outFolder, "mark_variogram.pdf"))+        tryCatch(plot(markvario(pat_hs)), error = function (e) message(paste("Could not plot mark_variogram:", e)))+        dev.off()+      }++      # Envelope.+      message("Plotting envelope")+      pdf(file = file.path(outFolder, "envelope.pdf"))+      tryCatch(plot(envelope(pat_hs)), error = function (e) message(paste("Could not plot envelope", e)))+      dev.off()++      # Mark cross correlation+      message("Plotting mark cross correlation")+      pdf(file = file.path(outFolder, "cross_correlation_function.pdf"))+      tryCatch(crossCorrCurveStats(outFolder), error = function (e) message(paste("Could not compute mark cross correlation statistics:", e)))+      dev.off()++      return()++    }++    mainRelationships(isLabelMarks_hs, outDir_hs)++  |]++  pure ()

@@ -0,0 +1,597 @@+{- TooManyCells.Spatial.SummaryPlot+Gregory W. Schwartz++Collects functions pertaining to plotting the summary of statistics for all+samples.+-}++{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE PackageImports #-}++module TooManyCells.Spatial.SummaryPlot+    ( plotSummary+    ) where++import Control.Exception (SomeException (..), try, evaluate)+import Control.Monad (when, mfilter, guard)+import Data.Bool (bool)+import Data.Function (on)+import Data.List (sortBy, groupBy, foldl', sort, maximumBy, zipWith4)+import Data.Maybe (fromMaybe, mapMaybe, isNothing)+import Ploterific.Plot.Plot (labelColorScale)+import Ploterific.Plot.Types (ColorLabel (..))+import Safe (headMay)+import System.Environment (getArgs)+import System.IO (hPutStrLn, stderr)+import TextShow (showt)+import qualified BirchBeer.Types as Birch+import qualified Control.Foldl as Fold+import qualified Control.Lens as L+import qualified Data.ByteString.Lazy.Char8 as BL+import qualified Data.ByteString.Char8 as B+import qualified Data.Csv.Streaming as CSVStream+import qualified Data.Csv as CSV+import qualified Data.Foldable as F+import qualified Data.Map.Strict as Map+import qualified Data.Sequence as Seq+import qualified Data.Set as Set+import qualified Data.Text as T+import qualified Data.Text.Read as T+import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector as V+import qualified Graphics.Vega.VegaLite as VL+import qualified Graphics.Vega.VegaLite.Theme as VL+import qualified Statistics.Test.KruskalWallis as S+import qualified "statistics" Statistics.Quantile as S+import qualified "statistics" Statistics.Types as S+import qualified Filesystem.Path as FP+import qualified Turtle as TU++import qualified TooManyCells.File.Types as Too+import qualified TooManyCells.Spatial.Types as Too++-- | Get all pairwise comparisons.+pairwise :: (Eq a, Ord a) => [a] -> [[a]]+pairwise xs = Set.toList+            . Set.fromList+            . filter (\[x,y] -> x /= y)+            . fmap sort+            $ (\x y -> [x,y]) <$> xs <*> xs++-- | Safer Kruskal-Wallis test.+safeKWTest :: (VU.Unbox a, Ord a) => [VU.Vector a] -> IO Double+safeKWTest xs = do+  res <- try ( evaluate+             . maybe 1 (S.pValue . S.testSignificance)+             $ S.kruskalWallisTest xs+             ) :: IO (Either SomeException Double)+  res' <- case res of+            Left ex -> do+              hPutStrLn stderr $ "Caught exception, continuing with p = 1: " <> show ex+              return 1+            Right val -> return val+  return res'++-- | Add on the relevant p-values+pVal :: Birch.Feature+     -> [Map.Map T.Text T.Text]+     -> IO ([Map.Map T.Text T.Text], [Map.Map T.Text T.Text])+pVal (Birch.Feature feature) ms = do+  statRows' <- statRows+  statSummary' <- statSummary+  return (statRows', statSummary')+  where+    statRows = do+      tests' <- tests+      pairwiseTests' <- pairwiseTests+      return+        . mconcat+        . zipWith3+            (\ x xs ys+            -> fmap (\ m+                    -> foldl' (\acc (!a, !b) -> Map.insert a b acc) m (("kwPValue", x):xs)+                    )+            . concatMap Too.unLabelList+            . Too.unVarList+            $ ys+            )+            tests'+            (fmap Too.unVarList pairwiseTests')+        $ grouped+    statSummary = do+      tests' <- tests+      pairwiseTests' <- pairwiseTests+      return+        . mconcat+        . zipWith4+              (\ s x xs ys+              -> fmap (\ m+                      -> foldl'+                          (\acc (!a, !b) -> Map.insert a b acc)+                          m+                          (("highestLabels", s):("kwPValue", x):xs)+                      )+              . concatMap Too.unLabelList+              . Too.unVarList+              $ ys+              )+              highestLabels+              tests'+              (fmap Too.unVarList pairwiseTests')+          $ summarizedVarList+    highestLabels :: [T.Text]+    highestLabels = fmap ( T.intercalate "/"+                         . fmap showt+                         . sortBy ((flip compare) `on` fst)  -- Descending+                         . concatMap+                             ( fmap+                                 (\ !x+                                 -> ( maybe+                                       (error "No median found")+                                       (either error fst . T.double)+                                    $ Map.lookup "median" x+                                    , Map.findWithDefault+                                       (error "No label found")+                                       "statLabel"+                                       x+                                    )+                                 )+                                 . Too.unLabelList+                             )+                         . Too.unVarList+                         )+                  $ summarizedVarList+    summarizedVarList :: [Too.VarList (Too.LabelList (Map.Map T.Text T.Text))]+    summarizedVarList = Too.VarList+                      . fmap ( Too.LabelList+                             . (:[])+                             . getStatSummaries+                             . Too.unLabelList+                             )+                      . Too.unVarList+                    <$> grouped+    getStatSummaries !ms = Map.fromList+                         $ ("median", T.pack . show $ getMedianFromMaps ms)+                         : ("feature", feature)+                         : foldl'+                            (\ acc x+                            -> ( x+                               , fromMaybe "" $ Map.lookup x =<< headMay ms+                               )+                               : acc+                            )+                            []+                            ["statLabel", "Var1", "Var2"]+    getMedianFromMaps = S.median S.s . VU.fromList . fmap getFeature+    kruskalOrBust = fmap (T.pack . show) . safeKWTest+    tests :: IO [T.Text]+    tests = mapM ( kruskalOrBust+                 . fmap (VU.fromList . fmap getFeature . Too.unLabelList)+                 . Too.unVarList+                 )+          $ grouped+    pairwiseTests :: IO [Too.VarList (T.Text, T.Text)]+    pairwiseTests = mapM ( fmap Too.VarList+                         . mapM+                            ( sequence+                            . L.over+                                L._2+                                ( kruskalOrBust+                                . fmap+                                    ( VU.fromList+                                    . fmap getFeature+                                    . Too.unLabelList+                                    )+                                )+                                . Too.unCompare+                            )+                         . setupComparisonVar+                         )+                   $ grouped+    setupComparisonVar :: Too.VarList (Too.LabelList (Map.Map T.Text T.Text))+                       -> [Too.Compare (Too.LabelList (Map.Map T.Text T.Text))]+    setupComparisonVar = fmap ( Too.Compare+                              . L.over L._1 (T.intercalate "/")+                              . unzip+                              )+                       . (\x -> pairwise x :: [[(T.Text, Too.LabelList (Map.Map T.Text T.Text))]])+                       . fmap ( L.over L._2 Too.LabelList+                              . L.over L._1 (fromMaybe "" . headMay)+                              . unzip+                              . fmap (\ !x -> (fromMaybe "" $ compLabel x, x))+                              . Too.unLabelList+                              )+                       . Too.unVarList+    getFeature = maybe 0 (either error fst . T.double) . Map.lookup feature+    grouped :: [Too.VarList (Too.LabelList (Map.Map T.Text T.Text))]+    grouped = fmap ( Too.VarList+                   . fmap Too.LabelList+                   . groupBy ((==) `on` compLabel)+                   . sortBy (compare `on` compLabel)+                   )+            . groupBy ((==) `on` compVar)+            . sortBy (compare `on` compVar)+            $ ms+    compLabel = Map.lookup "statLabel"+    compVar x = (Map.lookup "Var1" x, Map.lookup "Var2" x)++-- -- | Add marksOther-markOther label if present to every group to compare with+-- -- all others.+-- addOthersOthersToEachGroup+--   :: [Too.VarList (Too.LabelList (Map.Map T.Text T.Text))]+--   -> [Too.VarList (Too.LabelList (Map.Map T.Text T.Text))]+-- addOthersOthersToEachGroup varLists = fmap addToVarList noOthersVarList+--   where+--     noOthersVarList =+--       filter (not . (\x -> isOthersList (&&) x && isOthersList (||) x)) varLists+--     addToVarList xs = xs <> mconcat others+--     others :: [Too.VarList (Too.LabelList (Map.Map T.Text T.Text))]+--     others = filter (isOthersList (&&)) varLists+--     isOthersList :: (Bool -> Bool -> Bool)+--                  -> Too.VarList (Too.LabelList (Map.Map T.Text T.Text))+--                  -> Bool+--     isOthersList f xs = fromMaybe False $ do+--       var1 <- getVarFromVarList "Var1" xs+--       var2 <- getVarFromVarList "Var2" xs+--       return $ f (var1 == "marksOther") (var2 == "marksOther")++-- -- | Add markOther label if present to every group to compare with+-- -- all others.+-- addOthersToGroup :: [Too.VarList (Too.LabelList (Map.Map T.Text T.Text))]+--                  -> [Too.VarList (Too.LabelList (Map.Map T.Text T.Text))]+-- addOthersToGroup varLists = fmap addToVarList varLists+--   where+--     addToVarList xs = fromMaybe xs $ do+--       var1 <- getVarFromVarList "Var1" xs+--       others <- Map.lookup var1 otherMap+--       return $ xs <> others+--     otherMap = otherComparisonsMap varLists++-- -- | Get marksOther-marksOther label if exists.+-- otherComparisonsMap+--   :: [Too.VarList (Too.LabelList (Map.Map T.Text T.Text))]+--   -> Map.Map T.Text (Too.VarList (Too.LabelList (Map.Map T.Text T.Text)))+-- otherComparisonsMap = Map.unions . mapMaybe getVarListMap+--   where+--     getVarListMap xs = do+--       var1 <- getVarFromVarList "Var1" xs+--       var2 <- getVarFromVarList "Var2" xs+--       guard (var2 == "marksOther")+--       return $ Map.singleton var1 xs++-- | Return true if this is an other comparison row+isOther :: Map.Map T.Text T.Text -> Bool+isOther m = fromMaybe False $ do+      var1 <- Map.lookup "Var1" m+      var2 <- Map.lookup "Var2" m+      return (var1 == "marksOther" || var2 == "marksOther")++-- | Get a Var column value from a VarList+getVarFromVarList :: T.Text+                  -> Too.VarList (Too.LabelList (Map.Map T.Text T.Text))+                  -> Maybe T.Text+getVarFromVarList var xs = (headMay . Too.unVarList $ xs)+                       >>= headMay . Too.unLabelList+                       >>= Map.lookup var+-- | Parse a file into CSV.+parseFile :: Maybe Too.StateLabelMap -> TU.FilePath -> IO [Map.Map T.Text T.Text]+parseFile slm file = do+  contents <- BL.readFile . T.unpack . TU.format TU.fp $ file+  -- Force handles to close+  let !rows = either error ( fmap (insertSampleLabel slm file)+                           . insertOldVar2+                           . F.toList+                           . snd+                           )+            . CSVStream.decodeByName+            $ contents+  return rows++-- | Insert samples and labels in a row. Labels are for plotting niceness.+insertSampleLabel :: Maybe Too.StateLabelMap+                  -> TU.FilePath+                  -> Map.Map T.Text T.Text+                  -> Map.Map T.Text T.Text+insertSampleLabel slm file m =+  Map.insert "label" (fromMaybe "No label" label)+    . Map.insert "statLabel" statLabel+    . Map.insert "sample" sample+    $ m+  where+    statLabel = maybe (fromMaybe "No label" alternateLabel) id label+    label+      | isOldVar2OtherOnly m = (\x y -> x <> "-" <> y)+                           <$> state+                           <*> (Just "marksOther")+      | otherwise = state+      where+        state = fmap Birch.unLabel+              . (=<<) (Map.lookup (Birch.Id sample) . Too.unStateLabelMap)+              $ slm+    alternateLabel = (\x y -> x <> "-" <> y)+                 <$> Map.lookup "Var1" m+                 <*> Map.lookup "Var2" m+    sample = Birch.unSample $ getSample file++-- | Insert the original Var2 for reference in the other comparison.+insertOldVar2 :: [Map.Map T.Text T.Text] -> [Map.Map T.Text T.Text]+insertOldVar2 ms = mapMaybe insertOldVar2 ms+  where+    insertOldVar2 m = do+      var2 <- Map.lookup "Var2" m+      relevantVar <- getRelevantVar m+      return+        . bool id (Map.insert "Var2" relevantVar) (isVar2OtherOnly m)+        . Map.insert "oldVar2" var2+        $ m+    getRelevantVar m = do+      var1 <- Map.lookup "Var1" m+      var2 <- Map.lookup "Var2" m+      if (var1 == var2)+        then return var1+        else do+          -- If Nothing, then it is Var1 vs. Var1 Other, so use Var1.+          oldVar2 <- Just . fromMaybe var1 . headMay . filter (/= var1) $ vars+          return oldVar2+    vars = filter (/= "marksOther")+         . mapMaybe (Map.lookup "Var2")+         $ ms++-- | Check if Var2 is marksOther.+isVar2OtherOnly :: Map.Map T.Text T.Text -> Bool+isVar2OtherOnly m = fromMaybe False $ do+  var1 <- Map.lookup "Var1" m+  var2 <- Map.lookup "Var2" m+  return (var1 /= "marksOther" && var2 == "marksOther")++-- | Check if oldVar2 is marksOther.+isOldVar2OtherOnly :: Map.Map T.Text T.Text -> Bool+isOldVar2OtherOnly m = fromMaybe False $ do+  var1 <- Map.lookup "Var1" m+  var2 <- Map.lookup "oldVar2" m+  return (var1 /= "marksOther" && var2 == "marksOther")++-- | Check if Var1 is marksOther.+isVar1Other :: Map.Map T.Text T.Text -> Bool+isVar1Other m = fromMaybe False $ do+  var1 <- Map.lookup "Var1" m+  return (var1 == "marksOther")++-- | Get sample from file path.+getSample :: TU.FilePath -> Birch.Sample+getSample = Birch.Sample . TU.format TU.fp . TU.dirname . p . p . p+  where+    p = TU.parent++-- | Remove columns with "marksOther".+removeOthers :: [Map.Map T.Text T.Text] -> [Map.Map T.Text T.Text]+removeOthers =+  fmap (Map.filterWithKey (\k _ -> not . T.isInfixOf "marksOther" $ k))++-- | Remove rows with "marksOther".+removeOthersRows :: [Map.Map T.Text T.Text] -> [Map.Map T.Text T.Text]+removeOthersRows = filter (\ x -> (not . isOther $ x)+                               && (Map.lookup "oldVar2" x /= Just "marksOther")+                          )++-- | Make sure all columns are represented in a row+fillColumns :: [Map.Map T.Text T.Text] -> [Map.Map T.Text T.Text]+fillColumns xs = fmap addDummy xs+  where+    addDummy = flip Map.union dummyColMap  -- Prefer original value, not dummy.+    dummyColMap = Map.fromList+                . fmap (\ !x -> (x, "NA"))+                . Set.toList+                . Set.fromList+                . concatMap Map.keys+                $ xs++plotSummary :: Too.OutputDirectory+            -> Too.IncludeOthersFlag+            -> Maybe Too.StateLabelMap+            -> Birch.Feature+            -> IO ()+plotSummary (Too.OutputDirectory inDir) iof slm feature = do+  let inDir' = TU.fromText . T.pack $ inDir+  statsFiles <- TU.reduce Fold.list . TU.find (TU.suffix "stats.csv") $ inDir'+  statsParsed <- mapM (parseFile slm) statsFiles++  when (null statsParsed) $ print "No stats.csv files found, first generate spatial results"++  let outDir = inDir' FP.</> "summary"++  TU.mktree outDir++  let outputBasename = "summary_" <> Birch.unFeature feature+      outputPlot = outDir FP.</> TU.fromText (outputBasename <> ".html")+      outputStats = outDir FP.</> TU.fromText (outputBasename <> ".csv")+      rows' = bool+                removeOthersRows+                (filter (not . isVar1Other))+                (Too.unIncludeOthersFlag iof)+            . mconcat+            $ statsParsed+      noLabelFlag = isNothing slm++  (rows, stats) <- pVal feature rows'++  let pValComparisons = fmap (T.intercalate "/" . fmap (fromMaybe ""))+                      . pairwise+                      . Set.toList+                      . Set.fromList+                      . fmap (Map.lookup "label")+                      $ rows'+      colorLabels = fmap ColorLabel+                  . Set.toAscList+                  . Set.fromList+                  $ fmap (fromMaybe "" . Map.lookup "label") rows+      feature' = Birch.unFeature feature+      dataSet = VL.dataFromColumns []+              . VL.dataColumn+                  "var1"+                  (VL.Strings . fmap (Map.findWithDefault "" "Var1") $ rows)+              . VL.dataColumn+                  "var2"+                  (VL.Strings . fmap (Map.findWithDefault "" "Var2") $ rows)+              . VL.dataColumn+                  "label"+                  (VL.Strings . fmap (Map.findWithDefault "" "label") $ rows)+              . VL.dataColumn+                  "sample"+                  (VL.Strings . fmap (Map.findWithDefault "" "sample") $ rows)+              . VL.dataColumn+                  "kwPValue"+                  ( VL.Numbers+                  . fmap+                      ( maybe 0 (either error fst . T.double)+                      . Map.lookup "kwPValue"+                      )+                  $ rows+                  )+              . (\ y+                -> foldl'+                    (\ acc x+                    -> VL.dataColumn+                        x+                        ( VL.Numbers+                        . fmap ( maybe 0 (either error fst . T.double)+                               . Map.lookup x+                               )+                        $ rows+                        ) acc+                    )+                    y+                    pValComparisons+                )+              $ VL.dataColumn+                  feature'+                  ( VL.Numbers+                  . fmap+                      ( maybe 0 (either error fst . T.double)+                      . Map.lookup feature'+                      )+                  $ rows+                  )+                  []+      transDens = transMax+                . VL.density feature' [ VL.DnGroupBy ["var1", "var2", "label"]]+      transHref = VL.calculateAs ("'../' + datum.sample + '/projections/' + datum.sample + '_projection.html'") "url"+      transMax = VL.joinAggregate+                  [ VL.opAs VL.Max feature' "max" ]+                  [ VL.WGroupBy ["var1", "var2"] ]+               . VL.calculateAs ("datum." <> feature' <> "/datum.max") "maxTrans"+               . transHref+      picked = "picked"+      sel = VL.selection+          . VL.select+              picked+              VL.Interval+              [ VL.Encodings [ VL.ChX, VL.ChY ], VL.BindScales ]+      axisSize = VL.Axis [ VL.TickWidth 0.756, VL.DomainWidth 0.756 ]+      densitySpec = [ VL.mark VL.Area [ VL.MOpacity 0.5]+                    , sel []+                    , VL.encoding+                    . VL.position+                        VL.Y+                        [ VL.PName "density"+                        , VL.PmType VL.Quantitative+                        ]+                    . VL.position+                        VL.X+                        [ VL.PName "value"+                        , VL.PmType VL.Quantitative+                        ]+                    . VL.color+                        [ VL.MName "label"+                        , VL.MmType VL.Nominal+                        , labelColorScale colorLabels+                        ]+                    $ []+                    , VL.transform $ transDens []+                    ]+      tickSpec = [ VL.mark VL.Tick []+                 , VL.transform $ transMax []+                 , VL.encoding+                 . VL.position+                    VL.X+                    [VL.PName feature', VL.PmType VL.Quantitative]+                 . VL.tooltips ( fmap+                                  (\x ->  [VL.TName x, VL.TmType VL.Nominal])+                                  ( [ feature'+                                    , "sample"+                                    , "label"+                                    , "kwPValue"+                                    ]+                                 <> pValComparisons+                                  )+                               )+                 . VL.color [ VL.MName "label"+                            , VL.MmType VL.Nominal+                            , labelColorScale colorLabels+                            ]+                 . VL.hyperlink [ VL.HName "url", VL.HmType VL.Nominal ]+                 $ []+                 ]+      boxSpec = [ VL.mark+                    VL.Boxplot+                    [ VL.MNoOutliers+                    , VL.MOpacity 0.5+                    , VL.MMedian [ VL.MColor "black" ]+                    ]+                , VL.transform $ transMax []+                , VL.encoding+                . VL.position+                    VL.X+                    [VL.PName feature', VL.PmType VL.Quantitative]+                . VL.position VL.Y [ VL.PNumber 30 ]+                . VL.color+                    [ VL.MName "label"+                    , VL.MmType VL.Nominal+                    , labelColorScale colorLabels+                    ]+                $ []+                ]+      res = VL.resolve+          . VL.resolution ( VL.RAxis+                              [ (VL.ChX, VL.Independent)+                              , (VL.ChY, VL.Independent)+                              ]+                          )+      p = VL.toVegaLite [ dataSet+                        , VL.specification+                        $ VL.asSpec+                            [ VL.layer+                                (fmap VL.asSpec [densitySpec, tickSpec, boxSpec])+                            ]+                        , VL.facet [ VL.RowBy+                                      [VL.FName "var1", VL.FmType VL.Nominal]+                                   , VL.ColumnBy+                                      [VL.FName "var2", VL.FmType VL.Nominal]+                                   ]+                        , res []+                        , VL.theme+                            ( VL.defaultConfig+                                { VL.configHeight = Just 170+                                , VL.configWidth = Just 200+                                , VL.configFontSize = Just 8+                                , VL.configTitleFontSize = Just 9.33333015+                                }+                            )+                            (VL.configuration axisSize [])+                        ]++  -- Output plot+  VL.toHtmlFile (T.unpack . TU.format TU.fp $ outputPlot) p++  -- Output stats+  let header = V.fromList+             . fmap (B.pack . T.unpack)+             . Set.toList+             . Set.fromList+             . concatMap Map.keys+             $ stats+  BL.writeFile (T.unpack . TU.format TU.fp $ outputStats)+    . CSV.encodeByName header+    . fillColumns+    $ stats

@@ -0,0 +1,42 @@+{- TooManyCells.Spatial.Types+Gregory W. Schwartz++Collects the spatial types used in the program.+-}++{-# LANGUAGE StrictData #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}++module TooManyCells.Spatial.Types where++-- Remote+import Control.Applicative (Alternative (..))+import Control.Monad (Monad (..), MonadPlus (..))+import Data.Colour.Palette.BrewerSet (Kolor)+import qualified BirchBeer.Types as Birch+import qualified Data.Map.Strict as Map+import qualified Data.Text as T++-- Local++newtype ColorMap = ColorMap { unColorMap :: Map.Map Birch.Feature Kolor }+newtype StartCol = StartCol { unStartCol :: T.Text }+newtype EndCol = EndCol { unEndCol :: T.Text }+newtype AnnoSpatMarkerFile = AnnoSpatMarkerFile { unAnnoSpatMarkerFile :: String }+newtype AnnoSpatCommand = AnnoSpatCommand { unAnnoSpatCommand :: String }+newtype PCFCrossFlag = PCFCrossFlag { unPCFCrossFlag :: Bool }+newtype SkipFinishedFlag = SkipFinishedFlag { unSkipFinishedFlag :: Bool }+newtype IncludeOthersFlag = IncludeOthersFlag { unIncludeOthersFlag :: Bool }+newtype StateLabelsFile = StateLabelsFile { unStateLabelsFile :: String }+newtype StateLabelMap = StateLabelMap { unStateLabelMap :: Map.Map Birch.Id Birch.Label }+newtype LabelList a = LabelList { unLabelList :: [a] } deriving (Eq, Ord, Show, Functor, Monoid, Semigroup, Applicative, Alternative, Monad, MonadPlus)+newtype Compare a = Compare { unCompare :: (T.Text, [a]) } deriving (Show)+newtype VarList a = VarList { unVarList :: [a] } deriving (Show, Functor, Semigroup, Monoid, Applicative, Alternative, Monad, MonadPlus)++data Range = Range { minX :: Double+                   , maxX :: Double+                   , minY :: Double+                   , maxY :: Double+                   }+data Mark = MarkFeature Birch.Feature | MarkLabel Birch.Label deriving (Show, Read, Eq, Ord)+data TopDistances = TopQuantile Double | TopDistance Double deriving (Show, Read)

@@ -0,0 +1,146 @@+{- TooManyCells.Spatial.Utility+Gregory W. Schwartz++Collects helper functions in the program for spatial analysis.+-}++{-# LANGUAGE QuasiQuotes #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE TupleSections #-}++module TooManyCells.Spatial.Utility+    ( scToVLData+    , subsampleProjectionMap+    , scToRPat+    , markToText+    ) where++-- Remote+import BirchBeer.Types+import Control.Monad (liftM2, liftM3)+import Data.Bool (bool)+import Data.List (zipWith3, foldl')+import Data.Maybe (catMaybes, fromMaybe, isJust)+import Language.R as R+import Language.R.QQ (r)+import TooManyCells.Matrix.Types+import qualified Control.Lens as L+import qualified Data.Map.Strict as Map+import qualified Data.Sparse.Common as S+import qualified Data.Set as Set+import qualified Data.Text as T+import qualified Data.Vector as V+import qualified Graphics.Vega.VegaLite as VL++-- Local+import TooManyCells.Matrix.Utility (sparseMatToSparseRMat, subsetCellsSc)+import TooManyCells.Spatial.Types (Mark (..))++-- | SingleCells to Vega Data+scToVLData :: Maybe LabelMap -> ProjectionMap -> SingleCells -> VL.Data+scToVLData lm pm sc =+  VL.dataFromRows []+    . foldl' (\acc x -> VL.dataRow x acc) []+    . catMaybes  -- remove missing projections+    . joinAll  -- required+        projectionOrdered+        (fmap (Just . ("item",) . VL.Str) . V.toList . getRowNames $ sc)+    . fmap Just+    . maybe id (\x -> zipWith (:) (labelOrdered x)) lm  -- optional+    . fmap (zip (V.toList . getColNames $ sc) . fmap VL.Number . S.toDenseListSV)+    . S.toRowsL+    . getMatrix+    $ sc+  where+    joinAll :: [Maybe [(T.Text, VL.DataValue)]]+            -> [Maybe (T.Text, VL.DataValue)]+            -> [Maybe [(T.Text, VL.DataValue)]]+            -> [Maybe [(T.Text, VL.DataValue)]]+    joinAll = zipWith3 (\xs y z -> liftM3 (\as b c -> (as <> (b:c))) xs y z)+    labelOrdered :: LabelMap -> [(T.Text, VL.DataValue)]+    labelOrdered (LabelMap lm') = V.toList+                                . fmap ( ("label",)+                                       . VL.Str+                                       . unLabel+                                       . fromMaybe (Label "NA")+                                       . flip Map.lookup lm'+                                       . Id+                                       )+                                . getRowNames+                                $ sc+    projectionOrdered = fmap getProjectionInfo . V.toList . L.view rowNames $ sc+    getProjectionInfo :: Cell -> Maybe [(T.Text, VL.DataValue)]+    getProjectionInfo item = Map.lookup item (unProjectionMap pm)+                         >>= pure . projectionToDataValueInfo . snd+    projectionToDataValueInfo :: (X, Y) -> [(T.Text, VL.DataValue)]+    projectionToDataValueInfo (X x, Y y) = [("x", VL.Number x), ("y", VL.Number y)]++-- | Subsample a projection map.+subsampleProjectionMap :: Maybe Sample -> ProjectionMap -> ProjectionMap+subsampleProjectionMap sample = ProjectionMap+                              . Map.filter ((== sample) . fst)+                              . unProjectionMap++-- | SingleCells to easy input for R patterns+scToRPat :: Maybe LabelMap+         -> ProjectionMap+         -> SingleCells+         -> [Mark]+         -> R.R s (R.SomeSEXP s)+scToRPat lm pm sc marks = do+  [r| suppressMessages(library(spatstat)) |]++  let cells  = V.toList . L.view rowNames $ sc+      ps     = fmap (\x -> Map.lookup x $ unProjectionMap pm) cells+      getLabel (MarkFeature _) = error "Expected MarkLabel, not MarkFeature in scToRPat."+      getLabel (MarkLabel l) = l+      labels = Set.fromList . fmap getLabel $ marks+      labelOrOther l =+        bool "Other" (T.unpack . unLabel $ l) . Set.member l $ labels+      ls     =+        fmap+          (\ lm'+          -> fmap+              (\ (Cell x)+              -> labelOrOther+               . Map.findWithDefault (Label "NA") (Id x)+               $ unLabelMap lm'+              )+            . fmap fst+            . catMaybes+            . zipWith (liftM2 (,)) (fmap Just cells)+            $ ps+          )+          lm+      subsampledSc = flip subsetCellsSc sc+                   . catMaybes+                   . zipWith (liftM2 const) (fmap Just cells)+                   $ ps+      isLabel = isJust lm+      xs = fmap (unX . fst) . fmap snd . catMaybes $ ps+      ys = fmap (unY . snd) . fmap snd . catMaybes $ ps++  mat <- fmap unRMatObsRow $ sparseMatToSparseRMat subsampledSc+  pat <- [r|+          ppp(xs_hs, ys_hs, c(min(xs_hs), max(xs_hs)), c(min(ys_hs),max(ys_hs)))+         |]++  case ls of+    Nothing -> do+      let getFeature (MarkFeature f) = T.unpack . unFeature $ f+          getFeature (MarkLabel f) = error "Expected MarkFeature, not MarkLabel in scToRPat."+          features = fmap getFeature marks+      marksDf <- [r| as.data.frame(as.matrix(mat_hs[, features_hs])) |]+      [r| p = pat_hs+          marks(p) = marksDf_hs+          p+      |]+    (Just l) -> [r| p = pat_hs+                    marks(p) = as.factor(l_hs)+                    p+                |]++-- | Get the mark to the text.+markToText (MarkFeature (Feature x)) = x+markToText (MarkLabel (Label x)) = x

@@ -1,13 +1,13 @@ name:                too-many-cells-version:             2.1.1.0+version:             3.0.1.0 synopsis:            Cluster single cells and analyze cell clade relationships. description:         Different methods to cluster and analyze single cell data with diversity indices and differential expression. homepage:            http://github.com/GregorySchwartz/too-many-cells#readme license:             GPL-3 license-file:        LICENSE author:              Gregory W. Schwartz-maintainer:          gsch@pennmedicine.upenn.edu-copyright:           2020 Gregory W. Schwartz+maintainer:          gregory.schwartz@uhnresearch.ca+copyright:           2022 Gregory W. Schwartz category:            Bioinformatics build-type:          Simple -- extra-source-files:@@ -55,7 +55,14 @@                      , TooManyCells.Program.Options                      , TooManyCells.Program.Paths                      , TooManyCells.Program.Peaks+                     , TooManyCells.Program.Spatial                      , TooManyCells.Program.Utility+                     , TooManyCells.Spatial.AnnoSpat+                     , TooManyCells.Spatial.ProjectionPlot+                     , TooManyCells.Spatial.Relationships+                     , TooManyCells.Spatial.SummaryPlot+                     , TooManyCells.Spatial.Types+                     , TooManyCells.Spatial.Utility   build-depends:       base >= 4.7 && < 5                      , IntervalMap                      , SVGFonts@@ -86,6 +93,8 @@                      , hierarchical-spectral-clustering                      , hmatrix                      , hmatrix-svdlibc+                     , hvega+                     , hvega-theme                      , inline-r                      , lens                      , managed@@ -93,10 +102,11 @@                      , modularity                      , mtl                      , mwc-random-                     , optparse-generic+                     , optparse-applicative                      , palette                      , parallel                      , plots+                     , ploterific                      , process                      , resourcet                      , safe@@ -123,16 +133,18 @@                      , vector                      , vector-algorithms                      , zlib-  ghc-options:         -O2+  -- No specConstr due to memory limitations+  ghc-options:         -O2 -fno-spec-constr   default-language:    Haskell2010  executable too-many-cells   hs-source-dirs:      app   main-is:             Main.hs-  ghc-options:         -threaded -rtsopts -O2+  -- No specConstr due to memory limitations+  ghc-options:         -threaded -rtsopts -O2 -fno-spec-constr   build-depends:       base                      , too-many-cells-                     , optparse-generic+                     , optparse-applicative   default-language:    Haskell2010  source-repository head