too-many-cells 0.2.2.2 → 2.1.0.1
raw patch · 39 files changed
+4321/−1527 lines, 39 filesdep +IntervalMapdep +asyncdep +async-pooldep ~SVGFontsdep ~aesondep ~base
Dependencies added: IntervalMap, async, async-pool, attoparsec, hashable, hmatrix-svdlibc, mwc-random, process, resourcet, stm, streaming-commons, system-filepath, turtle, unordered-containers
Dependency ranges changed: SVGFonts, aeson, base, birch-beer, bytestring, cassava, colour, containers, deepseq, diagrams, diagrams-cairo, diagrams-graphviz, diagrams-lib, differential, directory, diversity, fgl, filepath, find-clumpiness, foldl, graphviz, hierarchical-clustering, hierarchical-spectral-clustering, hmatrix, inline-r, lens, managed, matrix-market-attoparsec, modularity, mtl, optparse-generic, palette, parallel, plots, safe, scientific, sparse-linear-algebra, split, statistics, streaming, streaming-bytestring, streaming-cassava, streaming-utils, streaming-with, temporary, terminal-progress-bar, text, text-show, transformers, vector, vector-algorithms, zlib
Files
- app/Main.hs +34/−1067
- src/TooManyCells/Classify/Classify.hs +86/−0
- src/TooManyCells/Classify/Types.hs +18/−0
- src/TooManyCells/Differential/Differential.hs +212/−73
- src/TooManyCells/Differential/Types.hs +7/−0
- src/TooManyCells/Diversity/Plot.hs +3/−0
- src/TooManyCells/File/Types.hs +8/−1
- src/TooManyCells/MakeTree/Adjacency.hs +3/−3
- src/TooManyCells/MakeTree/Cluster.hs +8/−51
- src/TooManyCells/MakeTree/Load.hs +1/−1
- src/TooManyCells/MakeTree/Plot.hs +4/−41
- src/TooManyCells/MakeTree/Print.hs +61/−4
- src/TooManyCells/MakeTree/Types.hs +7/−5
- src/TooManyCells/MakeTree/Utility.hs +51/−0
- src/TooManyCells/Matrix/AtacSeq.hs +99/−0
- src/TooManyCells/Matrix/Load.hs +188/−36
- src/TooManyCells/Matrix/Preprocess.hs +354/−58
- src/TooManyCells/Matrix/Types.hs +282/−36
- src/TooManyCells/Matrix/Utility.hs +134/−8
- src/TooManyCells/Motifs/FindMotif.hs +155/−0
- src/TooManyCells/Motifs/Types.hs +26/−0
- src/TooManyCells/Paths/Distance.hs +25/−8
- src/TooManyCells/Paths/Plot.hs +20/−5
- src/TooManyCells/Paths/Types.hs +1/−0
- src/TooManyCells/Peaks/ClusterPeaks.hs +448/−0
- src/TooManyCells/Peaks/Types.hs +29/−0
- src/TooManyCells/Program/Classify.hs +64/−0
- src/TooManyCells/Program/Differential.hs +201/−0
- src/TooManyCells/Program/Diversity.hs +129/−0
- src/TooManyCells/Program/Interactive.hs +88/−0
- src/TooManyCells/Program/LoadMatrix.hs +267/−0
- src/TooManyCells/Program/MakeTree.hs +501/−0
- src/TooManyCells/Program/MatrixOutput.hs +31/−0
- src/TooManyCells/Program/Motifs.hs +64/−0
- src/TooManyCells/Program/Options.hs +305/−0
- src/TooManyCells/Program/Paths.hs +102/−0
- src/TooManyCells/Program/Peaks.hs +99/−0
- src/TooManyCells/Program/Utility.hs +70/−0
- too-many-cells.cabal +136/−130
app/Main.hs view
@@ -4,1070 +4,37 @@ Clusters single cell data. -} -{-# LANGUAGE QuasiQuotes #-}-{-# LANGUAGE BangPatterns #-}-{-# LANGUAGE DataKinds #-}-{-# LANGUAGE DeriveGeneric #-}-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE PackageImports #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE TupleSections #-}--module Main where---- Remote-import Control.Concurrent-import BirchBeer.ColorMap-import BirchBeer.Interactive-import BirchBeer.Load-import BirchBeer.MainDiagram-import BirchBeer.Plot-import BirchBeer.Types-import BirchBeer.Utility-import Control.Monad (when, unless, join)-import Control.Monad.Trans (liftIO)-import Control.Monad.Trans.Maybe (MaybeT (..))-import Data.Bool (bool)-import Data.Colour.SRGB (sRGB24read)-import Data.Matrix.MatrixMarket (readMatrix, writeMatrix)-import Data.Maybe (fromMaybe, isJust, isNothing)-import Data.Monoid ((<>))-import Data.Tree (Tree (..))-import Language.R as R-import Language.R.QQ (r)-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)-import qualified "find-clumpiness" Types as Clump-import qualified Control.Lens as L-import qualified Data.Aeson as A-import qualified Data.ByteString.Lazy.Char8 as B-import qualified Data.Colour.Palette.BrewerSet as D-import qualified Data.Colour.Palette.Harmony as D-import qualified Data.Csv as CSV-import qualified Data.GraphViz as G-import qualified Data.Set as Set-import qualified Data.Text as T-import qualified Data.Text.Lazy.IO as T-import qualified Data.Vector as V-import qualified Diagrams.Backend.Cairo as D-import qualified Diagrams.Prelude as D-import qualified H.Prelude as H-import qualified Plots as D-import qualified System.Directory as FP-import qualified System.FilePath as FP-import qualified System.ProgressBar as Progress---- Local-import TooManyCells.Differential.Differential-import TooManyCells.Differential.Types-import TooManyCells.Diversity.Diversity-import TooManyCells.Diversity.Load-import TooManyCells.Diversity.Plot-import TooManyCells.Diversity.Types-import TooManyCells.File.Types-import TooManyCells.MakeTree.Clumpiness-import TooManyCells.MakeTree.Cluster-import TooManyCells.MakeTree.Load-import TooManyCells.MakeTree.Plot-import TooManyCells.MakeTree.Print-import TooManyCells.MakeTree.Types-import TooManyCells.Matrix.Load-import TooManyCells.Matrix.Preprocess-import TooManyCells.Matrix.Types-import TooManyCells.Matrix.Utility-import TooManyCells.Paths.Distance-import TooManyCells.Paths.Plot-import TooManyCells.Paths.Types---- | 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 gene row names and cell column names. If given as a list (--matrixPath input1 --matrixPath input2 etc.) then will join all matrices together. Assumes the same number and order of genes in each matrix, so only cells are added."- , projectionFile :: Maybe String <?> "([Nothing] | FILE) The input file containing positions of each cell for plotting. Format is \"barcode,x,y\" and matches column order in the matrix file. 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."- , delimiter :: Maybe Char <?> "([,] | CHAR) The delimiter for the csv file 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 :: Maybe String <?> "([TfIdfNorm] | UQNorm | MedNorm | TotalMedNorm | BothNorm | LogCPMNorm | NoneNorm) Type of normalization before clustering. 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. TotalMedNorm normalized first each observation by total count then by median of non-zero counts across features. BothNorm uses both normalizations (first TotalMedNorm for all analysis then additionally TfIdfNorm during clustering). LogCPMNorm normalizes by log2(CPM + 1). 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). Cannot use TfIdfNorm for any other process as NoneNorm will become the default."- , 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. --min-proportion distribution is log2 transformed."- , 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. See --matrix-output-transpose."- , matrixOutputTranspose :: 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. Differs from --matrix-output in that features are columns.columns."- , 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, DOUBLE)]), where each item is colored by the binary high / low expression of FEATURE based on DOUBLE 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 | PieNone | 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. PieNone only draws items, no pie rings or charts. (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) 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."- , drawFont :: Maybe String <?> "([Arial] | FONT) Specify the font to use for the labels when plotting."- , 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."- , noFilter :: Bool <?> "Whether to bypass filtering genes and cells by low counts."- , shiftPositive :: Bool <?> "Shift features to positive values. Positive values are shifted to allow modularity to work correctly."- , filterThresholds :: Maybe String <?> "([(250, 1)] | (DOUBLE, DOUBLE)) The minimum filter thresholds for (MINCELL, MINFEATURE) when filtering cells and features by low read counts. See also --no-filter."- , 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."- , 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 gene row names and cell column names. If given as a list (--matrixPath input1 --matrixPath input2 etc.) then will join all matrices together. Assumes the same number and order of genes in each matrix, so only cells are added."- , 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."- , delimiter :: Maybe Char <?> "([,] | CHAR) The delimiter for the csv file 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 :: Maybe String <?> "([TfIdfNorm] | UQNorm | MedNorm | TotalMedNorm | BothNorm | LogCPMNorm | NoneNorm) Type of normalization before clustering. 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. TotalMedNorm normalized first each observation by total count then by median of non-zero counts across features. BothNorm uses both normalizations (first TotalMedNorm for all analysis then additionally TfIdfNorm during clustering). LogCPMNorm normalizes by log2(CPM + 1). 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). Cannot use TfIdfNorm for any other process as NoneNorm will become the default."- , 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."- , shiftPositive :: Bool <?> "Shift features to positive values. Positive values are shifted to allow modularity to work correctly."- , noFilter :: Bool <?> "Whether to bypass filtering genes and cells by low counts."- , filterThresholds :: Maybe String <?> "([(250, 1)] | (DOUBLE, DOUBLE)) The minimum filter thresholds for (MINCELL, MINFEATURE) when filtering cells and features by low read counts. See also --no-filter."- , 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."}- | 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 gene row names and cell column names. If given as a list (--matrixPath input1 --matrixPath input2 etc.) then will join all matrices together. Assumes the same number and order of genes in each matrix, so only cells are added."- , 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."- , 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."- , noFilter :: Bool <?> "Whether to bypass filtering genes and cells by low counts."- , shiftPositive :: Bool <?> "Shift features to positive values. Positive values are shifted to allow modularity to work correctly."- , filterThresholds :: Maybe String <?> "([(250, 1)] | (DOUBLE, DOUBLE)) The minimum filter thresholds for (MINCELL, MINFEATURE) when filtering cells and features by low read counts. See also --no-filter."- , delimiter :: Maybe Char <?> "([,] | CHAR) The delimiter for the csv file if using a normal csv rather than cellranger output and for --labels-file."- , normalization :: Maybe String <?> "([TfIdfNorm] | UQNorm | MedNorm | TotalMedNorm | BothNorm | LogCPMNorm | NoneNorm) Type of normalization before clustering. 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. TotalMedNorm normalized first each observation by total count then by median of non-zero counts across features. BothNorm uses both normalizations (first TotalMedNorm for all analysis then additionally TfIdfNorm during clustering). LogCPMNorm normalizes by log2(CPM + 1). 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). Cannot use TfIdfNorm for any other process as NoneNorm will become the default."- , 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."- , 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 (node / other cells) using the Kruskal-Wallis Test (--genes 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.)."- , 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 genes."- , genes :: [T.Text] <?> "([Nothing] | GENE) List of genes to plot for all cells within selected nodes. Invoked by --genes CD4 --genes 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 \"--genes\" instead of plotting different distributions for each feature."- , 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 csv file 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."- , pathDistance :: Maybe String <?> "([PathModularity] | PathStep) 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 csv file if using a normal csv rather than cellranger output and for --labels-file."- , output :: Maybe String <?> "([out] | STRING) The folder containing output."}- deriving (Generic)--modifiers :: Modifiers-modifiers = lispCaseModifiers { shortNameModifier = short }- where- short "aggregate" = Nothing- short "customCut" = Nothing- short "clumpinessMethod" = Just 'u'- short "clusterNormalization" = Just 'C'- short "dendrogramOutput" = Just 'U'- 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 "eigenGroup" = Just 'B'- short "featureColumn" = Nothing- short "filterThresholds" = Just 'H'- short "labels" = Nothing- short "matrixOutput" = Nothing- short "matrixOutputTranspose" = Nothing- 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 "noFilter" = Just 'F'- short "normalization" = Just 'z'- short "numEigen" = Just 'G'- short "numRuns" = Nothing- short "order" = Just 'O'- short "pca" = Just 'a'- short "plotOutput" = Nothing- short "priors" = Just 'P'- short "projectionFile" = Just 'j'- short "rootCut" = Nothing- short "shiftPositive" = Nothing- short x = firstLetter x--instance ParseRecord Options where- parseRecord = parseRecordWithModifiers modifiers---- | Notify user of limitations and error out for incompatabilities. Empty for--- now.-limitationWarningsErrors :: Options -> IO ()-limitationWarningsErrors opts = do- return ()---- | Read or return an error.-readOrErr err = fromMaybe (error err) . readMaybe---- | Normalization defaults.-getNormalization :: Options -> NormType-getNormalization opts@(MakeTree{}) =- maybe- TfIdfNorm- (readOrErr "Cannot read --normalization.")- . unHelpful- . normalization- $ opts-getNormalization opts@(Interactive{}) =- maybe- TfIdfNorm- (readOrErr "Cannot read --normalization.")- . unHelpful- . normalization- $ opts-getNormalization opts@(Differential{}) =- maybe- NoneNorm- (readOrErr "Cannot read --normalization.")- . unHelpful- . normalization- $ opts---- | Load the single cell matrix.-loadSSM :: Options -> FilePath -> IO SingleCells-loadSSM opts matrixPath' = do- fileExist <- FP.doesFileExist matrixPath'- directoryExist <- FP.doesDirectoryExist matrixPath'- compressedFileExist <- FP.doesFileExist $ matrixPath' FP.</> "matrix.mtx.gz"-- let matrixFile' =- case (fileExist, directoryExist, compressedFileExist) of- (False, False, False) -> error "\nMatrix path does not exist."- (True, False, False) ->- Left . DecompressedMatrix . MatrixFile $ matrixPath'- (False, True, False) ->- Right . DecompressedMatrix . MatrixFile $ matrixPath' FP.</> "matrix.mtx"- (False, True, True) ->- Right . CompressedMatrix . MatrixFile $ matrixPath' FP.</> "matrix.mtx.gz"- _ -> error "Cannot determine matrix pointed to, are there too many matrices here?"- genesFile' = GeneFile- $ matrixPath'- FP.</> (bool "genes.tsv" "features.tsv.gz" compressedFileExist)- cellsFile' = CellFile- $ matrixPath'- FP.</> (bool "barcodes.tsv" "barcodes.tsv.gz" compressedFileExist)- delimiter' =- Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts- featureColumn' =- FeatureColumn . fromMaybe 1 . unHelpful . featureColumn $ opts- unFilteredSc =- case matrixFile' of- (Left (DecompressedMatrix file)) ->- loadSparseMatrixDataStream delimiter' file- (Right (DecompressedMatrix file)) ->- loadCellrangerData featureColumn' genesFile' cellsFile' file- (Right (CompressedMatrix file)) ->- loadCellrangerDataFeatures featureColumn' genesFile' cellsFile' file- _ -> error "Does not supported this matrix type. See too-many-cells -h for each entry point for more information"- unFilteredSc---- | Load all single cell matrices.-loadAllSSM :: Options -> IO (Maybe SingleCells)-loadAllSSM opts = runMaybeT $ do- let matrixPaths' = unHelpful . matrixPath $ opts- cellWhitelistFile' =- fmap CellWhitelistFile . unHelpful . cellWhitelistFile $ opts- normalization' = getNormalization opts- pca' = fmap PCADim . unHelpful . pca $ opts- noFilterFlag' = NoFilterFlag . unHelpful . noFilter $ opts- shiftPositiveFlag' =- ShiftPositiveFlag . unHelpful . shiftPositive $ opts- filterThresholds' = FilterThresholds- . maybe (250, 1) read- . unHelpful- . filterThresholds- $ opts-- liftIO $ when (isJust pca' && (elem normalization' [TfIdfNorm, BothNorm])) $- hPutStrLn stderr "\nWarning: PCA (creating negative numbers) with tf-idf\- \ normalization may lead to NaNs or 0s before spectral\- \ clustering (leading to svdlibc to hang or dense SVD\- \ to error out)! Continuing..."-- mats <- MaybeT- $ if null matrixPaths'- then return Nothing- else fmap Just . mapM (loadSSM opts) $ matrixPaths'- cellWhitelist <- liftIO . sequence $ fmap getCellWhitelist cellWhitelistFile'-- let whiteListFilter Nothing = id- whiteListFilter (Just wl) = filterWhitelistSparseMat wl- unFilteredSc = mconcat mats- sc =- ( bool (filterNumSparseMat filterThresholds') id- $ unNoFilterFlag noFilterFlag'- )- . whiteListFilter cellWhitelist- $ unFilteredSc- normMat TfIdfNorm = id -- Normalize during clustering.- normMat UQNorm = uqScaleSparseMat- normMat MedNorm = medScaleSparseMat- normMat TotalMedNorm = scaleSparseMat- normMat BothNorm = scaleSparseMat- normMat LogCPMNorm = logCPMSparseMat- normMat NoneNorm = id- processMat = normMat normalization' . _matrix- processedSc = ( bool id shiftPositiveSc- $ unShiftPositiveFlag shiftPositiveFlag'- )- . (\m -> maybe m (flip pcaDenseSc m) pca')- $ sc { _matrix = processMat sc }-- -- Check for empty matrix.- when (V.null . getRowNames $ processedSc) $ error "Matrix is empty. Check --filter-thresholds, --normalization, or the input matrix for over filtering or incorrect input format."-- liftIO . mapM_ (hPutStrLn stderr) . matrixValidity $ processedSc-- return processedSc--makeTreeMain :: Options -> IO ()-makeTreeMain opts = H.withEmbeddedR defaultConfig $ do- let readOrErr err = fromMaybe (error err) . readMaybe- matrixPaths' = unHelpful . matrixPath $ opts- labelsFile' =- fmap LabelFile . unHelpful . labelsFile $ opts- prior' =- fmap PriorPath . unHelpful . prior $ opts- delimiter' =- Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts- eigenGroup' =- maybe SignGroup (readOrErr "Cannot read eigen-group.")- . unHelpful- . eigenGroup- $ opts- dense' = DenseFlag . unHelpful . dense $ opts- normalization' = 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- 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- customCut' = CustomCut . Set.fromList . unHelpful . customCut $ opts- rootCut' = fmap RootCut . unHelpful . rootCut $ opts- dendrogramOutput' = DendrogramFile- . fromMaybe "dendrogram.svg"- . unHelpful- . dendrogramOutput- $ opts- matrixOutput' = fmap (getMatrixOutputType . (unOutputDirectory output' FP.</>))- . unHelpful- . matrixOutput- $ opts- matrixOutputTranspose' = fmap (getMatrixOutputType . (unOutputDirectory output' FP.</>))- . unHelpful- . matrixOutputTranspose- $ 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- 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- drawColors' = fmap ( CustomColors- . fmap sRGB24read- . (\x -> readOrErr "Cannot read draw-colors." x :: [String])- )- . unHelpful- . drawColors- $ opts- drawDiscretize' = (=<<) (\x -> either error Just- . either- (\ err -> either- (\y -> Left $ finalError err y)- (Right . SegmentColorMap)- (readEither x :: Either String Int)- )- (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- drawFont' = fmap DrawFont . unHelpful . drawFont $ 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-- -- Increment progress bar.- Progress.autoProgressBar- (Progress.msg "Loading matrix")- Progress.percentage- 80- $ Progress.Progress 0 10-- -- Load matrix once.- processedSc <- loadAllSSM opts-- -- Notify user of limitations.- limitationWarningsErrors opts-- -- Increment progress bar.- Progress.autoProgressBar- (Progress.msg "Planning leaf colors")- Progress.percentage- 80- $ Progress.Progress 1 10-- -- 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- _ -> mapM (loadLabelData delimiter') $ labelsFile'-- -- Increment progress bar.- Progress.autoProgressBar- (Progress.msg "Planting tree")- Progress.percentage- 80- $ Progress.Progress 2 10-- --R.withEmbeddedR R.defaultConfig $ R.runRegion $ do- -- For r clustering.- -- mat <- scToRMat processedSc- -- clusterRes <- hdbscan mat- -- clusterList <- clustersToClusterList sc clusterRes-- -- For agglomerative clustering.- --let clusterResults = fmap hClust processedSc-- -- Load previous results or calculate results if first run.- originalClusterResults <- case prior' of- Nothing -> do- (fullCr, _) <-- hSpecClust dense' eigenGroup' normalization' numEigen' minModularity' numRuns'- . extractSc- $ processedSc-- return fullCr :: IO ClusterResults- (Just x) -> do- let clInput = unPriorPath x FP.</> "cluster_list.json"- treeInput = unPriorPath x FP.</> "cluster_tree.json"-- -- Strict loading in order to avoid locked file.- loadClusterResultsFiles clInput treeInput-- -- Increment progress bar.- Progress.autoProgressBar- (Progress.msg "Measuring roots")- Progress.percentage- 80- $ Progress.Progress 3 10-- let birchMat = processedSc- birchSimMat =- case (not . null . unHelpful . matrixPath $ opts, drawCollection') of- (True, CollectionGraph{} ) ->- Just- . B2Matrix- . L.over matrix (MatObsRow . unB2 . b1ToB2 . B1 . unMatObsRow)- . extractSc- $ processedSc- _ -> Nothing-- -- Increment progress bar.- Progress.autoProgressBar- (Progress.msg "Sketching tree")- Progress.percentage- 80- $ Progress.Progress 4 10-- let config :: BirchBeer.Types.Config CellInfo SingleCells- config = BirchBeer.Types.Config- { _birchLabelMap = labelMap- , _birchMinSize = minSize'- , _birchMaxStep = maxStep'- , _birchMaxProportion = maxProportion'- , _birchMinDistance = minDistance'- , _birchMinDistanceSearch = minDistanceSearch'- , _birchSmartCutoff = smartCutoff'- , _birchCustomCut = customCut'- , _birchRootCut = rootCut'- , _birchOrder = Just order'- , _birchDrawLeaf = drawLeaf'- , _birchDrawCollection = drawCollection'- , _birchDrawMark = drawMark'- , _birchDrawNodeNumber = drawNodeNumber'- , _birchDrawMaxNodeSize = drawMaxNodeSize'- , _birchDrawMaxLeafNodeSize = drawMaxLeafNodeSize'- , _birchDrawNoScaleNodes = drawNoScaleNodes'- , _birchDrawLegendSep = drawLegendSep'- , _birchDrawLegendAllLabels = drawLegendAllLabels'- , _birchDrawPalette = drawPalette'- , _birchDrawColors = drawColors'- , _birchDrawDiscretize = drawDiscretize'- , _birchDrawScaleSaturation = drawScaleSaturation'- , _birchDrawFont = drawFont'- , _birchTree = _clusterDend originalClusterResults- , _birchMat = birchMat- , _birchSimMat = birchSimMat- }-- (plot, labelColorMap, itemColorMap, markColorMap, tree', gr') <- mainDiagram config-- -- Increment progress bar.- Progress.autoProgressBar- (Progress.msg "Recording tree measurements")- Progress.percentage- 80- $ Progress.Progress 5 10-- -- Write results.- clusterResults <- case prior' of- Nothing -> do- let clusterList' = treeToClusterList tree'- cr' = ClusterResults clusterList' tree'-- return cr'-- (Just x) -> do- let clusterList' = treeToClusterList tree'- cr' = ClusterResults clusterList' tree'-- return cr'-- B.writeFile- (unOutputDirectory output' FP.</> "cluster_list.json")- . A.encode- . _clusterList- $ clusterResults- B.writeFile- (unOutputDirectory output' FP.</> "cluster_tree.json")- . A.encode- . _clusterDend- $ clusterResults- T.writeFile- (unOutputDirectory output' FP.</> "graph.dot")- . G.printDotGraph- . G.graphToDot G.nonClusteredParams- . unClusterGraph- $ gr'- B.writeFile- (unOutputDirectory output' FP.</> "cluster_info.csv")- . printClusterInfo- $ gr'- -- Write matrix- mapM_ (\x -> writeMatrixLike (MatrixTranspose False) x . extractSc $ processedSc) matrixOutput'- -- Write matrix transpose- mapM_ (\x -> writeMatrixLike (MatrixTranspose True) x . extractSc $ processedSc) matrixOutputTranspose'- -- Write node info- B.writeFile- (unOutputDirectory output' FP.</> "node_info.csv")- . printNodeInfo labelMap- $ gr'- case labelMap of- Nothing -> return ()- (Just lm) ->- -- Write cluster diversity- case clusterDiversity order' lm clusterResults of- (Left err) -> hPutStrLn stderr- $ err- <> "\nError in diversity, skipping cluster_diversity.csv output."- (Right result) ->- B.writeFile- ( unOutputDirectory output'- FP.</> "cluster_diversity.csv"- )- . printClusterDiversity- $ result-- -- Increment progress bar.- Progress.autoProgressBar- (Progress.msg "Counting leaves")- Progress.percentage- 80- $ Progress.Progress 6 10-- -- Header- B.putStrLn $ "cell,cluster,path"-- -- Body- B.putStrLn- . CSV.encode- . fmap (\ (!ci, !(c:cs))- -> ( unCell . _barcode $ ci- , showt $ unCluster c- , T.intercalate "/" . fmap (showt . unCluster) $ c:cs- )- )- . _clusterList- $ clusterResults-- -- Increment progress bar.- Progress.autoProgressBar- (Progress.msg "Painting sketches")- Progress.percentage- 80- $ Progress.Progress 7 10-- -- Plot only if needed and ignore non-tree analyses if dendrogram is- -- supplied.- H.runRegion $ do- -- Calculations with plotting the label map (clumpiness).- case labelMap of- Nothing ->- H.io $ hPutStrLn stderr "\nClumpiness requires labels for cells, skipping..."- (Just lm) -> do- -- Get clumpiness.- case treeToClumpList clumpinessMethod' lm . _clusterDend $ clusterResults of- (Left err) -> H.io- . hPutStrLn stderr- $ err- <> "\nError in clumpiness, skipping clumpiness.* output."- (Right clumpList) -> do- -- Save clumpiness to a file.- H.io- . B.writeFile (unOutputDirectory output' FP.</> "clumpiness.csv")- . clumpToCsv- $ clumpList-- -- Plot clumpiness.- either (H.io . hPutStrLn stderr) id- $ plotClumpinessHeatmapR- (unOutputDirectory output' FP.</> "clumpiness.pdf")- clumpList-- -- View cutting location for modularity.- case minDistanceSearch' of- Nothing -> return ()- (Just _) -> plotRankedModularityR- (unOutputDirectory output' FP.</> "modularity_rank.pdf")- . L.view clusterDend- $ clusterResults-- -- Increment progress bar.- H.io $ Progress.autoProgressBar- (Progress.msg "Painting tree")- Progress.percentage- 80- $ Progress.Progress 8 10-- -- Plot.- H.io $ do- -- cr <- clusterResults- -- gr <- graph- -- cm <- itemColorMap-- -- plot <- if drawDendrogram'- -- then return . plotDendrogram legend drawLeaf' cm . _clusterDend $ cr- -- else do- -- plotGraph legend drawConfig cm markColorMap gr-- D.renderCairo- ( unOutputDirectory output'- FP.</> unDendrogramFile dendrogramOutput'- )- (D.mkHeight 1000)- plot-- -- Increment progress bar.- H.io $ Progress.autoProgressBar- (Progress.msg "Squishing tree")- Progress.percentage- 80- $ Progress.Progress 9 10-- -- Plot clustering of the projections.- case projectionFile' of- Nothing -> return ()- (Just f) -> do- -- Load projection map if provided.- projectionMap <- H.io $ loadProjectionMap f-- plotClustersR- (unOutputDirectory output' FP.</> "projection.pdf")- projectionMap- . _clusterList- $ clusterResults-- -- Plot clustering with labels.- case (labelMap, itemColorMap) of- (Just lm, Just icm) ->- plotLabelClustersR- (unOutputDirectory output' FP.</> "label_projection.pdf")- projectionMap- lm- icm- (_clusterList clusterResults)- _ -> return ()-- -- Increment progress bar.- H.io $ Progress.autoProgressBar- (Progress.msg "Packing up")- Progress.percentage- 80- $ Progress.Progress 10 10-- return ()---- | Interactive tree interface.-interactiveMain :: Options -> IO ()-interactiveMain opts = H.withEmbeddedR defaultConfig $ do- let labelsFile' =- fmap LabelFile . unHelpful . labelsFile $ opts- prior' = maybe (error "\nRequires --prior") PriorPath- . unHelpful- . prior- $ opts- delimiter' =- Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts- normalization' = getNormalization opts-- mat <- loadAllSSM opts- labelMap <- sequence . fmap (loadLabelData delimiter') $ labelsFile'-- tree <- fmap (either error id . A.eitherDecode)- . B.readFile- . (FP.</> "cluster_tree.json")- . unPriorPath- $ prior' :: IO (Tree (TreeNode (V.Vector CellInfo)))-- interactiveDiagram- tree- labelMap- mat- . fmap ( B2Matrix- . L.over matrix (MatObsRow . unB2 . b1ToB2 . B1 . unMatObsRow)- )- $ mat-- return ()---- | Differential path.-differentialMain :: Options -> IO ()-differentialMain opts = do- let readOrErr err = fromMaybe (error err) . readMaybe- delimiter' =- Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts- labelsFile' =- fmap LabelFile . unHelpful . labelsFile $ 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- genes' = fmap Gene . unHelpful . genes $ opts- aggregate' = Aggregate . unHelpful . aggregate $ opts- labels' = fmap ( DiffLabels- . L.over L.both ( (\x -> bool (Just x) Nothing . Set.null $ x)- . Set.fromList- . fmap Label- )- . readOrErr "Cannot read --labels."- )- . unHelpful- . labels- $ opts- (combined1, combined2) = combineNodesLabels nodes' labels'- plotOutputR = fromMaybe "out.pdf" . unHelpful . plotOutput $ opts-- processedSc <- loadAllSSM opts-- labelMap <- mapM (loadLabelData delimiter') $ labelsFile'-- let clInput = (FP.</> "cluster_list.json") . unPriorPath $ prior'- treeInput = (FP.</> "cluster_tree.json") . unPriorPath $ prior'-- cr :: IO ClusterResults- cr = loadClusterResultsFiles clInput treeInput-- gr <- fmap (treeToGraph . _clusterDend) cr-- H.withEmbeddedR defaultConfig $ H.runRegion $ do- case genes' of- [] -> do- case nodes' of- (DiffNodes ([], [])) -> do- let res = getAllDEGraphKruskalWallis- topN'- labelMap- (maybe (DiffLabels (Nothing, Nothing)) id labels')- (extractSc processedSc)- $ gr-- H.io . B.putStrLn . getAllDEStringKruskalWallis $ res- (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- res <- getDEGraph- topN'- labelMap- (extractSc processedSc)- combined1- combined2- gr-- H.io . B.putStrLn . getDEString $ res- _ -> do- let outputCsvR = FP.replaceExtension plotOutputR ".csv"-- diffPlot <- getSingleDiff- False- aggregate'- labelMap- (extractSc processedSc)- combined1- combined2- genes'- gr- [r| suppressMessages(write.csv(diffPlot_hs[[2]], file = outputCsvR_hs, row.names = FALSE, quote = FALSE)) |]- [r| suppressMessages(ggsave(diffPlot_hs[[1]], file = plotOutputR_hs)) |]-- let normOutputR = FP.replaceBaseName- plotOutputR- (FP.takeBaseName plotOutputR <> "_scaled")- normOutputCsvR = FP.replaceExtension normOutputR ".csv"-- diffNormPlot <- getSingleDiff- True- aggregate'- labelMap- (extractSc processedSc)- combined1- combined2- genes'- gr- [r| suppressMessages(write.csv(diffNormPlot_hs[[2]], file = normOutputCsvR_hs, row.names = FALSE, quote = FALSE)) |]- [r| suppressMessages(ggsave(diffNormPlot_hs[[1]], file = normOutputR_hs)) |]-- return ()---- | Diversity path.-diversityMain :: Options -> IO ()-diversityMain opts = do- let priors' =- fmap PriorPath . unHelpful . priors $ opts- delimiter' =- Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts- labelsFile' =- fmap LabelFile . unHelpful . labelsFile $ 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-- -- Where to place output files.- FP.createDirectoryIfMissing True . unOutputDirectory $ output'-- labelMap <- sequence . fmap (loadLabelData delimiter') $ labelsFile'-- pops <- fmap ( either- (\err -> error $ err <> "\nEncountered error in population loading, aborting process.")- id- . sequence- )- . mapM (\x -> (fmap . fmap) (Label . T.pack . unPriorPath $ x,)- . loadPopulation labelMap- $ x- )- $ priors'-- popDiversities <-- mapM- (\ (l, pop) -> getPopulationDiversity- l- order'- start'- interval'- endMay'- pop- )- pops-- -- Output quantifications.- B.writeFile (unOutputDirectory output' FP.</> "diversity.csv")- . CSV.encodeDefaultOrderedByName- $ popDiversities-- -- D.renderCairo (unOutputDirectory output' FP.</> "diversity.pdf") D.absolute- -- . plotDiversity- -- $ popDiversities-- -- D.renderCairo (unOutputDirectory output' FP.</> "chao1.pdf") D.absolute- -- . plotChao1- -- $ popDiversities-- -- D.renderCairo (unOutputDirectory output' FP.</> "rarefaction.pdf") D.absolute- -- . plotRarefaction- -- $ popDiversities-- -- Output plots.- let colors = D.colorRamp (length pops) . D.brewerSet D.Set1 $ 9-- H.withEmbeddedR defaultConfig $ H.runRegion $ do- let divFile = unOutputDirectory output' FP.</> "diversity.pdf"- divPlot <- plotDiversityR colors popDiversities- [r| suppressMessages(ggsave(divPlot_hs, file = divFile_hs)) |]-- -- let chao1File = unOutputDirectory output' FP.</> "chao_r.pdf"- -- chao1Plot <- plotChao1R colors popDiversities- -- [r| suppressMessages(ggsave(chao1Plot_hs, file = chao1File_hs)) |]-- let rarefactionFile = unOutputDirectory output' FP.</> "rarefaction.pdf"- rarefactionPlot <- plotRarefactionR colors popDiversities- [r| suppressMessages(ggsave(rarefactionPlot_hs, file = rarefactionFile_hs)) |]-- return ()-- return ()---- | Paths path.-pathsMain :: Options -> IO ()-pathsMain opts = do- let labelsFile' =- maybe (error "\nNeed a label file.") LabelFile- . unHelpful- . labelsFile- $ opts- prior' =- maybe (error "\nNeed a prior path containing tree.") PriorPath- . unHelpful- . prior- $ opts- delimiter' =- Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts- bandwidth' = Bandwidth . fromMaybe 1 . unHelpful . bandwidth $ opts- direction' = FlipFlag . unHelpful . flipDirection $ opts- pathDistance' =- maybe PathStep read . unHelpful . pathDistance $ opts- output' =- OutputDirectory . fromMaybe "out" . unHelpful . output $ opts-- -- Where to place output files.- FP.createDirectoryIfMissing True . unOutputDirectory $ output'-- -- Get the label map from a file.- labelMap <- loadLabelData delimiter' labelsFile'-- -- Load previous results or calculate results if first run.- tree <- fmap (either error id . A.eitherDecode)- . B.readFile- . (FP.</> "cluster_tree.json")- . unPriorPath- $ prior'-- let gr = treeToGraph tree- pathDistances :: [(CellInfo, Double)]- pathDistances = linearItemDistance direction' pathDistance' gr- labeledPathDistances =- labelItemDistance labelMap pathDistances-- H.withEmbeddedR defaultConfig $ H.runRegion $ do- plotPathDistanceR- (unOutputDirectory output' FP.</> "path_distances.pdf")- bandwidth'- labeledPathDistances- return ()-- return ()--main :: IO ()-main = do- opts <- getRecord "too-many-cells, Gregory W. Schwartz.\- \ Clusters and analyzes single cell data."-- case opts of- MakeTree{} -> makeTreeMain opts- Interactive{} -> interactiveMain opts- Differential{} -> differentialMain opts- Main.Diversity{} -> diversityMain opts- Paths{} -> pathsMain opts+{-# LANGUAGE OverloadedStrings #-}++module Main where++-- Remote+import Options.Generic++-- Local+import TooManyCells.Program.Differential+import TooManyCells.Program.Diversity+import TooManyCells.Program.Interactive+import TooManyCells.Program.MakeTree+import TooManyCells.Program.MatrixOutput+import TooManyCells.Program.Classify+import TooManyCells.Program.Motifs+import TooManyCells.Program.Options+import TooManyCells.Program.Paths+import TooManyCells.Program.Peaks++main :: IO ()+main = do+ opts <- getRecord "too-many-cells, Gregory W. Schwartz.\+ \ Clusters and analyzes single cell data."++ 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
+ src/TooManyCells/Classify/Classify.hs view
@@ -0,0 +1,86 @@+{- TooManyCells.Classify.Classify+Gregory W. Schwartz++Collects functions pertaining to classifying cells from reference populations.+-}++{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE BangPatterns #-}++module TooManyCells.Classify.Classify+ ( classifyCells+ , classifyMat+ ) where++-- Remote+import BirchBeer.Types (getMatrix, Label (..))+import Data.Function (on)+import Data.List (maximumBy)+import qualified Control.Lens as L+import qualified Data.Sparse.Common as S+import qualified Data.Vector as V+import qualified Math.Clustering.Spectral.Sparse as Spectral++-- Local+import TooManyCells.Matrix.Types+import TooManyCells.Matrix.Utility+import TooManyCells.MakeTree.Adjacency++-- | Classify cells in a SingleCells type from a list of references.+classifyCells :: SingleCells -> [AggReferenceMat] -> [(Cell, (Label, Double))]+classifyCells sc refs = zip (V.toList . L.view rowNames $ sc')+ . fmap (getMaxLabel . zip refNames . S.toDenseListSV)+ . S.toRowsL+ $ classifyMat sc' refs'+ where+ (sc', refs') = unifyFeatures sc refs+ refNames = fmap (Label . unCell)+ . V.toList+ . L.view rowNames+ . unAggSingleCells+ $ unAggReferenceMat refs'+ getMaxLabel = maximumBy (compare `on` snd)++-- | Classify matrices using cosine similarity. Results in a matrix with scores+-- (cosine similarity) for each column and observations per row.+classifyMat :: SingleCells -> AggReferenceMat -> S.SpMatrix Double+classifyMat sc refs = scMat S.## S.transpose refMat+ where+ refMat = normScMat . unAggSingleCells . unAggReferenceMat $ refs+ scMat = normScMat sc+ normScMat = Spectral.unB . Spectral.b2ToB . Spectral.B2 . getMatrix++-- | Force unification of features with reference features.+unifyFeatures :: SingleCells -> [AggReferenceMat] -> (SingleCells, AggReferenceMat)+unifyFeatures sc refs = (newSc, newRefs)+ where+ newSc = L.over rowNames (V.take nRowsSc)+ . L.over matrix (extractObsRows 0 (nRowsSc - 1))+ $ unified+ newRefs = AggReferenceMat+ . AggSingleCells+ . L.over rowNames (V.drop nRowsSc)+ . L.over matrix (extractObsRows nRowsSc (nRows - 1))+ $ unified+ nRefs = length refs+ (nRows, nCols) = S.dim . getMatrix $ unified+ (nRowsSc, nColsSc) = S.dim . getMatrix $ sc+ unified = mconcat+ $ sc : fmap (unAggSingleCells . unAggReferenceMat) refs+ extractObsRows lb ub (MatObsRow x) =+ MatObsRow $ S.extractSubmatrixRebalanceKeys x (lb, ub) (0, nCols - 1)++-- | Get the name of a reference from an AggReferenceMat.+getRefName :: AggReferenceMat -> Label+getRefName = maybe (error "unifyFeatures: no reference data") (Label . unCell)+ . flip (V.!?) 0+ . L.view rowNames+ . unAggSingleCells+ . unAggReferenceMat++-- | Get the reference vector from an AggReferenceMat.+getRefVec :: AggReferenceMat -> S.SpVector Double+getRefVec = flip S.extractRow 0+ . getMatrix+ . unAggSingleCells+ . unAggReferenceMat
+ src/TooManyCells/Classify/Types.hs view
@@ -0,0 +1,18 @@+{- TooManyCells.Classify.Types+Gregory W. Schwartz++Collects the types used in classification of cells.+-}++{-# LANGUAGE StrictData #-}++module TooManyCells.Classify.Types where++-- Remote++-- Local++-- Basic+newtype SingleMatrixFlag = SingleMatrixFlag+ { unSingleMatrixFlag :: Bool+ } deriving (Read,Show)
src/TooManyCells/Differential/Differential.hs view
@@ -12,7 +12,9 @@ module TooManyCells.Differential.Differential ( scToTwoD , getDEGraph+ , getDEGraphKruskalWallis , getDEString+ , getDEStringKruskalWallis , getSingleDiff , combineNodesLabels , getAllDEGraphKruskalWallis@@ -20,11 +22,12 @@ ) where -- Remote+import Data.Bool (bool) import BirchBeer.Types import BirchBeer.Utility (getGraphLeaves, getGraphLeafItems) import Control.Monad (join, mfilter) import Data.Function (on)-import Data.List (sort, sortBy, genericLength)+import Data.List (sort, sortBy, groupBy, genericLength, partition, foldl') import Data.Maybe (fromMaybe, catMaybes, isJust) import Data.Monoid ((<>)) import Language.R as R@@ -34,18 +37,22 @@ import qualified "differential" Differential as Diff import qualified "differential" Plot as Diff import qualified "differential" Types as Diff+import qualified Control.Concurrent.Async as Async import qualified Control.Foldl as Fold import qualified Control.Lens as L import qualified Data.ByteString.Lazy.Char8 as B import qualified Data.Csv as CSV import qualified Data.Foldable as F import qualified Data.Graph.Inductive as G-import qualified Data.Map as Map+import qualified Data.Map.Strict as Map import qualified Data.Set as Set import qualified Data.Sparse.Common as S 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 System.Random.MWC as MWC+import qualified System.Random.MWC.Distributions as MWC -- Local import TooManyCells.Differential.Types@@ -57,26 +64,27 @@ scToTwoD cellGroups sc = Diff.TwoDMat rNames cNames statuses nRows nCols . S.toListSM $ filteredMat where- rNames = fmap (Diff.Name . unGene) . V.toList . _colNames $ sc+ rNames = fmap (Diff.Name . unFeature) . V.toList . _colNames $ sc cNames = fmap (Diff.Name . unCell . L.view L._2) cellGroups -- We flip row and column because cells are columns here statuses = fmap (Diff.Status . showt . L.view (L._3 . L._1)) cellGroups nRows = S.nrows filteredMat nCols = S.ncols filteredMat filteredMat = S.fromColsL -- Here the columns should be observations. . fmap (S.extractRow (unMatObsRow . _matrix $ sc) . L.view L._1)+ . filter ((>) (S.nrows . unMatObsRow $ _matrix sc) . L.view L._1) $ cellGroups --- | Get the indices and statuses for two lists of nodes.+-- | Get the indices and statuses for a list of groups of nodes. getStatuses :: Maybe LabelMap- -> ([G.Node], Maybe (Set.Set Label))- -> ([G.Node], Maybe (Set.Set Label))+ -> [([G.Node], Maybe (Set.Set Label))] -> ClusterGraph CellInfo -> [(Int, Cell, (Int, Diff.Status))]-getStatuses lm (v1, l1) (v2, l2) (ClusterGraph gr) =+getStatuses lm gs (ClusterGraph gr) = sort- . F.toList- $ mappend (collapseStatus (1 :: Int) v1 l1) (collapseStatus (2 :: Int) v2 l2)+ . concatMap F.toList+ . zipWith (\x (a, b) -> collapseStatus x a b) [1..]+ $ gs where collapseStatus s vs ls = fmap (\ !x -> (unRow . _cellRow $ x, _barcode x, (s, Diff.Status $ statusName vs ls)))@@ -89,6 +97,18 @@ statusName vs (Just ls) = (T.intercalate " " . fmap unLabel . Set.toAscList $ ls) <> " " <> showt vs +-- | Get the indices and statuses for a list of groups of nodes and subsample if+-- desired.+subsampleGetStatuses+ :: Seed+ -> Maybe Subsample+ -> Maybe LabelMap+ -> [([G.Node], Maybe (Set.Set Label))]+ -> ClusterGraph CellInfo+ -> IO [(Int, Cell, (Int, Diff.Status))]+subsampleGetStatuses seed subN lm gs gr =+ maybe return (subsampleGroups seed) subN $ getStatuses lm gs gr+ -- | Filter barcodes by labels. validCellInfo :: Maybe LabelMap -> Maybe (Set.Set Label) -> CellInfo -> Bool validCellInfo Nothing _ = const True@@ -111,77 +131,136 @@ . unCell -- | Get the differential expression of two sets of cells, filtered by labels.-getDEGraph :: TopN+getDEGraph :: Seed+ -> Maybe Subsample+ -> TopN -> Maybe LabelMap -> SingleCells -> ([G.Node], Maybe (Set.Set Label)) -> ([G.Node], Maybe (Set.Set Label)) -> ClusterGraph CellInfo -> R.R s [(Diff.Name, Double, Diff.PValue, Diff.FDR)]-getDEGraph (TopN topN) lm sc v1 v2 gr = do- let cellGroups = getStatuses lm v1 v2 gr- mat = scToTwoD cellGroups sc+getDEGraph seed subN (TopN topN) lm sc v1 v2 gr = do+ cellGroups <- H.io $ subsampleGetStatuses seed subN lm [v1, v2] gr + let mat = scToTwoD cellGroups sc+ Diff.edgeR topN mat --- | Get the differential expression of each cluster to each other cluster using--- KruskalWallis.-getAllDEGraphKruskalWallis :: TopN- -> Maybe LabelMap- -> DiffLabels- -> SingleCells- -> ClusterGraph CellInfo- -> [(G.Node, Gene, Diff.Log2Diff, Maybe Diff.PValue, Maybe Diff.FDR)]-getAllDEGraphKruskalWallis topN lm ls sc gr =+-- | Get the differential expression using Kruskall-Wallis of two sets of cells,+-- filtered by labels.+getDEGraphKruskalWallis+ :: Seed+ -> Maybe Subsample+ -> TopN+ -> Maybe LabelMap+ -> SingleCells+ -> ([G.Node], Maybe (Set.Set Label))+ -> ([G.Node], Maybe (Set.Set Label))+ -> ClusterGraph CellInfo+ -> IO [ ( Feature+ , Diff.Log2Diff+ , Maybe Diff.PValue+ , Maybe Diff.FDR+ , Maybe Diff.QValue+ )+ ]+getDEGraphKruskalWallis seed subN (TopN topN) lm sc v1 v2 gr = do+ cellGroups <- subsampleGetStatuses seed subN lm [v1, v2] gr++ let fastFiveCheck = (< 5) . length . take 5+ res = filter (isJust . L.view L._3)+ . zipWith+ (\name (!a, !b, !c, !d) -> (name, a, b, c, d))+ (V.toList . L.view colNames $ sc)+ . Diff.differentialMatrixFeatRow as bs+ . S.transpose+ . unMatObsRow+ . L.view matrix+ $ sc+ (as, bs) = L.over L.both (fmap (L.view L._1))+ . partition ((== 1) . L.view (L._3 . L._1))+ $ cellGroups+ if fastFiveCheck as || fastFiveCheck bs+ then error "Less than five cells in one node to compare."+ else return . take topN . sortBy (compare `on` L.view L._3) $ res++-- | Get the differential expression of each cluster to all other cells in the+-- data set using KruskalWallis.+getAllDEGraphKruskalWallis+ :: Seed+ -> Maybe Subsample+ -> TopN+ -> Maybe LabelMap+ -> DiffLabels+ -> SingleCells+ -> ClusterGraph CellInfo+ -> IO [(G.Node, Feature, Diff.Log2Diff, Maybe Diff.PValue, Maybe Diff.FDR, Maybe Diff.QValue)]+getAllDEGraphKruskalWallis seed subN topN lm ls sc gr = mconcat . catMaybes- . parMap rdeepseq (\n -> compareClusterToOthersKruskalWallis n topN lm ls sc mat gr)- $ nodes+ <$> Async.mapConcurrently (\n -> compareClusterToOthersKruskalWallis n seed subN topN lm ls sc mat gr)+ nodes where nodes = filter (/= 0) . G.nodes . unClusterGraph $ gr -- Don't want to look at root. mat = S.transpose . unMatObsRow . L.view matrix $ sc --- | Get the differential expression of a cluster (n) to all other clusters (ns)--- using KruskalWallis such that n / ns.+-- | Get the differential expression of a cluster (n) to all other cells in the+-- data set (ns) using KruskalWallis such that n / ns. compareClusterToOthersKruskalWallis :: G.Node+ -> Seed+ -> Maybe Subsample -> TopN -> Maybe LabelMap -> DiffLabels -> SingleCells -> S.SpMatrix Double -> ClusterGraph CellInfo- -> Maybe [(G.Node, Gene, Diff.Log2Diff, Maybe Diff.PValue, Maybe Diff.FDR)]-compareClusterToOthersKruskalWallis n (TopN topN) lm (DiffLabels (ls1, ls2)) sc mat gr- | fastFiveCheck nCells || fastFiveCheck nsCells = Nothing- | otherwise = Just . take topN . sortBy (compare `on` (L.view L._4)) $ res- where- fastFiveCheck = (< 5) . length . take 5- nCells' = F.toList $ getGraphLeafItems gr n- nCellsSet = Set.fromList . fmap (L.view barcode) $ nCells'- nCells = fmap (unRow . L.view cellRow)- . filter (validCellInfo lm ls2)- $ nCells' -- All cells from node and labels- nsCells =- fmap fst- . filter (\ (_, !x) -> validCell lm ls1 x && not (Set.member x nCellsSet)) -- All cells outside of node and from labels- . zip [0..]- . V.toList- . L.view rowNames- $ sc- res = filter (isJust . L.view L._4)- . ( zipWith- (\name (!x, !y, !z) -> (n, name, x, y, z))- (V.toList . L.view colNames $ sc)- )- $ Diff.differentialMatrixFeatRow nsCells nCells mat -- Here the matrix rows are features+ -> IO (Maybe [(G.Node, Feature, Diff.Log2Diff, Maybe Diff.PValue, Maybe Diff.FDR, Maybe Diff.QValue)])+compareClusterToOthersKruskalWallis n (Seed seed) subN (TopN topN) lm (DiffLabels (ls1, ls2)) sc mat gr = do+ let fastFiveCheck = (< 5) . length . take 5+ nCells' = F.toList $ getGraphLeafItems gr n+ nCellsSet = Set.fromList . fmap (L.view barcode) $ nCells'+ nCells = fmap (unRow . L.view cellRow)+ . filter (validCellInfo lm ls2)+ $ nCells' -- All cells from node and labels+ nsCells =+ fmap fst+ . filter (\ (_, !x) -> validCell lm ls1 x && not (Set.member x nCellsSet)) -- All cells outside of node and from labels+ . zip [0..]+ . V.toList+ . L.view rowNames+ $ sc+ sample n x =+ maybe+ return+ (\s -> fmap (take s . V.toList) . flip MWC.uniformShuffle x . V.fromList)+ n+ subN' = bool (min (length nCells) (length nsCells)) (maybe 0 unSubsample subN)+ . (/= Subsample 0)+ <$> subN+ g <- MWC.restore . MWC.toSeed $ V.fromList [fromIntegral seed]+ nCellsFinal <- sample subN' g nCells+ nsCellsFinal <- sample subN' g nsCells + let res = filter (isJust . L.view L._4)+ . ( zipWith+ (\name (!a, !b, !c, !d) -> (n, name, a, b, c, d))+ (V.toList . L.view colNames $ sc)+ )+ $ Diff.differentialMatrixFeatRow nsCellsFinal nCellsFinal mat -- Here the matrix rows are features++ if fastFiveCheck nCellsFinal || fastFiveCheck nsCellsFinal+ then return Nothing+ else return . Just . take topN . sortBy (compare `on` (L.view L._4)) $ res+ -- | Get the differential expression of two sets of cells. getDEString :: [(Diff.Name, Double, Diff.PValue, Diff.FDR)] -> B.ByteString getDEString xs = header <> "\n" <> body where- header = "gene,logFC,pVal,FDR"+ header = "feature,log2FC,pVal,FDR" body = CSV.encode . fmap ( L.over L._1 Diff.unName . L.over L._3 Diff.unPValue@@ -192,68 +271,108 @@ -- | Get the differential expression of each node to all other nodes using -- KruskalWallis. getAllDEStringKruskalWallis- :: [(G.Node, Gene, Diff.Log2Diff, Maybe Diff.PValue, Maybe Diff.FDR)] -> B.ByteString+ :: [(G.Node, Feature, Diff.Log2Diff, Maybe Diff.PValue, Maybe Diff.FDR, Maybe Diff.QValue)]+ -> B.ByteString getAllDEStringKruskalWallis xs = header <> "\n" <> body where- header = "node,gene,log2FC,pVal,FDR"+ header = "node,feature,log2FC,pVal,FDR" body = CSV.encode- . fmap ( L.over L._5 (maybe "NA" (showt . Diff.unFDR))+ . fmap ( L.over L._6 (maybe "NA" (showt . Diff.unQValue))+ . L.over L._5 (maybe "NA" (showt . Diff.unFDR)) . L.over L._4 (maybe "NA" (showt . Diff.unPValue)) . L.over L._3 Diff.unLog2Diff- . L.over L._2 unGene+ . L.over L._2 unFeature ) $ xs +-- | Get the differential expression string between two sets of nodes using+-- KruskalWallis.+getDEStringKruskalWallis+ :: [(Feature, Diff.Log2Diff, Maybe Diff.PValue, Maybe Diff.FDR, Maybe Diff.QValue)]+ -> B.ByteString+getDEStringKruskalWallis xs = header <> "\n" <> body+ where+ header = "feature,log2FC,pVal,FDR,qVal"+ body = CSV.encode+ . fmap ( 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+ . L.over L._1 unFeature+ )+ $ xs+ -- | Convert a single cell matrix to a list of Entities with the specified--- features. Also aggregates genes by average value or not.+-- features. Also aggregates features by average value or not. scToEntities :: Aggregate- -> [Gene]+ -> [Feature] -> [(Int, Cell, (Int, Diff.Status))] -> SingleCells -> [Diff.Entity]-scToEntities aggregate genes cellGroups sc =- concatMap (\x -> toEntity aggregate x geneIdxs) cellGroups+scToEntities aggregate features cellGroups sc =+ concatMap (\x -> toEntity aggregate x featureIdxs) cellGroups where mat = getMatrix sc toEntity (Aggregate False) (cellIdx, (Cell b), (_, status)) =- fmap (\ (Gene gene, idx) -> Diff.Entity (Diff.Name gene) status (Diff.Id b)+ fmap (\ (Feature feature, idx) -> Diff.Entity (Diff.Name feature) status (Diff.Id b) $ S.lookupWD_SM mat (cellIdx, idx) ) toEntity (Aggregate True) (cellIdx, (Cell b), (_, status)) = (:[]) . Diff.Entity- (Diff.Name . T.intercalate " " . fmap (unGene . fst) $ geneIdxs)+ (Diff.Name . T.intercalate " " . fmap (unFeature . fst) $ featureIdxs) status (Diff.Id b) . (/ n)- . sum+ . foldl' (+) 0 . fmap (\(_, idx) -> S.lookupWD_SM mat (cellIdx, idx))- n = genericLength geneIdxs- geneIdxs :: [(Gene, Int)]- geneIdxs = fmap (\ !x -> ( x+ n = genericLength featureIdxs+ featureIdxs :: [(Feature, Int)]+ featureIdxs = fmap (\ !x -> ( x , fromMaybe (err x)- $ V.elemIndex (unGene x) (getColNames sc)+ $ V.elemIndex (unFeature x) (getColNames sc) )- ) genes+ ) features err x = error $ "Feature " <> show x <> " not found for differential." -- | Get the differential expression plot of features (or aggregate of features -- by average) over statuses, filtered by labels.-getSingleDiff :: Bool+getSingleDiff :: Seed+ -> Maybe Subsample+ -> Bool+ -> ViolinFlag+ -> NoOutlierFlag -> Aggregate+ -> SeparateNodes+ -> SeparateLabels -> Maybe LabelMap -> SingleCells -> ([G.Node], Maybe (Set.Set Label)) -> ([G.Node], Maybe (Set.Set Label))- -> [Gene]+ -> [Feature] -> ClusterGraph CellInfo -> R.R s (R.SomeSEXP s)-getSingleDiff normalize aggregate lm sc v1 v2 genes gr = do- let cellGroups = getStatuses lm v1 v2 gr- entities = scToEntities aggregate genes cellGroups sc+getSingleDiff seed subN normalize (ViolinFlag vf) (NoOutlierFlag noOutlierF) aggregate sn sl lm sc v1 v2 features gr = do+ let splitNodeGroup (!ns, !ls) = fmap (\ !x -> ([x], ls)) ns+ splitLabelGroup (!ns, !ls) =+ maybe+ [(ns, ls)]+ (fmap (\ !l -> (ns, Just $ Set.singleton l)) . Set.toAscList)+ ls+ groupsAssign' = case (unSeparateNodes sn, unSeparateLabels sl) of+ (False, False) -> [v1, v2]+ (True, False) -> concatMap splitNodeGroup [v1, v2]+ (False, True) -> concatMap splitLabelGroup [v1, v2]+ (True, True) -> concatMap splitNodeGroup+ . concatMap splitLabelGroup+ $ [v1, v2] - Diff.plotSingleDiff normalize entities+ cellGroups <- H.io $ subsampleGetStatuses seed subN lm groupsAssign' gr + let entities = scToEntities aggregate features cellGroups sc++ Diff.plotSingleDiff normalize vf noOutlierF entities+ -- | Combine nodes and labels. combineNodesLabels :: DiffNodes@@ -262,3 +381,23 @@ combineNodesLabels (DiffNodes (v1, v2)) Nothing = ((v1, Nothing), (v2, Nothing)) combineNodesLabels (DiffNodes (v1, v2)) (Just (DiffLabels (l1, l2))) = ((v1, l1), (v2, l2))++-- | Subsample a cell group list with provided subsampling number or the+-- smallest of the two groups if not provided.+subsampleGroups :: Seed+ -> Subsample+ -> [(Int, Cell, (Int, Diff.Status))]+ -> IO [(Int, Cell, (Int, Diff.Status))]+subsampleGroups (Seed seed) subN xs = do+ g <- MWC.restore . MWC.toSeed $ V.fromList [fromIntegral seed]++ let sample n x = fmap (take n . V.toList) . flip MWC.uniformShuffle x+ grouped = fmap V.fromList+ . groupBy ((==) `on` L.view (L._3 . L._1))+ . sortBy (compare `on` L.view (L._3 . L._1))+ $ xs+ subN' = bool (minimum . fmap V.length $ grouped) (unSubsample subN)+ . (/= 0)+ $ unSubsample subN++ fmap mconcat $ mapM (sample subN' g) grouped
src/TooManyCells/Differential/Types.hs view
@@ -17,7 +17,14 @@ -- Basic newtype TopN = TopN { unTopN :: Int }+newtype NoEdger = NoEdger { unNoEdger :: Bool } newtype DiffNodes = DiffNodes {unDiffNodes :: ([G.Node], [G.Node])} newtype DiffLabels = DiffLabels { unDiffLabels :: (Maybe (Set.Set Label), Maybe (Set.Set Label)) } newtype Aggregate = Aggregate { unAggregate :: Bool }+newtype SeparateNodes = SeparateNodes { unSeparateNodes :: Bool }+newtype SeparateLabels = SeparateLabels { unSeparateLabels :: Bool }+newtype ViolinFlag = ViolinFlag { unViolinFlag :: Bool }+newtype NoOutlierFlag = NoOutlierFlag { unNoOutlierFlag :: Bool }+newtype Seed = Seed { unSeed :: Int }+newtype Subsample = Subsample { unSubsample :: Int } deriving (Eq)
src/TooManyCells/Diversity/Plot.hs view
@@ -96,6 +96,7 @@ ylab("Diversity") + scale_fill_manual(values = as.character(colorsR_hs)) + guides(fill = "none") ++ theme_cowplot() + theme(aspect.ratio = 1, axis.text.x = element_text(angle = 315, hjust = 0)) |] @@ -117,6 +118,7 @@ ylab("Chao1") + scale_fill_manual(values = as.character(colorsR_hs)) + guides(fill = "none") ++ theme_cowplot() + theme(aspect.ratio = 1, axis.text.x = element_text(angle = 315, hjust = 0)) |] @@ -149,5 +151,6 @@ ylab("Estimated richness") + scale_color_manual(values = as.character(colorsR_hs)) + guides(color = guide_legend(title = "")) ++ theme_cowplot() + theme(aspect.ratio = 0.5) |]
src/TooManyCells/File/Types.hs view
@@ -36,7 +36,7 @@ -- Basic newtype DendrogramFile = DendrogramFile { unDendrogramFile :: FilePath } newtype CellFile = CellFile { unCellFile :: FilePath }-newtype GeneFile = GeneFile { unGeneFile :: FilePath }+newtype FeatureFile = FeatureFile { unFeatureFile :: FilePath } newtype ProjectionFile = ProjectionFile { unProjectionFile :: FilePath } newtype CellWhitelistFile = CellWhitelistFile { unCellWhitelistFile :: FilePath@@ -48,5 +48,12 @@ -- Advanced data MatrixFileFolder = MatrixFile FilePath | MatrixFolder FilePath+ deriving (Read, Show)+data FragmentsFile = FragmentsFile FilePath+ deriving (Read, Show) data MatrixFileType = DecompressedMatrix MatrixFileFolder | CompressedMatrix MatrixFileFolder+ | CompressedFragments FragmentsFile+ | BedGraph FilePath+ | BigWig FilePath+ deriving (Read, Show)
src/TooManyCells/MakeTree/Adjacency.hs view
@@ -13,6 +13,7 @@ ) where -- Remote+import Data.List (foldl') import qualified Numeric.LinearAlgebra as H import qualified Data.Sparse.Common as S @@ -27,10 +28,9 @@ -- | Get the cosine similarity between two vectors. cosineSimilaritySparse :: S.SpVector Double -> S.SpVector Double -> Double-cosineSimilaritySparse v w = dot v w / (norm2 v * norm2 w)+cosineSimilaritySparse v w = S.dot v w / (norm2 v * norm2 w) where- dot x = sum . S.liftU2 (*) x- norm2 = sqrt . sum . fmap (** 2)+ norm2 = sqrt . foldl' (+) 0 . fmap (** 2) -- | Get an adjacency matrix based on a matrix where each row is an observation -- and the adjacencies are cosine similarities.
src/TooManyCells/MakeTree/Cluster.hs view
@@ -24,7 +24,7 @@ import BirchBeer.Utility (getGraphLeaves, getGraphLeavesWithParents, dendrogramToGraph, dendToTree, clusteringTreeToTree, treeToGraph) import Control.Monad (join) import Data.Function (on)-import Data.List (sortBy, groupBy, zip4, genericLength)+import Data.List (sortBy, groupBy, zip4, genericLength, foldl') import Data.Int (Int32) import Data.Maybe (fromMaybe, catMaybes, mapMaybe) import Math.Modularity.Types (Q (..))@@ -37,7 +37,7 @@ import Math.Clustering.Hierarchical.Spectral.Sparse (hierarchicalSpectralCluster, B (..)) import Math.Clustering.Hierarchical.Spectral.Types (clusteringTreeToDendrogram, getClusterItemsDend, EigenGroup (..)) import Math.Diversity.Diversity (diversity)-import Statistics.Quantile (continuousBy, s)+import Statistics.Quantile (quantile, s) import System.IO (hPutStrLn, stderr) import Safe (headMay) import TextShow (showt)@@ -93,7 +93,7 @@ $ dend dend = HC.dendrogram HC.CLINK items euclDist euclDist x y =- sqrt . sum . fmap (** 2) $ S.liftU2 (-) (L.view L._2 y) (L.view L._2 x)+ sqrt . foldl' (+) 0 . fmap (** 2) $ (L.view L._2 y) S.^-^ (L.view L._2 x) items = (\ fs -> zip3 (V.toList $ _rowNames sc)@@ -113,7 +113,7 @@ -- | Find cut value. findCut :: HC.Dendrogram a -> HC.Distance-findCut = continuousBy s 9 10 . VU.fromList . F.toList . flattenDist+findCut = quantile s 9 10 . VU.fromList . F.toList . flattenDist where flattenDist (HC.Leaf _) = Seq.empty flattenDist (HC.Branch !d !l !r) =@@ -134,38 +134,17 @@ -- | Hierarchical spectral clustering. hSpecClust :: DenseFlag -> EigenGroup- -> NormType -> Maybe NumEigen -> Maybe Q -> Maybe NumRuns -> SingleCells -> IO (ClusterResults, ClusterGraph CellInfo)-hSpecClust (DenseFlag isDense) eigenGroup norm numEigen minModMay runsMay sc = do+hSpecClust (DenseFlag isDense) eigenGroup numEigen minModMay runsMay sc = do let items = V.zipWith (\x y -> CellInfo x y) (_rowNames sc) (fmap Row . flip V.generate id . V.length . _rowNames $ sc)- hSpecCommand TfIdfNorm False =- hierarchicalSpectralCluster- eigenGroup- True- (fmap unNumEigen numEigen)- Nothing- minModMay- (fmap unNumRuns runsMay)- items- . Left- hSpecCommand BothNorm False =- hierarchicalSpectralCluster- eigenGroup- True- (fmap unNumEigen numEigen)- Nothing- minModMay- (fmap unNumRuns runsMay)- items- . Left- hSpecCommand _ False =+ hSpecCommand False = hierarchicalSpectralCluster eigenGroup False@@ -175,29 +154,7 @@ (fmap unNumRuns runsMay) items . Left- hSpecCommand TfIdfNorm True =- HSD.hierarchicalSpectralCluster- eigenGroup- True- (fmap unNumEigen numEigen)- Nothing- minModMay- (fmap unNumRuns runsMay)- items- . Left- . sparseToHMat- hSpecCommand BothNorm True =- HSD.hierarchicalSpectralCluster- eigenGroup- True- (fmap unNumEigen numEigen)- Nothing- minModMay- (fmap unNumRuns runsMay)- items- . Left- . sparseToHMat- hSpecCommand _ True =+ hSpecCommand True = HSD.hierarchicalSpectralCluster eigenGroup False@@ -209,7 +166,7 @@ . Left . sparseToHMat - tree <- hSpecCommand norm isDense . unMatObsRow . _matrix $ sc+ tree <- hSpecCommand isDense . unMatObsRow . _matrix $ sc let clustering :: [(CellInfo, [Cluster])] clustering =
src/TooManyCells/MakeTree/Load.hs view
@@ -32,7 +32,7 @@ import qualified Data.ByteString.Streaming.Char8 as BS import qualified Data.Csv as CSV import qualified Data.Foldable as F-import qualified Data.Map as Map+import qualified Data.Map.Strict as Map import qualified Data.Sparse.Common as HS import qualified Data.Text as T import qualified Data.Text.Read as T
src/TooManyCells/MakeTree/Plot.hs view
@@ -92,17 +92,6 @@ return $ (unX . fst $ p, unY . snd $ p, l, c) --- -- | Plot clusters on a 2D axis.--- plotClusters :: [(CellInfo, Cluster)] -> Axis B V2 Double--- plotClusters vs = r2Axis &~ do---- forM vs $ \(CellInfo { _projection = (X !x, Y !y)}, Cluster c) ->--- scatterPlot [(x, y)] $ do--- let color = (cycle colours2) !! c--- plotMarker .= circle 1 # fc color # lwO 1 # lc color---- hideGridLines- -- | Plot clusters. plotClustersR :: String -> ProjectionMap -> [(CellInfo, [Cluster])] -> R s () plotClustersR outputPlot pm clusterList = do@@ -119,6 +108,7 @@ xlab("Projection 1") + ylab("Projection 2") + scale_color_discrete(guide = guide_legend(title = "Cluster", ncol = 3)) ++ theme_cowplot() + theme(aspect.ratio = 1) suppressMessages(ggsave(p, file = outputPlot_hs))@@ -153,6 +143,7 @@ scale_color_manual( guide = guide_legend(title = "Label") , values = colorMap ) ++ theme_cowplot() + theme(aspect.ratio = 1) suppressMessages(ggsave(p, file = outputPlot_hs))@@ -169,12 +160,6 @@ dev.off() |] - -- Plot flat hierarchy.- -- [r| pdf(paste0(outputPlot_hs, "_flat_hierarchy.pdf", sep = ""))- -- plot(clustering_hs)- -- dev.off()- -- |]- -- Plot clustering. [r| colors = rainbow(length(unique(clustering_hs$cluster))) names(colors) = unique(clustering_hs$cluster)@@ -196,30 +181,6 @@ dev.off() |] - -- [r| library(tsne)-- -- colors = rainbow(length(unique(clustering_hs$cluster)))- -- names(colors) = unique(clustering_hs$cluster)-- -- tsneMat = tsne(mat_hs, perplexity=50)-- -- pdf(paste0(outputPlot_hs, "_tsne.pdf", sep = ""))-- -- plot(tsneMat- -- , col=clustering_hs$cluster+1- -- , pch=ifelse(clustering_hs$cluster == 0, 8, 1) # Mark noise as star- -- , cex=ifelse(clustering_hs$cluster == 0, 0.5, 0.75) # Decrease size of noise- -- , xlab=NA- -- , ylab=NA- -- )- -- colors = sapply(1:length(clustering_hs$cluster)- -- , function(i) adjustcolor(palette()[(clustering_hs$cluster+1)[i]], alpha.f = clustering_hs$membership_prob[i])- -- )- -- points(tsneMat, col=colors, pch=20)-- -- dev.off()- -- |]- return () -- | Plot ranked modularity.@@ -242,6 +203,7 @@ # geom_vline(xintercept = cutoff_hs, linetype = "dashed", color = "red") + xlab("Ranking") + ylab("Modularity") ++ theme_cowplot() + theme(aspect.ratio = 1) suppressMessages(ggsave(p, file = outputPlot_hs))@@ -286,6 +248,7 @@ xlab("") + ylab("") + scale_fill_gradient2(guide = guide_colorbar(title = "Clumpiness"), midpoint = 0.5, low = muted("blue"), high = muted("red")) ++ theme_cowplot() + theme(axis.text.x = element_text(angle = 315, hjust = 0)) suppressMessages(ggsave(p, file = outputPlot_hs))
src/TooManyCells/MakeTree/Print.hs view
@@ -12,31 +12,42 @@ ( printClusterDiversity , printClusterInfo , printNodeInfo+ , printLabelMap+ , saveFragments ) where -- Remote import BirchBeer.Types import BirchBeer.Utility (getGraphLeaves, getGraphLeavesWithParents) import Control.Monad (join)-import Data.List (genericLength, intercalate)-import Data.Maybe (fromMaybe, mapMaybe, catMaybes)+import Data.List (genericLength, intercalate, foldl')+import Data.Maybe (fromMaybe, mapMaybe, catMaybes, isJust) import Data.Monoid ((<>))+import Data.Streaming.Zlib (WindowBits (..)) import Safe (headMay) import TextShow (showt)+import Turtle hiding (Size) import qualified Control.Lens as L import qualified Data.ByteString.Lazy.Char8 as B import qualified Data.Csv as CSV import qualified Data.Foldable as F import qualified Data.Graph.Inductive as G+import qualified Data.HashMap.Strict as HMap 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.Encoding as T+import qualified Data.Vector as V+import qualified Filesystem.Path as FP+import qualified Turtle.Bytes as TB -- Local import TooManyCells.Diversity.Types+import TooManyCells.File.Types import TooManyCells.MakeTree.Types import TooManyCells.MakeTree.Cluster+import TooManyCells.Matrix.Types -- | Print the diversity of each leaf cluster. printClusterDiversity :: [(Cluster, Diversity, Size)] -> B.ByteString@@ -59,7 +70,7 @@ getSizes :: ([G.Node], a) -> [Int] getSizes = fmap getSize . fst getSize :: G.Node -> Int- getSize = sum . fmap (maybe 0 Seq.length . snd) . getGraphLeaves gr+ getSize = foldl' (+) 0 . fmap (maybe 0 Seq.length . snd) . getGraphLeaves gr getQs :: ([G.Node], a) -> [Double] getQs = mapMaybe getQ . fst getQ :: G.Node -> Maybe Double@@ -107,7 +118,7 @@ (fmap (\a -> getComposition a (ClusterGraph gr) x) lm) (getNodeChildren gr x) getSize :: G.Node -> Int- getSize = sum . fmap (maybe 0 Seq.length . snd) . getGraphLeaves gr+ getSize = foldl' (+) 0 . fmap (maybe 0 Seq.length . snd) . getGraphLeaves gr getQ :: G.Node -> Maybe Double getQ = join . fmap (L.view edgeDistance . snd) . headMay . G.lsuc gr getSignificance :: G.Node -> Maybe Double@@ -141,3 +152,49 @@ ) ) . nodeInfo lm++-- | Print the label map to a string.+printLabelMap :: LabelMap -> B.ByteString+printLabelMap = (<>) "item,label\n" . CSV.encode + . fmap ( L.over L._2 unLabel+ . L.over L._1 unId+ )+ . Map.toAscList+ . unLabelMap++-- | Print the fragments to a file. Will reassign original barcodes from the+-- original file to those in the single cell matrix by splitting the barcode on+-- "BARCODE-LABEL". Mainly for use with --custom-label.+saveFragments :: OutputDirectory+ -> SingleCells+ -> [Either MatrixFileType MatrixFileType]+ -> IO ()+saveFragments (OutputDirectory output') sc files = sh $ do+ let out = (fromText . T.pack $ output') FP.</> "fragments.tsv.gz"+ newBarcodes = HMap.fromList+ . fmap (\ (Cell !x) -> (fst . T.breakOn "-" $ x, x))+ . V.toList+ . L.view rowNames+ $ sc+ assignRow :: T.Text -> Maybe T.Text+ assignRow = fmap ((<> "\n") . T.intercalate "\t")+ . ( (L.ix 3)+ L.%%~ (flip HMap.lookup newBarcodes . fst . T.breakOn "-")+ )+ . T.splitOn "\t"+ getFile (Left (CompressedFragments (FragmentsFile file))) = file+ getFile x = error $ "saveFragments: Not a fragments.tsv.gz file: "+ <> show x+ loadFile =+ TB.decompress (WindowBits 31) . TB.input . fromText . T.pack . getFile++ TB.output out+ . TB.compress 6 (WindowBits 31)+ . fmap (maybe "" T.encodeUtf8)+ . mfilter isJust+ . fmap (assignRow . lineToText)+ . toLines+ . TB.toUTF8+ . msum+ . fmap loadFile+ $ files
src/TooManyCells/MakeTree/Types.hs view
@@ -53,6 +53,13 @@ newtype IsLeaf = IsLeaf {unIsLeaf :: Bool} deriving (Eq, Ord, Read, Show) newtype DenseFlag = DenseFlag { unDenseFlag :: Bool } deriving (Eq, Ord, Read, Show)+newtype LabelMapOutputFlag = LabelMapOutputFlag { unLabelMapOutputFlag :: Bool }+ deriving (Eq, Ord, Read, Show)+newtype FragmentsOutputFlag =+ FragmentsOutputFlag { unFragmentsOutputFlag :: Bool }+ deriving (Eq, Ord, Read, Show)+newtype UpdateTreeRowsFlag = UpdateTreeRowsFlag { unUpdateTreeRowsFlag :: Bool }+ deriving (Read, Show) newtype AdjacencyMat = AdjacencyMat { unAdjacencyMat :: H.Matrix H.R } deriving (Read,Show)@@ -113,11 +120,6 @@ instance TreeItem CellInfo where getId = Id . unCell . _barcode--instance MatrixLike SingleCells where- getMatrix = unMatObsRow . _matrix- getRowNames = fmap unCell . _rowNames- getColNames = fmap unGene . _colNames instance A.ToJSON Q where toEncoding = A.genericToEncoding A.defaultOptions
+ src/TooManyCells/MakeTree/Utility.hs view
@@ -0,0 +1,51 @@+{- TooManyCells.MakeTree.Utility+Gregory W. Schwartz++Collects utility functions for the tree.+-}++{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE OverloadedStrings #-}++module TooManyCells.MakeTree.Utility+ ( updateTreeRowIndex+ , updateTreeRowBool+ ) where++-- Remote+import BirchBeer.Types+import Data.Maybe (fromMaybe)+import Data.Tree (Tree)+import qualified Control.Lens as L+import qualified Data.HashMap.Strict as HMap+import qualified Data.Vector as V++-- Local+import TooManyCells.MakeTree.Types+import TooManyCells.Matrix.Types++-- | Update the CellInfo row index with new matrix row index values.+updateTreeRowIndex :: SingleCells+ -> Tree (TreeNode (V.Vector CellInfo))+ -> Tree (TreeNode (V.Vector CellInfo))+updateTreeRowIndex sc = fmap (L.over item (fmap (V.mapMaybe updateCell)))+ where+ cellRowMap =+ HMap.fromList . flip zip (fmap Row [0..]) . V.toList . getRowNames $ sc+ updateCell c = do+ newRow <- flip HMap.lookup cellRowMap+ . unCell+ . L.view barcode+ $ c+ return $ L.set cellRow newRow c++-- | updateTreeRowIndex meant to be called from the main program with these+-- options.+updateTreeRowBool :: UpdateTreeRowsFlag+ -> Maybe SingleCells+ -> Tree (TreeNode (V.Vector CellInfo))+ -> Tree (TreeNode (V.Vector CellInfo))+updateTreeRowBool (UpdateTreeRowsFlag False) _ tree = tree+updateTreeRowBool (UpdateTreeRowsFlag True) Nothing tree = tree+updateTreeRowBool (UpdateTreeRowsFlag True) (Just sc) tree =+ updateTreeRowIndex sc tree
+ src/TooManyCells/Matrix/AtacSeq.hs view
@@ -0,0 +1,99 @@+{- TooManyCells.MakeTree.AtacSeq+Gregory W. Schwartz++Collects functions pertaining to converting a feature count matrix to an AtacSeq+matrix, where each feature is a bin range.+-}++{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE TupleSections #-}++module TooManyCells.Matrix.AtacSeq+ ( rangeToBin+ , rangeToBinSc+ , binarizeSc+ ) where++-- Remote+import BirchBeer.Types+import Data.Bool (bool)+import Data.Maybe (fromMaybe)+import Data.List (sort, foldl')+import TextShow (showt)+import qualified Control.Lens as L+import qualified Data.Attoparsec.Text as A+import qualified Data.Foldable as F+import qualified Data.IntMap as IMap+import qualified Data.IntervalMap.Interval as I+import qualified Data.Map.Strict as Map+import qualified Data.Sequence as Seq+import qualified Data.Set as Set+import qualified Data.Sparse.Common as S+import qualified Data.Text as T+import qualified Data.Vector as V++-- Local+import TooManyCells.Matrix.Types++-- | Convert an input range to a bin.+rangeToBin :: BinWidth -> ChrRegion -> ChrRegionBin+rangeToBin (BinWidth !b) region@(ChrRegion (!l, _)) =+ ChrRegionBin $ ChrRegion (l, I.ClosedInterval lb (lb + b - 1))+ where+ lb = ((div (midpoint region) b) * b) + 1++-- | Find midpoint of a range.+midpoint :: ChrRegion -> Int+midpoint (ChrRegion (_, !i)) =+ I.lowerBound i + div (I.upperBound i - I.lowerBound i) 2++-- | Get the map of ranges to their bins, keeping bin order.+binsToRangeBinMapWithOrder :: BinWidth+ -> [ChrRegion]+ -> (RangeBinMap, [ChrRegionBin])+binsToRangeBinMapWithOrder b rs =+ (rangeBinMap, Set.toAscList . Set.fromList $ bins)+ where+ rangeBinMap = RangeBinMap . IMap.fromList . zip [0..] . fmap fst $ rangeBins+ rangeBins = fmap (\ !v -> (Map.findWithDefault (error "Bin missing index in binsToRangeBinMapWithOrder") v binMap, v)) bins -- Insert correct idx+ binMap = Map.fromList+ . flip zip (fmap BinIdx [0..])+ . Set.toAscList+ . Set.fromList+ $ bins+ bins = fmap (rangeToBin b) rs++-- | Convert a range matrix to a bin matrix.+rangeToBinMat :: RangeBinMap -> MatObsRow -> MatObsRow+rangeToBinMat (RangeBinMap bm) (MatObsRow mat) =+ MatObsRow . foldl' addToMat init . S.toListSM $ mat+ where+ addToMat !m val = m S.^+^ S.fromListSM (S.dimSM init) [rangeToBinVal val] -- Possible space leak with ^+^+ init = S.zeroSM (S.nrows mat) . Set.size . Set.fromList . IMap.elems $ bm+ rangeToBinVal all@(!i, !j, !v) = (i, unBinIdx $ rangeToBinCol all, v)+ rangeToBinCol all@(_, !j, _) = IMap.findWithDefault (error $ "Missing range index in rangeToBinMat for: " <> show all)+ j+ bm++-- | Convert a range SingleCells matrix to a bin SingleCells matrix.+rangeToBinSc :: BinWidth -> SingleCells -> SingleCells+rangeToBinSc b sc =+ SingleCells { _matrix = rangeToBinMat rangeBinMap . L.view matrix $ sc+ , _rowNames = L.view rowNames sc+ , _colNames =+ V.fromList . fmap (Feature . showt . unChrRegionBin) $ bins+ }+ where+ (rangeBinMap, bins) =+ binsToRangeBinMapWithOrder b+ . fmap ( either (error . (<>) "Cannot parse region format `chrN:START-END` in: ") id+ . parseChrRegion+ . unFeature+ )+ . V.toList+ . L.view colNames+ $ sc++-- | Binarize a matrix.+binarizeSc :: SingleCells -> SingleCells+binarizeSc = L.over matrix (MatObsRow . fmap (bool 0 1 . (> 0)) . unMatObsRow)
src/TooManyCells/Matrix/Load.hs view
@@ -16,42 +16,62 @@ , loadHMatrixData , loadSparseMatrixData , loadSparseMatrixDataStream+ , loadFragments , loadProjectionMap+ , loadBdgBW ) where -- Remote import BirchBeer.Types import Codec.Compression.GZip (decompress)-import Control.DeepSeq (force)+import Control.DeepSeq (force, deepseq) import Control.Exception (evaluate) import Control.Monad.Except (runExceptT, ExceptT (..)) import Control.Monad.Managed (with, liftIO, Managed (..))+import Control.Monad.Trans.Resource (runResourceT, MonadResource) import Data.Char (ord)-import Data.Matrix.MatrixMarket (readMatrix)+import Data.Function (on)+import Data.List (sortBy, sort, foldl')+import Data.Matrix.MatrixMarket (readMatrix, readMatrix') import Data.Maybe (fromMaybe, maybe) import Data.Monoid ((<>))+import Data.Streaming.Zlib (WindowBits (..)) import Data.Vector (Vector)-import Safe+import Safe (atMay, headMay) import System.IO.Temp (withSystemTempFile)+import TextShow (showt)+import qualified Control.Foldl as Fold import qualified Control.Lens as L-import qualified Data.ByteString.Lazy.Char8 as B+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 Data.Csv as CSV import qualified Data.Foldable as F-import qualified Data.Map as Map+import qualified Data.HashMap.Strict as HMap+import qualified Data.HashSet as HSet+import qualified Data.IntMap.Strict as IMap+import qualified Data.IntervalMap.Interval as I+import qualified Data.Map.Strict as Map+import qualified Data.Set as Set import qualified Data.Sparse.Common as HS import qualified Data.Text as T+import qualified Data.Text.Encoding as T import qualified Data.Text.Read as T import qualified Data.Vector as V import qualified Numeric.LinearAlgebra as H import qualified Streaming as S import qualified Streaming.Cassava as S import qualified Streaming.Prelude as S-import qualified Streaming.With.Lifted as SW+import qualified Streaming.With.Lifted as S+import qualified Streaming.Zip as S import qualified System.IO as IO+import qualified Turtle as Turtle+import qualified Turtle.Bytes as TB+import qualified Turtle.Line as Turtle -- Local import TooManyCells.File.Types+import TooManyCells.Matrix.AtacSeq import TooManyCells.Matrix.Types import TooManyCells.Matrix.Utility @@ -64,7 +84,7 @@ -- | Load output of cellranger. loadCellrangerData :: FeatureColumn- -> GeneFile+ -> FeatureFile -> CellFile -> MatrixFileFolder -> IO SingleCells@@ -75,23 +95,23 @@ m <- fmap (MatObsRow . HS.transposeSM . matToSpMat) -- We want observations as rows . readMatrix $ mf- g <- fmap (\ x -> either error (fmap (Gene . getFeature fc)) ( CSV.decodeWith csvOptsTabs CSV.NoHeader x+ g <- fmap (\ x -> either error (fmap (Feature . getFeature fc)) ( CSV.decodeWith csvOptsTabs CSV.NoHeader x :: Either String (Vector [T.Text]) ) )- . B.readFile- . unGeneFile+ . BL.readFile+ . unFeatureFile $ gf c <- fmap (\ x -> either error (fmap (Cell . head)) ( CSV.decodeWith csvOptsTabs CSV.NoHeader x :: Either String (Vector [T.Text]) ) )- . B.readFile+ . BL.readFile . unCellFile $ cf return $- SingleCells { _matrix = m -- We want observations as rows.+ SingleCells { _matrix = MatObsRow . HS.sparsifySM . unMatObsRow $ m -- We want observations as rows. , _rowNames = c , _colNames = g }@@ -99,36 +119,34 @@ -- | Load output of cellranger >= 3.0.0 loadCellrangerDataFeatures :: FeatureColumn- -> GeneFile+ -> FeatureFile -> CellFile -> MatrixFileFolder -> IO SingleCells loadCellrangerDataFeatures _ _ _ (MatrixFolder mf) = error "Expected matrix.mtx.gz, impossible error."-loadCellrangerDataFeatures fc gf cf (MatrixFile mf) = withSystemTempFile "temp_mat.mtx" $ \tempMatFile h -> do+loadCellrangerDataFeatures fc gf cf (MatrixFile mf) = do let csvOptsTabs = CSV.defaultDecodeOptions { CSV.decDelimiter = fromIntegral (ord '\t') } - B.readFile mf >>= B.hPut h . decompress >> IO.hClose h-- m <- fmap (MatObsRow . HS.transposeSM . matToSpMat) -- We want observations as rows- . readMatrix- $ tempMatFile- g <- fmap (\ x -> either error (fmap (Gene . getFeature fc)) ( CSV.decodeWith csvOptsTabs CSV.NoHeader (decompress x)+ m <- BL.readFile mf+ >>= readMatrix' . decompress+ >>= pure . MatObsRow . HS.transposeSM . matToSpMat -- We want observations as rows+ g <- fmap (\ x -> either error (fmap (Feature . getFeature fc)) ( CSV.decodeWith csvOptsTabs CSV.NoHeader (decompress x) :: Either String (Vector [T.Text]) ) )- . B.readFile- . unGeneFile+ . BL.readFile+ . unFeatureFile $ gf c <- fmap (\ x -> either error (fmap (Cell . head)) ( CSV.decodeWith csvOptsTabs CSV.NoHeader (decompress x) :: Either String (Vector [T.Text]) ) )- . B.readFile+ . BL.readFile . unCellFile $ cf return $- SingleCells { _matrix = m -- We want observations as rows.+ SingleCells { _matrix = MatObsRow . HS.sparsifySM . unMatObsRow $ m -- We want observations as rows. , _rowNames = c , _colNames = g }@@ -146,11 +164,11 @@ :: Either String (Vector (Vector T.Text)) ) )- . B.readFile+ . BL.readFile $ mf let c = fmap Cell . V.drop 1 . V.head $ all- g = fmap (Gene . V.head) . V.drop 1 $ all+ g = fmap (Feature . V.head) . V.drop 1 $ all m = MatObsRow . hToSparseMat . H.tr -- We want observations as rows@@ -161,7 +179,7 @@ $ all return $- SingleCells { _matrix = m+ SingleCells { _matrix = MatObsRow . HS.sparsifySM . unMatObsRow $ m , _rowNames = c , _colNames = g }@@ -175,14 +193,14 @@ loadSparseMatrixData (Delimiter delim) (MatrixFile mf) = do let csvOpts = CSV.defaultDecodeOptions { CSV.decDelimiter = fromIntegral (ord delim) }- strictRead path = (evaluate . force) =<< B.readFile path+ strictRead path = (evaluate . force) =<< BL.readFile path all <- fmap (\x -> either error id $ (CSV.decodeWith csvOpts CSV.NoHeader x :: Either String (V.Vector [T.Text]))) . strictRead $ mf let c = V.fromList . fmap Cell . drop 1 . V.head $ all- g = fmap (Gene . head) . V.drop 1 $ all+ g = fmap (Feature . head) . V.drop 1 $ all m = MatObsRow . HS.sparsifySM . HS.fromColsV -- We want observations as rows@@ -191,7 +209,7 @@ $ all return $- SingleCells { _matrix = m+ SingleCells { _matrix = MatObsRow . HS.sparsifySM . unMatObsRow $ m , _rowNames = c , _colNames = g }@@ -207,17 +225,17 @@ { S.decDelimiter = fromIntegral (ord delim) } cS = fmap (S.first (fmap Cell . drop 1 . fromMaybe (error "\nNo header."))) . S.head- gS = S.toList . S.map (Gene . head) . S.drop 1+ gS = S.toList . S.map (Feature . head) . S.drop 1 mS = S.toList . S.map (HS.sparsifySV . HS.vr . fmap (either error fst . T.double) . drop 1) . S.drop 1 res <- flip with return $ do - contents <- SW.withBinaryFileContents mf+ contents <- S.withBinaryFileContents mf (c S.:> g S.:> m S.:> _) <-- fmap (either (error . show) id)+ fmap (either (\x -> error $ show x <> " Expecting csv format. Is this the correct file name and file format?") id) . runExceptT . cS . (S.store gS)@@ -228,13 +246,147 @@ let finalMat = MatObsRow . HS.sparsifySM . HS.fromColsL $ m -- We want observations as rows return $- SingleCells { _matrix = finalMat+ SingleCells { _matrix = MatObsRow . HS.sparsifySM . unMatObsRow $ finalMat , _rowNames = V.fromList c , _colNames = V.fromList g } return res +-- | Load a range feature list streaming in TSV 10x fragments.tsv.gz format.+loadFragments :: Maybe CellWhitelist+ -> Maybe BlacklistRegions+ -> Maybe ExcludeFragments+ -> Maybe BinWidth+ -> FragmentsFile+ -> IO SingleCells+loadFragments whitelist blacklist excludeFragments binWidth (FragmentsFile mf) = do+ let readDecimal = either error fst . T.decimal+ readDouble = either error fst . T.double+ labelsFold = (,) <$> cellsFold <*> featuresFold+ cellsFold = Fold.premap (\(!c, _, _) -> c) hashNub+ featuresFold = Fold.premap (\(_, !r, _) -> r) hashNub+ preprocessStream = maybe id filterCells whitelist+ . Turtle.mfilter (\(_, _, !x) -> x /= 0) -- Ignore missing regions+ . fmap parseLine+ filterCells (CellWhitelist wl) =+ Turtle.mfilter (\(!b, _, _) -> HSet.member b wl)+ filterFragments (ExcludeFragments ef) = Turtle.mfilter (T.isInfixOf ef)+ filterBlacklist (BlacklistRegions br) =+ Turtle.inproc "bedtools" ["subtract", "-A", "-a", "stdin", "-b", br]+ parseLine (chrom:start:end:barcode:duplicateCount:_) =+ ( barcode+ , maybe showt (\x -> showt . unChrRegionBin . rangeToBin x) binWidth+ . either (\x -> error $ "Cannot read region in file: " <> x) id+ . parseChrRegion+ $ mconcat [chrom, ":", start, "-", end]+ , readDouble duplicateCount+ )+ parseLine xs = error $ "loadFragments: Unexpected number of columns, did you use the fragments.tsv.gz 10x format? Input: " <> show xs+ 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)+ . fmap Turtle.lineToText+ . (maybe id filterBlacklist blacklist)+ . Turtle.inproc "bedtools" ["sort", "-i", "stdin"]+ . Turtle.toLines+ . TB.toUTF8+ . TB.decompress (WindowBits 31)+ . TB.input+ . Turtle.fromText+ . T.pack+ $ mf++ (cellsList, featuresList) <- Turtle.reduce labelsFold stream++ let cells = V.fromList cellsList+ features = V.fromList featuresList+ cellMap = HMap.fromList . flip zip ([0..] :: [Int]) $ cellsList+ featureMap = HMap.fromList . flip zip ([0..] :: [Int]) $ featuresList+ findErr x = fromMaybe (error $ "loadFragments: No indices found: " <> show x)+ . HMap.lookup x+ matFold = Fold.Fold addToMat init MatObsRow+ addToMat !m val = m HS.^+^ HS.fromListSM (HS.dimSM init) [val]+ init = HS.zeroSM (V.length cells) (V.length features)+ getIndices (!c, !r, !v) =+ (findErr c cellMap, findErr r featureMap, v)++ -- Get map second.+ mat <- Turtle.reduce matFold+ . fmap getIndices+ $ stream++ return $+ SingleCells { _matrix = MatObsRow . HS.sparsifySM . unMatObsRow $ mat+ , _rowNames = fmap Cell cells+ , _colNames = fmap Feature features+ }++-- | Load a bedGraph or bigWig file as a reference to classify cells.+loadBdgBW :: Maybe BlacklistRegions+ -> Maybe ExcludeFragments+ -> Maybe BinWidth+ -> MatrixFileType+ -> IO SingleCells+loadBdgBW blacklist excludeFragments binWidth inFile = do+ let readDecimal = either error fst . T.decimal+ readDouble = either error fst . T.double+ featuresFold = Fold.premap (\(!r, _) -> r) hashNub+ preprocessStream = Turtle.mfilter (\(_, !x) -> x /= 0) -- Ignore missing regions+ . fmap parseLine+ filterFragments (ExcludeFragments ef) = Turtle.mfilter (T.isInfixOf ef)+ filterBlacklist (BlacklistRegions br) =+ Turtle.inproc "bedtools" ["subtract", "-A", "-a", "stdin", "-b", br]+ parseLine (chrom:start:end:duplicateCount:_) =+ ( maybe showt (\x -> showt . unChrRegionBin . rangeToBin x) binWidth+ . either (\x -> error $ "Cannot read region in file: " <> x) id+ . parseChrRegion+ $ mconcat [chrom, ":", start, "-", end]+ , readDouble duplicateCount+ )+ parseLine xs = error $ "loadBW: Unexpected number of columns, did you use the bigWig format? Input: " <> show xs+ fromFile (BigWig file) =+ Turtle.inproc "bigWigToBedGraph" [T.pack file, "stdout"] mempty+ fromFile (BedGraph file) = Turtle.input . Turtle.fromText . T.pack $ file+ fromFile file =+ error $ "Only accepts BedGraph or BigWig for now: " <> show file+ getFileName (BigWig file) = file+ getFileName (BedGraph file) = file+ getFileName file =+ error $ "Only accepts BedGraph or BigWig for now: " <> show file+ 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)+ . fmap Turtle.lineToText+ . (maybe id filterBlacklist blacklist)+ . Turtle.inproc "bedtools" ["sort", "-i", "stdin"]+ $ fromFile inFile++ featuresList <- Turtle.reduce featuresFold stream++ let features = V.fromList featuresList+ featureMap = HMap.fromList . flip zip ([0..] :: [Int]) $ featuresList+ findErr x = fromMaybe (error $ "loadBW: No indices found: " <> show x)+ . HMap.lookup x+ matFold = Fold.Fold addToMat init MatObsRow+ addToMat !m (!i, !j, !x) = HS.insertSpMatrix i j x m+ init = HS.zeroSM 1 (V.length features)+ getIndices (!r, !v) =+ (0, findErr r featureMap, v)++ -- Get map second.+ mat <- Turtle.reduce matFold+ . fmap getIndices+ $ stream++ return $+ SingleCells { _matrix = MatObsRow . HS.sparsifySM . unMatObsRow $ mat+ , _rowNames = V.singleton . Cell . T.pack $ getFileName inFile+ , _colNames = fmap Feature features+ }+ -- | Load a projection file to get a vector of projections. loadProjectionFile :: ProjectionFile -> IO (Vector (X, Y)) loadProjectionFile =@@ -243,7 +395,7 @@ (fmap getProjection . V.drop 1) (CSV.decode CSV.NoHeader x :: Either String (Vector [T.Text])) )- . B.readFile+ . BL.readFile . unProjectionFile where getProjection [_, x, y] = ( X . either error fst . T.double $ x@@ -263,7 +415,7 @@ (fmap getProjection . V.drop 1) $ (CSV.decode CSV.NoHeader x :: Either String (Vector [T.Text])) )- . B.readFile+ . BL.readFile . unProjectionFile where getProjection [b, x, y] = ( Cell b
src/TooManyCells/Matrix/Preprocess.hs view
@@ -7,14 +7,19 @@ {-# LANGUAGE BangPatterns #-} {-# LANGUAGE QuasiQuotes #-} {-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TupleSections #-} module TooManyCells.Matrix.Preprocess ( scaleRMat , scaleDenseMat , scaleSparseMat+ , totalScaleSparseMat , logCPMSparseMat , uqScaleSparseMat , medScaleSparseMat+ , quantileScaleSparseMat+ , centerScaleSparseCell+ , tfidfScaleSparseMat , filterRMat , filterDenseMat , filterNumSparseMat@@ -24,20 +29,45 @@ , removeCorrelated , pcaRMat , pcaDenseSc+ , pcaSparseSc+ , lsaSparseSc+ , svdSparseSc , shiftPositiveSc+ , pcaDenseMat+ , pcaSparseMat+ , lsaSparseMat+ , svdSparseMat+ , shiftPositiveMat+ , emptyMatErr+ , labelRows+ , labelCols+ , fastBinJoinRows+ , fastBinJoinCols+ , transformChrRegions ) where -- Remote+import BirchBeer.Types (LabelMap (..), Id (..), Label (..), Feature (..)) import Data.Bool (bool)-import Data.List (sort)+import Data.List (sort, foldl', transpose) import Data.Monoid ((<>))-import Data.Maybe (fromMaybe)+import Data.Maybe (fromMaybe, isJust) import H.Prelude (io) import Language.R as R import Language.R.QQ (r)-import Statistics.Quantile (continuousBy, s)+import Math.Clustering.Spectral.Sparse (B1 (..), B2 (..), b1ToB2)+import Statistics.Quantile (quantile, s)+import Statistics.Sample (mean, stdDev) import TextShow (showt)+import qualified Control.Foldl as Fold import qualified Control.Lens as L+import qualified Data.HashMap.Strict as HMap+import qualified Data.HashSet as HSet+import qualified Data.IntMap as IMap+import qualified Data.IntSet as ISet+import qualified Data.IntervalMap.Interval as IntervalMap+import qualified Data.IntervalMap.Strict as IntervalMap+import qualified Data.Map.Strict as Map import qualified Data.Set as Set import qualified Data.Sparse.Common as S import qualified Data.Text as T@@ -45,14 +75,25 @@ import qualified Data.Vector as V import qualified Data.Vector.Algorithms.Radix as V import qualified Data.Vector.Storable as VS+import qualified Math.Clustering.Spectral.Sparse as S import qualified Numeric.LinearAlgebra as H+import qualified Numeric.LinearAlgebra.Devel as H import qualified Numeric.LinearAlgebra.HMatrix as H+import qualified Numeric.LinearAlgebra.SVD.SVDLIBC as SVD -- Local import TooManyCells.File.Types import TooManyCells.Matrix.Types import TooManyCells.Matrix.Utility +-- | Empty matrix error.+emptyMatErr :: String -> String+emptyMatErr checkType = "Matrix is empty in " <> checkType <> ". Check --filter-thresholds, --normalization, or the input matrix for over filtering or incorrect input format."++-- | Check if valid features or cells are empty, error if so.+emptyMatCheckErr :: String -> [a] -> [a]+emptyMatCheckErr checkType xs = bool xs (error $ emptyMatErr checkType) . null $ xs+ -- | Scale a matrix. scaleRMat :: RMatObsRow s -> R s (RMatObsRow s) scaleRMat (RMatObsRow mat) = do@@ -62,7 +103,8 @@ t(mat[,colSums(!is.na(mat)) > 0]) |] --- | Scale a matrix based on the library size.+-- | Scale a matrix based on the library size for each cell and median feature+-- size. scaleDenseMat :: MatObsRow -> MatObsRow scaleDenseMat (MatObsRow mat) = MatObsRow . hToSparseMat@@ -75,7 +117,8 @@ . sparseToHMat $ mat --- | Scale a matrix based on the library size.+-- | Scale a matrix based on the library size for each cell and median feature+-- size. scaleSparseMat :: MatObsRow -> MatObsRow scaleSparseMat (MatObsRow mat) = MatObsRow . S.sparsifySM@@ -87,36 +130,78 @@ . S.toRowsL $ mat +-- | Scale a matrix based on the library size.+totalScaleSparseMat :: MatObsRow -> MatObsRow+totalScaleSparseMat = MatObsRow+ . S.sparsifySM+ . S.fromRowsL+ . fmap scaleSparseCell+ . S.toRowsL+ . unMatObsRow+ -- | Scale a matrix based on the upper quartile. uqScaleSparseMat :: MatObsRow -> MatObsRow uqScaleSparseMat (MatObsRow mat) = MatObsRow . S.sparsifySM- . S.transposeSM- . S.fromColsL+ . S.fromRowsL . fmap uqScaleSparseCell . S.toRowsL $ mat -- | Scale a matrix based on log(CPM + 1).-logCPMSparseMat :: MatObsRow -> MatObsRow-logCPMSparseMat (MatObsRow mat) = MatObsRow- . S.sparsifySM- . S.transposeSM- . S.fromColsL- . fmap logCPMSparseCell- . S.toRowsL- $ mat+logCPMSparseMat :: Double -> MatObsRow -> MatObsRow+logCPMSparseMat b (MatObsRow mat) = MatObsRow+ . S.sparsifySM+ . S.fromRowsL+ . fmap (logCPMSparseCell b)+ . S.toRowsL+ $ mat -- | Scale a matrix based on the median. medScaleSparseMat :: MatObsRow -> MatObsRow medScaleSparseMat (MatObsRow mat) = MatObsRow . S.sparsifySM- . S.transposeSM- . S.fromColsL+ . S.fromRowsL . fmap medScaleSparseCell . S.toRowsL $ mat +-- | Scale a matrix based quantile normalization. Will ignore 0s if missing from+-- input.+quantileScaleSparseMat :: MatObsRow -> MatObsRow+quantileScaleSparseMat (MatObsRow mat) =+ MatObsRow+ . S.sparsifySM+ . fmap (flip (IMap.findWithDefault 0) totalRankMap)+ $ rankMat+ where+ totalRankMap = IMap.fromList+ . zip [1,2..]+ . fmap (Fold.fold avg)+ . transpose+ . fmap (sort . fmap snd . S.toListSV)+ . S.toRowsL+ $ mat+ avg = (/) <$> Fold.sum <*> Fold.genericLength+ rankMat :: S.SpMatrix Int+ rankMat = S.fromRowsL . fmap rankTransform . S.toRowsL $ mat+ rankTransform :: S.SpVector Double -> S.SpVector Int+ rankTransform xs = fmap (flip (Map.findWithDefault 0) rankMap) xs+ where+ rankMap :: Map.Map Double Int+ rankMap = Map.fromList+ . flip zip [1,2..]+ . Set.toList+ . Set.fromList+ . sort+ . fmap snd+ . S.toListSV+ $ xs++-- | TF-IDF normalization.+tfidfScaleSparseMat :: MatObsRow -> MatObsRow+tfidfScaleSparseMat = MatObsRow . unB2 . b1ToB2 . B1 . unMatObsRow+ -- | Scale a cell by the library size. scaleDenseCell :: H.Vector H.R -> H.Vector H.R scaleDenseCell xs = H.cmap (/ total) xs@@ -127,20 +212,20 @@ scaleSparseCell :: S.SpVector Double -> S.SpVector Double scaleSparseCell xs = fmap (/ total) xs where- total = sum xs+ total = foldl' (+) 0 xs -- | log(CPM + 1) normalization for cell.-logCPMSparseCell :: S.SpVector Double -> S.SpVector Double-logCPMSparseCell xs = fmap (logBase 2 . (+ 1) . cpm) xs+logCPMSparseCell :: Double -> S.SpVector Double -> S.SpVector Double+logCPMSparseCell b xs = fmap (logBase b . (+ 1) . cpm) xs where cpm x = x / tpm- tpm = sum xs / 1000000+ tpm = foldl' (+) 0 xs / 1000000 -- | Upper quartile scale cells. uqScaleSparseCell :: S.SpVector Double -> S.SpVector Double uqScaleSparseCell xs = fmap (/ uq) xs where- uq = continuousBy s 3 4+ uq = quantile s 3 4 . VS.fromList . filter (/= 0) . fmap snd@@ -151,31 +236,49 @@ medScaleSparseCell :: S.SpVector Double -> S.SpVector Double medScaleSparseCell xs = fmap (/ med) xs where- med = continuousBy s 2 4+ med = quantile s 2 4 . VS.fromList . filter (/= 0) . fmap snd . S.toListSV $ xs +-- | Center and scale a cell.+centerScaleSparseCell :: S.SpVector Double -> S.SpVector Double+centerScaleSparseCell xs = fmap standard xs+ where+ standard x+ | sigma == 0 = 0+ | otherwise = (x - mu) / sigma+ mu = mean v+ sigma = stdDev v+ v = S.toVector xs++-- | Center a sparse vector.+centerSparseVector :: S.SpVector Double -> S.SpVector Double+centerSparseVector xs = fmap center xs+ where+ center x = x - mu+ mu = mean $ S.toVector xs+ -- | Median scale molecules across cells. scaleDenseMol :: H.Vector H.R -> H.Vector H.R scaleDenseMol xs = H.cmap (/ med) xs where- med = continuousBy s 2 4 . VS.filter (/= 0) $ xs+ med = quantile s 2 4 . VS.filter (/= 0) $ xs -- | Median scale molecules across cells. scaleSparseMol :: S.SpVector Double -> S.SpVector Double scaleSparseMol xs = fmap (/ med) xs where- med = continuousBy s 2 4+ med = quantile s 2 4 . VS.fromList . filter (/= 0) . fmap snd . S.toListSV $ xs --- | Filter a matrix to remove low count cells and genes.+-- | Filter a matrix to remove low count cells and features. filterDenseMat :: FilterThresholds -> SingleCells -> SingleCells filterDenseMat (FilterThresholds (rowThresh, colThresh)) sc = SingleCells { _matrix = m@@ -187,30 +290,32 @@ rowFilter = (>= rowThresh) . H.sumElements colFilter = (>= colThresh) . H.sumElements mat = sparseToHMat . unMatObsRow . _matrix $ sc- validRows = Set.fromList+ validRows = ISet.fromList+ . emptyMatCheckErr "cells" . fmap fst . filter (rowFilter . snd) . zip [0..] . H.toRows $ mat- validCols = Set.fromList+ validCols = ISet.fromList+ . emptyMatCheckErr "features" . fmap fst . filter (colFilter . snd) . zip [0..] . H.toColumns $ mat filteredMat = mat- H.?? ( H.Pos (H.idxs (Set.toAscList validRows))- , H.Pos (H.idxs (Set.toAscList validCols))+ H.?? ( H.Pos (H.idxs (ISet.toAscList validRows))+ , H.Pos (H.idxs (ISet.toAscList validCols)) )- r = V.ifilter (\i _ -> Set.member i validRows)+ r = V.ifilter (\i _ -> ISet.member i validRows) . _rowNames $ sc- c = V.ifilter (\i _ -> Set.member i validCols)+ c = V.ifilter (\i _ -> ISet.member i validCols) . _colNames $ sc --- | Filter a matrix to remove low count cells and genes.+-- | Filter a matrix to remove low count cells and features. filterNumSparseMat :: FilterThresholds -> SingleCells -> SingleCells filterNumSparseMat (FilterThresholds (rowThresh, colThresh)) sc = SingleCells { _matrix = m@@ -219,41 +324,41 @@ } where m = MatObsRow colFilteredMat- rowFilter = (>= rowThresh) . sum- colFilter = (>= colThresh) . sum+ rowFilter = (>= rowThresh) . foldl' (+) 0+ colFilter = (>= colThresh) . foldl' (+) 0 mat = unMatObsRow . _matrix $ sc- mat' = S.transposeSM mat- validRows = Set.fromList+ validRows = ISet.fromList+ . emptyMatCheckErr "cells" . fmap fst . filter (rowFilter . snd) . zip [0..] . S.toRowsL $ mat- validCols = Set.fromList+ validCols = ISet.fromList+ . emptyMatCheckErr "features" . fmap fst . filter (colFilter . snd) . zip [0..] . S.toRowsL . S.transposeSM -- toRowsL is much faster. $ mat- rowFilteredMat = S.transposeSM- . S.fromColsL -- fromRowsL still broken.+ rowFilteredMat = S.fromRowsL -- fromRowsL fixed. . fmap snd- . filter (flip Set.member validRows . fst)+ . filter (flip ISet.member validRows . fst) . zip [0..] . S.toRowsL $ mat colFilteredMat = S.fromColsL . fmap snd- . filter (flip Set.member validCols . fst)+ . filter (flip ISet.member validCols . fst) . zip [0..] . S.toRowsL -- Rows of transpose are faster. . S.transposeSM $ rowFilteredMat- r = V.ifilter (\i _ -> Set.member i validRows)+ r = V.ifilter (\i _ -> ISet.member i validRows) . _rowNames $ sc- c = V.ifilter (\i _ -> Set.member i validCols)+ c = V.ifilter (\i _ -> ISet.member i validCols) . _colNames $ sc @@ -270,13 +375,13 @@ mat = unMatObsRow . _matrix $ sc validIdx = V.modify V.sort . V.map fst- . V.filter (\(_, !c) -> Set.member c wl)+ . V.filter (\(_, (Cell !c)) -> HSet.member c wl) . V.imap (\i v -> (i, v)) . _rowNames $ sc- rowFilteredMat = S.transposeSM- . S.fromColsL -- fromRowsL still broken.+ rowFilteredMat = S.fromRowsL -- fromRowsL fixed. . fmap (S.extractRow mat)+ . emptyMatCheckErr "cells" . V.toList $ validIdx r = V.map (\x -> fromMaybe (error $ "\nWhitelist row index out of bounds (do the whitelist barcodes match the data?): " <> show x <> " out of " <> (show . length . _rowNames $ sc))@@ -291,8 +396,7 @@ contents <- T.readFile file let whiteList = CellWhitelist- . Set.fromList- . fmap Cell+ . HSet.fromList . filter (not . T.null) . T.lines $ contents@@ -332,9 +436,69 @@ [r| mat = prcomp(t(mat_hs), tol = 0.95)$x |] --- | Conduct PCA on a matrix, taking the first principal components.-pcaDenseMat :: PCADim -> MatObsRow -> MatObsRow-pcaDenseMat (PCADim pcaDim) (MatObsRow matObs) =+-- | Conduct PCA on a SingleCells, taking the first principal components.+pcaDenseSc :: DropDimensionFlag -> PCADim -> SingleCells -> SingleCells+pcaDenseSc df p@(PCADim n) sc =+ L.set colNames (V.fromList . fmap (\x -> Feature $ "PCA " <> showt x) $ [begin..end])+ . L.over matrix (pcaDenseMat df p)+ $ sc+ where+ begin = bool 1 (n + 1) . unDropDimensionFlag $ df+ end = bool n ((S.ncols . unMatObsRow . L.view matrix $ sc) - n + 1)+ . unDropDimensionFlag+ $ df++-- | Conduct PCA on a SingleCells, taking the first principal components.+pcaSparseSc :: DropDimensionFlag -> PCADim -> SingleCells -> SingleCells+pcaSparseSc df p@(PCADim n) sc =+ L.set colNames (V.fromList . fmap (\x -> Feature $ "PCA " <> showt x) $ [begin..end])+ . L.over matrix (pcaSparseMat df p)+ $ sc+ where+ begin = bool 1 (n + 1) . unDropDimensionFlag $ df+ end = bool n ((S.ncols . unMatObsRow . L.view matrix $ sc) - n + 1)+ . unDropDimensionFlag+ $ df++-- | Conduct LSA on a SingleCells, taking the first singular values.+lsaSparseSc :: DropDimensionFlag -> LSADim -> SingleCells -> SingleCells+lsaSparseSc df l@(LSADim n) sc =+ L.set colNames (V.fromList . fmap (\x -> Feature $ "LSA " <> showt x) $ [begin..end])+ . L.over matrix (lsaSparseMat df l)+ $ sc+ where+ begin = bool 1 (n + 1) . unDropDimensionFlag $ df+ end = bool n ((S.ncols . unMatObsRow . L.view matrix $ sc) - n + 1)+ . unDropDimensionFlag+ $ df++-- | Conduct SVD on a SingleCells, taking the first singular values.+svdSparseSc :: DropDimensionFlag -> SVDDim -> SingleCells -> SingleCells+svdSparseSc df s@(SVDDim n) sc =+ L.set colNames (V.fromList . fmap (\x -> Feature $ "SVD " <> showt x) $ [begin..end])+ . L.over matrix (svdSparseMat df s)+ $ sc+ where+ begin = bool 1 (n + 1) . unDropDimensionFlag $ df+ end = bool n ((S.ncols . unMatObsRow . L.view matrix $ sc) - n + 1)+ . unDropDimensionFlag+ $ df++-- | Obtain the right singular vectors from N to E on of a sparse+-- matrix.+sparseRightSVD :: Int -> Int -> S.SpMatrix Double -> [S.SpVector Double]+sparseRightSVD n e =+ fmap (S.sparsifySV . S.fromListDenseSV e . drop (n - 1) . H.toList)+ . H.toColumns+ . (\(_, _, !z) -> z)+ . SVD.sparseSvd (e + (n - 1))+ . H.mkCSR+ . fmap (\(!i, !j, !x) -> ((i, j), x))+ . S.toListSM++-- | Conduct PCA on a matrix, retaining the specified number of dimensions.+pcaDenseMat :: DropDimensionFlag -> PCADim -> MatObsRow -> MatObsRow+pcaDenseMat df (PCADim pcaDim) (MatObsRow matObs) = MatObsRow . hToSparseMat . H.takeColumns pcaDim@@ -347,13 +511,57 @@ $ mat where mat = sparseToHMat matObs+ begin = bool 1 (pcaDim + 1) . unDropDimensionFlag $ df+ end = bool pcaDim (S.ncols matObs - pcaDim)+ . unDropDimensionFlag+ $ df --- | Conduct PCA on a SingleCells, taking the first principal components.-pcaDenseSc :: PCADim -> SingleCells -> SingleCells-pcaDenseSc p@(PCADim n) =- L.set colNames (V.fromList . fmap (\x -> Gene $ "PCA " <> showt x) $ [1..n])- . L.over matrix (pcaDenseMat p)+-- | Conduct SVD on a matrix, retaining the specified number of dimensions.+pcaSparseMat :: DropDimensionFlag -> PCADim -> MatObsRow -> MatObsRow+pcaSparseMat df (PCADim pcaDim) (MatObsRow mat) =+ MatObsRow+ . (S.#~#) scaled+ . S.fromRowsL+ . sparseRightSVD begin end+ $ scaled+ where+ scaled = S.fromColsL+ . fmap centerScaleSparseCell+ . S.toRowsL -- Scale features for PCA+ . S.transpose+ $ mat+ begin = bool 1 (pcaDim + 1) . unDropDimensionFlag $ df+ end = bool pcaDim (S.ncols mat - pcaDim) . unDropDimensionFlag $ df +-- | Conduct LSA on a matrix, retaining the specified number of dimensions.+lsaSparseMat :: DropDimensionFlag -> LSADim -> MatObsRow -> MatObsRow+lsaSparseMat df (LSADim lsaDim) (MatObsRow mat) =+ MatObsRow+ . S.fromRowsL+ . S.secondLeft begin end+ . unB2+ . b1ToB2+ . B1+ $ mat+ where+ begin = bool 1 (lsaDim + 1) . unDropDimensionFlag $ df+ end = bool lsaDim (S.ncols mat - lsaDim) . unDropDimensionFlag $ df++-- | Conduct SVD on a matrix, retaining the specified number of dimensions.+svdSparseMat :: DropDimensionFlag -> SVDDim -> MatObsRow -> MatObsRow+svdSparseMat df (SVDDim svdDim) (MatObsRow mat) =+ MatObsRow+ . S.fromRowsL+ . S.secondLeft begin end+ . S.fromRowsL+ . fmap centerScaleSparseCell+ . S.toRowsL+ $ mat+ where+ begin = bool 1 (svdDim + 1) . unDropDimensionFlag $ df+ end = bool svdDim (S.ncols mat - svdDim) . unDropDimensionFlag $ df++ -- | Shift features to positive values. shiftPositiveMat :: MatObsRow -> MatObsRow shiftPositiveMat = MatObsRow@@ -363,7 +571,8 @@ . S.toDenseListSV $ xs )- . S.toColsL+ . S.toRowsL -- toRowsL much faster than toColsL+ . S.transpose . unMatObsRow where shift xs = S.sparsifySV . S.vr . fmap (+ (abs $ minimum xs)) $ xs@@ -371,3 +580,90 @@ -- | Shift features to positive values for SingleCells. shiftPositiveSc :: SingleCells -> SingleCells shiftPositiveSc = L.over matrix shiftPositiveMat++-- | Optionally give the filepath as a label to the rows.+labelRows :: Maybe CustomLabel -> SingleCells -> LabelMat+labelRows Nothing sc = LabelMat (sc, Nothing)+labelRows (Just (CustomLabel l)) sc =+ LabelMat (L.set rowNames newRowNames sc, Just labelMap)+ where+ labelMap = LabelMap+ . Map.fromList+ . flip zip (repeat (Label l))+ . fmap (Id . unCell)+ . V.toList+ $ newRowNames+ newRowNames = fmap labelRow . L.view rowNames $ sc+ labelRow (Cell x) = Cell $ x <> "-" <> l++-- | Optionally give the filepath as a label to the columns when transposing+-- matrix.+labelCols :: Maybe CustomLabel -> SingleCells -> LabelMat+labelCols Nothing sc = LabelMat (sc, Nothing)+labelCols (Just (CustomLabel l)) sc =+ LabelMat (L.set colNames newColNames sc, Just labelMap)+ where+ labelMap = LabelMap+ . Map.fromList+ . flip zip (repeat (Label l))+ . fmap (Id . unFeature)+ . V.toList+ $ newColNames+ newColNames = fmap labelCol . L.view colNames $ sc+ labelCol (Feature x) = Feature $ x <> "-" <> l++-- | Add or remove a single character to speed up joining if region data is+-- binned across rows before transposition.+fastBinJoinRows :: Maybe BinWidth -> Bool -> SingleCells -> SingleCells+fastBinJoinRows binWidth joining =+ bool id (L.over rowNames (fmap (Cell . j . unCell))) $ isJust binWidth+ where+ j = if joining then T.cons 'B' else T.drop 1++-- | Add or remove a single character to speed up joining if region data is+-- binned across columns.+fastBinJoinCols :: Maybe BinWidth -> Bool -> SingleCells -> SingleCells+fastBinJoinCols binWidth joining =+ bool id (L.over colNames (fmap (Feature . j . unFeature)))+ $ isJust binWidth+ where+ j = if joining then T.cons 'B' else T.drop 1++-- | Transform a list of chromosome region features to a list of custom features.+transformChrRegions :: CustomRegions -> SingleCells -> SingleCells+transformChrRegions (CustomRegions customRegions) sc =+ L.set matrix mat . L.set colNames features $ sc+ where+ features = V.fromList . fmap (Feature . showt) $ customRegions+ knownFeatureMap = foldl' (HMap.unionWith IntervalMap.union) HMap.empty+ . fmap (\ (!i, (!c, !interval))+ -> HMap.singleton c+ $ IntervalMap.singleton interval i+ )+ . zip ([0..] :: [Int])+ . fmap unChrRegion+ $ customRegions+ findErrKnown :: T.Text -> [Int]+ findErrKnown x = fromMaybe []+ . (\ (!c, !i) -> HMap.lookup c knownFeatureMap+ >>= Just+ . IntervalMap.elems+ . flip IntervalMap.intersecting i+ )+ . either error unChrRegion+ . parseChrRegion+ $ x+ mat = MatObsRow+ . foldl' addToMat init+ . concatMap (\ (!i, !j, !v)+ -> fmap (i, , v)+ . findErrKnown+ . unFeature+ . fromMaybe (error $ "transformChrRegions/mat: No indices found: " <> show j)+ $ (V.!?) (L.view colNames sc) j+ )+ . S.toListSM+ $ oldMat+ init = S.zeroSM (S.nrows oldMat) $ V.length features+ addToMat !m x = m S.^+^ S.fromListSM (S.dim init) [x]+ oldMat = unMatObsRow . L.view matrix $ sc
src/TooManyCells/Matrix/Types.hs view
@@ -4,59 +4,99 @@ Collects the types used in the program concerning the matrix. -} -{-# LANGUAGE StrictData #-}-{-# LANGUAGE DeriveAnyClass #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DerivingStrategies #-} {-# LANGUAGE DeriveGeneric #-}-{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE DeriveAnyClass #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE QuasiQuotes #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE StrictData #-} {-# LANGUAGE TemplateHaskell #-} module TooManyCells.Matrix.Types where -- Remote-import Control.DeepSeq (NFData (..))+import BirchBeer.Types+import Control.DeepSeq (NFData (..), force)+import Control.Monad (join)+import Control.Parallel (par)+import Control.Parallel.Strategies (withStrategy, rdeepseq, Strategy (..))+import Data.Either (isRight)+import Data.Maybe (fromMaybe) import Data.Monoid (Monoid (..), mempty) import Data.Colour.Palette.BrewerSet (Kolor) import Data.Colour.SRGB (Colour (..), RGB (..), toSRGB)+import Data.Function (on)+import Data.List (sortBy) import Data.Map.Strict (Map) import Data.Text (Text)-import Data.Vector (Vector)+import Data.Vector (Vector (..)) import GHC.Generics (Generic) import Language.R as R import Language.R.QQ (r) import Math.Clustering.Hierarchical.Spectral.Sparse (ShowB) import Math.Clustering.Hierarchical.Spectral.Types (ClusteringTree, ClusteringVertex) import Math.Modularity.Types (Q (..))+import Safe (headMay)+import TextShow (TextShow (..))+import qualified Control.DeepSeq as DS import qualified Control.Lens as L import qualified Data.Aeson as A+import qualified Data.Attoparsec.Text as Atto+import qualified Data.ByteString as B import qualified Data.Clustering.Hierarchical as HC import qualified Data.Graph.Inductive as G-import qualified Data.Map as Map+import qualified Data.HashMap.Strict as HMap+import qualified Data.HashSet as HSet+import qualified Data.IntMap as IMap+import qualified Data.IntervalMap.Interval as IntervalMap+import qualified Data.IntervalMap.Strict as IntervalMap+import qualified Data.Map.Strict as Map import qualified Data.Sequence as Seq import qualified Data.Set as Set import qualified Data.Sparse.Common as S+import qualified Data.Sparse.SpMatrix as S+import qualified Data.Text as T import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U import qualified Numeric.LinearAlgebra as H+import qualified TextShow as TS -- Local -- Basic newtype Cell = Cell { unCell :: Text- } deriving (Eq,Ord,Read,Show,Generic,A.ToJSON, A.FromJSON)+ } deriving (Eq,Ord,Read,Show,Generic) deriving anyclass (A.ToJSON, A.FromJSON) newtype Cols = Cols { unCols :: [Double] }-newtype Gene = Gene { unGene :: Text } deriving (Eq, Ord, Read, Show, Generic)-instance NFData Gene newtype FeatureColumn = FeatureColumn { unFeatureColumn :: Int }+newtype CustomLabel = CustomLabel { unCustomLabel :: Text}+ deriving (Eq,Ord,Read,Show)+newtype BinWidth = BinWidth { unBinWidth :: Int}+newtype NoBinarizeFlag = NoBinarizeFlag { unNoBinarizeFlag :: Bool }+newtype BinarizeFlag = BinarizeFlag { unBinarizeFlag :: Bool }+newtype TransposeFlag = TransposeFlag { unTransposeFlag :: Bool }+newtype DropDimensionFlag = DropDimensionFlag { unDropDimensionFlag :: Bool }+newtype BinIdx = BinIdx { unBinIdx :: Int} deriving (Eq, Ord, Read, Show, Generic) deriving newtype (Num)+newtype RangeIdx = RangeIdx { unRangeIdx :: Int} deriving (Eq, Ord, Read, Show, Generic) deriving newtype (Num) newtype CellWhitelist = CellWhitelist- { unCellWhitelist :: Set.Set Cell- } deriving (Eq,Ord,Read,Show)+ { unCellWhitelist :: HSet.HashSet Text+ } deriving (Eq,Ord,Show)+newtype ExcludeFragments = ExcludeFragments { unExcludeFragments :: T.Text }+newtype BlacklistRegions = BlacklistRegions { unBlacklistRegions :: T.Text } newtype PCADim = PCADim { unPCADim :: Int } deriving (Eq,Ord,Read,Show)-newtype NoFilterFlag = NoFilterFlag- { unNoFilterFlag :: Bool- } deriving (Read,Show)+newtype LSADim = LSADim+ { unLSADim :: Int+ } deriving (Eq,Ord,Read,Show)+newtype SVDDim = SVDDim+ { unSVDDim :: Int+ } deriving (Eq,Ord,Read,Show) newtype ShiftPositiveFlag = ShiftPositiveFlag { unShiftPositiveFlag :: Bool } deriving (Read,Show)@@ -66,12 +106,25 @@ newtype MatrixTranspose = MatrixTranspose { unMatrixTranspose :: Bool } deriving (Read,Show)+newtype ChrRegion = ChrRegion+ { unChrRegion :: (T.Text, IntervalMap.Interval Int)+ } deriving (Read,Show,Eq,Ord)+newtype ChrRegionBin = ChrRegionBin+ { unChrRegionBin :: ChrRegion+ } deriving (Read,Show,Eq,Ord)+newtype ChrRegionMat = ChrRegionMat+ { unChrRegionMat+ :: Map.Map T.Text (IntervalMap.IntervalMap Int (S.SpVector Double))+ } deriving (Show)+newtype CustomRegions = CustomRegions+ { unCustomRegions :: [ChrRegion]+ } deriving (Read,Show) newtype X = X { unX :: Double- } deriving (Eq,Ord,Read,Show,Num,Generic,A.ToJSON,A.FromJSON)+ } deriving (Eq,Ord,Read,Show,Generic) deriving newtype (Num) deriving anyclass (A.ToJSON,A.FromJSON) newtype Y = Y { unY :: Double- } deriving (Eq,Ord,Read,Show,Num,Generic,A.ToJSON,A.FromJSON)+ } deriving (Eq,Ord,Read,Show,Generic) deriving newtype (Num) deriving anyclass (A.ToJSON,A.FromJSON) newtype ProjectionMap = ProjectionMap { unProjectionMap :: Map.Map Cell (X, Y) }@@ -81,47 +134,240 @@ newtype RMatScaled s = RMatScaled { unRMatScaled :: R.SomeSEXP s } newtype Row = Row { unRow :: Int- } deriving (Eq,Ord,Read,Show,Generic,A.ToJSON,A.FromJSON)+ } deriving (Eq,Ord,Read,Show,Generic) deriving anyclass (A.ToJSON,A.FromJSON) newtype Rows = Rows { unRows :: [Double] } newtype Vals = Vals { unVals :: [Double] }++newtype LabelMat = LabelMat { unLabelMat :: (SingleCells, Maybe LabelMap) }+newtype AggSingleCells = AggSingleCells { unAggSingleCells :: SingleCells }+ deriving (Show)+newtype AggReferenceMat = AggReferenceMat { unAggReferenceMat :: AggSingleCells }+ deriving (Show)+instance Semigroup LabelMat where+ (<>) (LabelMat (!m1, !l1)) (LabelMat (!m2, !l2)) =+ LabelMat (mappend m1 m2, mappend l1 l2)+instance Monoid LabelMat where+ mempty = LabelMat (mempty, mempty)+ newtype MatObsRow = MatObsRow { unMatObsRow :: S.SpMatrix Double } deriving (Show)+instance Semigroup MatObsRow where+ (<>) (MatObsRow x) (MatObsRow y) = MatObsRow $ S.vertStackSM x y+instance Monoid MatObsRow where+ mempty = MatObsRow $ S.zeroSM 0 0 -- Advanced+-- | Parse a chromosome region feature into an interval map.+parseChrRegion :: T.Text -> Either String ChrRegion+parseChrRegion x = either (\a -> Left $ a <> " (did you follow the format chrN:LOWERBOUND-UPPERBOUND ?): " <> T.unpack x) Right+ . Atto.parseOnly parserChrRegion+ . T.filter (/= ',')+ $ x++parserChrRegion :: Atto.Parser ChrRegion+parserChrRegion = do+ Atto.string "chr"+ chr <- Atto.takeWhile1 (/= ':')+ Atto.char ':'+ start <- Atto.decimal+ Atto.char '-'+ end <- Atto.decimal++ return . ChrRegion $ (chr, IntervalMap.ClosedInterval start end)++-- | Check if a matrix is a ChrRegion matrix based on the first feature.+isChrRegionMat :: (MatrixLike m) => m -> Bool+isChrRegionMat =+ maybe False (isRight . parseChrRegion) . flip (V.!?) 0 . getColNames++instance TextShow ChrRegion where+ showb (ChrRegion (chr, i)) = "chr"+ <> TS.fromText chr+ <> ":"+ <> showb (IntervalMap.lowerBound i)+ <> "-"+ <> showb (IntervalMap.upperBound i)++instance Monoid ChrRegionMat where+ mempty = ChrRegionMat Map.empty++instance Semigroup ChrRegionMat where+ (<>) (ChrRegionMat x) (ChrRegionMat y) =+ ChrRegionMat+ . Map.unionWith+ ( \x y+ -> IntervalMap.flattenWith (S.^+^)+ . IntervalMap.unionWith (S.^+^) x+ $ y+ )+ x+ $ y+ data SingleCells = SingleCells { _matrix :: MatObsRow , _rowNames :: Vector Cell- , _colNames :: Vector Gene+ , _colNames :: Vector Feature } deriving (Show) L.makeLenses ''SingleCells -data CellInfo = CellInfo- { _barcode :: Cell- , _cellRow :: Row- } deriving (Eq,Ord,Read,Show,Generic,A.ToJSON,A.FromJSON)-L.makeLenses ''CellInfo+instance Semigroup SingleCells where+ (<>) x y+ | (null . L.view colNames $ x)+ || (null . L.view colNames $ y)+ || (L.view colNames x == L.view colNames y) = -- If either is null or they share features, then do a normal append+ SingleCells { _matrix = mappend (L.view matrix x) (L.view matrix y)+ , _rowNames = (V.++) (L.view rowNames x) (L.view rowNames y)+ , _colNames = _colNames x+ }+ | isChrRegionMat x =+ either (\a -> error $ "Invalid chromosome region notation for feature, must be in format `chrN:START-END`: " <> a) id+ $ mergeChrFeaturesSingleCells x y+ | otherwise = mergeFeaturesSingleCells x y -data NormType = TfIdfNorm | UQNorm | MedNorm | TotalMedNorm | BothNorm | LogCPMNorm | NoneNorm- deriving (Read, Show, Eq)+-- | Join two SingleCells with different features.+mergeFeaturesSingleCells :: SingleCells -> SingleCells -> SingleCells+mergeFeaturesSingleCells sc1 sc2 =+ SingleCells { _matrix = MatObsRow+ . force+ $ (newMat id sc1)+ S.^+^ (newMat (+ sc1NRows) sc2)+ , _rowNames = (V.++) (L.view rowNames sc1) (L.view rowNames sc2)+ , _colNames = V.fromList newCols+ }+ where+ newMat rowF sc = S.modifyKeysSM+ rowF+ (lookupCol (fmap unFeature $ L.view colNames sc))+ . S.SM (newNRows, HMap.size newColMap)+ . S.immSM+ $ getMatrix sc+ newNRows = (V.length . L.view rowNames $ sc1)+ + (V.length . L.view rowNames $ sc2)+ sc1NRows = S.nrows . getMatrix $ sc1+ lookupCol :: V.Vector T.Text -> Int -> Int+ lookupCol olds x = fromMaybe (error $ "mergeFeaturesSingleCells: No new index when joining cols: " <> show x)+ . (>>=) (olds V.!? x)+ $ flip HMap.lookup newColMap+ newColMap :: HMap.HashMap T.Text Int+ newColMap = HMap.fromList . flip zip [0..] . fmap unFeature $ newCols+ newCols = Set.toAscList . Set.union (colToSet sc1) $ colToSet sc2+ colToSet = Set.fromList . V.toList . L.view colNames -instance (Generic a) => Generic (Vector a)+-- | Join two SingleCells with different chromosome region features. Merge the+-- features such that overlapping regions become one region.+mergeChrFeaturesSingleCells :: SingleCells -> SingleCells -> Either String SingleCells+mergeChrFeaturesSingleCells sc1 sc2 = do+ let getScMap sc = ChrRegionMat+ . fmap (IntervalMap.flattenWith (S.^+^))+ -- . Map.unionsWith (\x y -> force x `par` force y `par` IntervalMap.unionWith (S.^+^) x y)+ . Map.unionsWith (IntervalMap.unionWith (S.^+^))+ . zipWith (\r (ChrRegion (c, i))+ -> Map.singleton c $ IntervalMap.singleton i r+ )+ (S.toRowsL . S.transpose $ getMatrix sc)+ <$> mapM parseChrRegion (V.toList $ getColNames sc)+ sc1NRows = (V.length . L.view rowNames $ sc1)+ newNRows = sc1NRows + (V.length . L.view rowNames $ sc2)+ updateRows rowF mat = MatObsRow+ . S.modifyKeysSM+ rowF+ id+ . S.SM (newNRows, S.ncols $ unMatObsRow mat)+ . S.immSM+ . unMatObsRow+ $ mat+ newSc1 = L.over matrix (updateRows id) sc1+ newSc2 = L.over matrix (updateRows (+ sc1NRows)) sc2 -instance Semigroup MatObsRow where- (<>) (MatObsRow x) (MatObsRow y) = MatObsRow $ S.vertStackSM x y+ sc1Map <- getScMap newSc1+ sc2Map <- getScMap newSc2 -instance Monoid MatObsRow where- mempty = MatObsRow $ S.zeroSM 0 0+ let newMap = sc1Map <> sc2Map+ newMat = S.fromColsL+ $ Map.elems (unChrRegionMat newMap) >>= IntervalMap.elems+ newCols =+ Map.toAscList (unChrRegionMat newMap)+ >>= ( \(!k, !v) -> zipWith+ (\ x y+ -> Feature+ . showt+ . ChrRegion+ $ (x, y)+ )+ (repeat k)+ $ IntervalMap.keys v+ ) -instance Semigroup SingleCells where- (<>) x y =- SingleCells { _matrix = mappend (_matrix x) (_matrix y)- , _rowNames = (V.++) (_rowNames x) (_rowNames y)- , _colNames = _colNames x- }+ return $+ SingleCells { _matrix = MatObsRow . force $ newMat+ , _rowNames = (V.++)+ (L.view rowNames sc1)+ (L.view rowNames sc2)+ , _colNames = V.fromList newCols+ } instance Monoid SingleCells where mempty = SingleCells { _matrix = mempty , _rowNames = V.empty , _colNames = V.empty }++instance MatrixLike SingleCells where+ getMatrix = unMatObsRow . _matrix+ getRowNames = fmap unCell . _rowNames+ getColNames = fmap unFeature . _colNames+++data CellInfo = CellInfo+ { _barcode :: Cell+ , _cellRow :: Row+ } deriving (Eq,Ord,Read,Show,Generic) deriving anyclass (A.ToJSON,A.FromJSON)+L.makeLenses ''CellInfo++-- | A range with a label, most ranges can vary with each other.+data Range = Range { rangeIdx :: Maybe RangeIdx+ , _rangeLabel :: Text+ , _rangeLowerBound :: Int+ , _rangeUpperBound :: Int+ }+ deriving (Eq, Ord, Read, Show)+L.makeFields ''Range+instance TextShow Range where+ showb (Range _ l lb ub) = mconcat [TS.fromText l, ":", showb lb, "-", showb ub]++-- | A bin with a label, bins are standard across all observations.+data Bin = Bin { binIdx :: Maybe BinIdx, _binLabel :: Text, _binLowerBound :: Int, _binUpperBound :: Int }+ deriving (Eq, Ord, Read, Show)+L.makeFields ''Bin+instance TextShow Bin where+ showb (Bin _ l lb ub)= mconcat [TS.fromText l, ":", showb lb, "-", showb ub]++-- | Map of ranges to their bins.+newtype RangeBinMap = RangeBinMap+ { unRangeBinMap :: IMap.IntMap BinIdx+ }++data NormType = TfIdfNorm+ | UQNorm+ | MedNorm+ | TotalMedNorm+ | TotalNorm+ | LogCPMNorm Double+ | QuantileNorm+ | NoneNorm+ deriving (Read, Show, Eq)++-- | Map of bins to their idxs in unsorted and un-uniqued list.+newtype BinIdxMap = BinIdxMap+ { unBinIdxMap :: Map.Map Bin BinIdx+ }++deriving stock instance Generic Feature+instance NFData Feature where rnf x = seq x ()++instance (NFData a) => NFData (S.SpMatrix a) where+ rnf all@(S.SM (!x, !y) !m) = seq all ()++instance (NFData a) => NFData (S.SpVector a) where+ rnf all@(S.SV !x !m) = seq all ()
src/TooManyCells/Matrix/Utility.hs view
@@ -23,35 +23,55 @@ , isCsvFile , getMatrixOutputType , matrixValidity+ , hashNub+ , decompressStreamAll+ , transposeSc+ , extractCellSc+ , extractCellsSc+ , removeCellsSc+ , extractCellV+ , aggSc ) where -- Remote import BirchBeer.Types+import Codec.Compression.GZip (compress) import Control.Monad.Managed (runManaged) import Control.Monad.State (MonadState (..), State (..), evalState, execState, modify)+import Control.Monad.IO.Class (MonadIO) import Data.Bool (bool) import Data.Char (toUpper) import Data.Function (on)-import Data.List (maximumBy)+import Data.Hashable (Hashable)+import Data.List (maximumBy, isInfixOf, foldl') import Data.Maybe (fromMaybe)-import Data.Matrix.MatrixMarket (Matrix(RMatrix, IntMatrix), Structure (..), writeMatrix)+import Data.Matrix.MatrixMarket (Matrix(RMatrix, IntMatrix), Structure (..), writeMatrix') import Data.Scientific (toRealFloat, Scientific)+import Data.Streaming.Zlib (WindowBits) import Language.R as R import Language.R.QQ (r) import System.FilePath ((</>)) import TextShow (showt)+import qualified Control.Foldl as Fold import qualified Control.Lens as L+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 Data.Clustering.Hierarchical as HC import qualified Data.Graph.Inductive as G+import qualified Data.HashSet as HSet+import qualified Data.IntMap.Strict as IMap+import qualified Data.IntervalMap.Strict as IntervalMap import qualified Data.Map.Strict as Map import qualified Data.Sequence as Seq+import qualified Data.Set as Set import qualified Data.Sparse.Common as S import qualified Data.Text as T import qualified Data.Text.Encoding as T import qualified Data.Text.IO as T import qualified Data.Text.Lazy as TL+import qualified Data.Text.Lazy.Encoding as TL import qualified Data.Text.Lazy.IO as TL import qualified Data.Text.Lazy.Read as TL import qualified Data.Vector as V@@ -60,6 +80,7 @@ import qualified Streaming.Cassava as Stream import qualified Streaming.Prelude as Stream import qualified Streaming.With.Lifted as SW+import qualified Streaming.Zip as S import qualified System.Directory as FP -- Local@@ -84,7 +105,7 @@ [r| library(Matrix) |] let rowNamesR = fmap (T.unpack . unCell) . V.toList . _rowNames $ sc- colNamesR = fmap (T.unpack . unGene) . V.toList . _colNames $ sc+ colNamesR = fmap (T.unpack . unFeature) . V.toList . _colNames $ sc mat <- fmap unRMatObsRow . matToRMat . _matrix $ sc @@ -169,18 +190,25 @@ -- Where to place output files. FP.createDirectoryIfMissing True folder - writeMatrix (folder </> "matrix.mtx")+ let writeCompressedFile x = BL.writeFile x . compress++ (=<<) (writeCompressedFile (folder </> "matrix.mtx.gz"))+ . writeMatrix' . spMatToMat- . (bool S.transposeSM id mt) -- To have cells as columns+ . (bool S.transposeSM id mt) -- whether to transpose . getMatrix $ mat- T.writeFile (folder </> "genes.tsv")+ writeCompressedFile (folder </> "features.tsv.gz")+ . TL.encodeUtf8+ . TL.fromStrict . T.unlines- . fmap (\x -> T.intercalate "\t" [x, x])+ . fmap (\x -> T.intercalate "\t" [x, x, "NA"]) . V.toList . getColNames $ mat- T.writeFile (folder </> "barcodes.tsv")+ writeCompressedFile (folder </> "barcodes.tsv.gz")+ . TL.encodeUtf8+ . TL.fromStrict . T.unlines . V.toList . getRowNames@@ -254,3 +282,101 @@ (rows, cols) = S.dimSM . getMatrix $ mat numCells = V.length . getRowNames $ mat numFeatures = V.length . getColNames $ mat++-- | Same as Control.Foldl.nub but using HashSet.+hashNub :: (Hashable a, Eq a) => Fold.Fold a [a]+hashNub = Fold.Fold step (HSet.empty, id) fin+ where+ step (!s, !r) a =+ if HSet.member a s+ then (s, r)+ else (HSet.insert a s, r . (a :))+ fin (_, !r) = r []+{-# INLINABLE hashNub #-}++-- | Keep decompressing a compressed bytestream until exhaused.+decompressStreamAll :: (MonadIO m) => WindowBits -> BS.ByteString m r -> BS.ByteString m r+decompressStreamAll w bs = S.decompress' w bs >>= go+ where+ go (Left bs) = S.decompress' w bs >>= go+ go (Right r) = return r+{-# INLINABLE decompressStreamAll #-}++-- | Transpose a SingleCells type.+transposeSc :: SingleCells -> SingleCells+transposeSc (SingleCells (MatObsRow mat) rows cols) =+ SingleCells+ (MatObsRow . S.transposeSM $ mat)+ (fmap (Cell . unFeature) cols)+ (fmap (Feature . unCell) rows)++-- | Get a single cell from a SingleCells type.+extractCellSc :: Cell -> SingleCells -> SingleCells+extractCellSc cell sc = L.set rowNames (V.singleton cell)+ . L.set matrix newMatrix+ $ sc+ where+ newMatrix = MatObsRow+ . S.fromRowsL+ . (:[])+ . flip S.extractRow (getCellIdx sc cell)+ . unMatObsRow+ . L.view matrix+ $ sc++-- | Get a list of single cells for each single cell from a SingleCells type.+extractCellsSc :: SingleCells -> [SingleCells]+extractCellsSc sc = zipWith+ (\ v c -> sc { _rowNames = V.singleton c+ , _matrix = MatObsRow $ S.fromRowsL [v]+ }+ )+ (S.toRowsL . unMatObsRow . L.view matrix $ sc)+ . V.toList+ . L.view rowNames+ $ sc++-- | Get a single cell vector from a SingleCells type.+extractCellV :: Cell -> SingleCells -> S.SpVector Double+extractCellV cell sc =+ flip S.extractRow (getCellIdx sc cell) . unMatObsRow . L.view matrix $ sc++-- | Remove single cells from a SingleCells type.+removeCellsSc :: [Cell] -> SingleCells -> SingleCells+removeCellsSc cells sc = L.set rowNames newRowNames+ . L.set matrix newMatrix+ $ 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+ newMatrix = MatObsRow+ . S.fromRowsL+ . fmap (\x -> flip S.extractRow x . unMatObsRow . L.view matrix $ sc)+ $ whitelistIdxs++-- | Get the index of a cell in a SingleCells type.+getCellIdx :: SingleCells -> Cell -> Int+getCellIdx sc x =+ fromMaybe (error $ "extractCellIdx: cell not found in matrix: " <> show x)+ . V.findIndex (== x)+ . L.view rowNames+ $ sc++-- | Aggregate SingleCells using average.+aggSc :: SingleCells -> AggSingleCells+aggSc sc = AggSingleCells+ . L.over matrix (MatObsRow . avgMat . unMatObsRow)+ . L.over rowNames ( fmap (Cell . (<> " Aggregate") . unCell)+ . V.take 1+ )+ $ sc+ where+ avgMat = S.fromListDenseSM numRows+ . fmap ((/ fromIntegral numCols) . foldl' (+) 0)+ . S.toRowsL+ (numRows, numCols) = S.dim . unMatObsRow . L.view matrix $ sc
+ src/TooManyCells/Motifs/FindMotif.hs view
@@ -0,0 +1,155 @@+{- TooManyCells.Motifs.FindMotif+Gregory W. Schwartz++Collects functions pertaining to finding motifs per cluster.+-}++{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TupleSections #-}++module TooManyCells.Motifs.FindMotif+ ( getMotif+ , getNodes+ ) where++-- Remote+import Control.Monad (mfilter)+import Data.Maybe (catMaybes, fromMaybe, isJust)+import System.Directory (getTemporaryDirectory)+import TextShow (showt)+import qualified Control.Foldl as Fold+import qualified Data.ByteString.Lazy.Char8 as B+import qualified Data.Csv as CSV+import qualified Data.Map.Strict as Map+import qualified Data.Text as T+import qualified Data.Text.Read as T+import qualified Text.Printf as TP+import qualified Turtle as TU+import qualified Turtle.Line as TU++-- Local+import TooManyCells.File.Types+import TooManyCells.Motifs.Types++newtype ColMatch = ColMatch T.Text+newtype Match = Match T.Text++-- | Get the temporary directory.+getTmpDir :: IO TU.FilePath+getTmpDir = fmap (TU.fromText . T.pack) getTemporaryDirectory++-- | Read a CSV to a shell.+readCsv :: TU.FilePath -> TU.Shell (Map.Map T.Text T.Text)+readCsv file = do+ contents <- TU.liftIO . B.readFile . T.unpack . TU.format TU.fp $ file++ TU.select+ . either error snd+ . CSV.decodeByName+ $ contents++-- | Cut a shell csv.+csvCut :: T.Text -> TU.Shell (Map.Map T.Text T.Text) -> TU.Shell (Maybe TU.Line)+csvCut c = fmap (fmap TU.unsafeTextToLine . Map.lookup c)++-- | Sort a csv.+csvSort :: T.Text+ -> TU.Shell (Map.Map T.Text T.Text)+ -> TU.Shell (Map.Map T.Text T.Text)+csvSort c = (=<<) TU.select . TU.sortOn (Map.lookup c)++-- | Sort a csv by a numeric column.+csvNumSort :: T.Text+ -> TU.Shell (Map.Map T.Text T.Text)+ -> TU.Shell (Map.Map T.Text T.Text)+csvNumSort c =+ (=<<) TU.select . TU.sortOn (fmap readNum . Map.lookup c)+ where+ readNum "Infinity" = 1 / 0 :: Double+ readNum "-Infinity" = -1 / 0 :: Double+ readNum x = either (error . (<> ": " <> T.unpack x)) fst . T.double $ x++-- | Match a column for a shell csv.+csvMatch :: ColMatch+ -> Match+ -> TU.Shell (Map.Map T.Text T.Text)+ -> TU.Shell (Map.Map T.Text T.Text)+csvMatch (ColMatch c) (Match m) =+ mfilter ((== m) . Map.findWithDefault m c)++-- | Get the nodes from a differential file+getNodes :: DiffFile -> IO [Node]+getNodes (DiffFile df)+ = fmap catMaybes+ . TU.reduce Fold.nub+ . (fmap . fmap) (either error (Node . fst) . T.decimal . TU.lineToText)+ . csvCut "node"+ . readCsv+ $ df++-- | Make a temporary fasta file for input into the motif command.+mkTmpFasta :: TU.FilePath+ -> DiffFile+ -> GenomeFile+ -> TopN+ -> Maybe Node+ -> TU.Shell ()+mkTmpFasta tmpFasta diffFile gf (TopN topN) node = TU.sh $ do+ tmpDir <- TU.liftIO getTmpDir+ tmpBed <- TU.mktempfile tmpDir "motif_input.bed"++ (=<<) (TU.output tmpBed . TU.select)+ . TU.reduce (Fold.lastN topN)+ . TU.sed ("-" *> pure "\t")+ . TU.sed (":" *> pure "\t")+ . fmap (fromMaybe (error "feature column not found."))+ . csvCut "feature"+ . csvNumSort "log2FC"+ . csvMatch (ColMatch "node") (Match $ maybe "0" (showt . unNode) node)+ . readCsv+ . unDiffFile+ $ diffFile++ TU.output tmpFasta+ . TU.inproc "bedtools" [ "getfasta"+ , "-fi", unGenomeFile gf+ , "-bed", TU.format TU.fp tmpBed+ ]+ $ mempty++getMotif :: DiffFile+ -> Maybe BackgroundDiffFile+ -> OutputPath+ -> MotifCommand+ -> GenomeFile+ -> TopN+ -> Maybe Node+ -> IO ()+getMotif diffFile bgDiffFile outPath mc gf topN node = TU.sh $ do+ tmpDir <- TU.liftIO getTmpDir+ tmpFasta <- TU.mktempfile tmpDir "motif_input.fasta"+ mkTmpFasta tmpFasta diffFile gf topN node++ tmpBgFasta <- TU.mktempfile tmpDir "motif_bg_input.fasta"+ maybe+ (return ())+ (\x -> mkTmpFasta tmpBgFasta (DiffFile . unBackgroundDiffFile $ x) gf topN node)+ bgDiffFile++ let cmd = if isJust bgDiffFile+ then+ T.pack $+ TP.printf+ (unMotifCommand mc)+ (TU.format TU.fp tmpFasta)+ (TU.format TU.fp . unOutputPath $ outPath)+ (TU.format TU.fp tmpBgFasta)+ else+ T.pack $+ TP.printf+ (unMotifCommand mc)+ (TU.format TU.fp tmpFasta)+ (TU.format TU.fp . unOutputPath $ outPath)++ TU.stdout . TU.inshell cmd $ mempty
+ src/TooManyCells/Motifs/Types.hs view
@@ -0,0 +1,26 @@+{- TooManyCells.Motifs.Types+Gregory W. Schwartz++Collects the motif types used in the program.+-}++{-# LANGUAGE StrictData #-}++module TooManyCells.Motifs.Types where++-- Remote+import qualified Data.Text as T+import qualified Turtle as TU++-- Local++-- Basic++newtype DiffFile = DiffFile { unDiffFile :: TU.FilePath }+newtype BackgroundDiffFile =+ BackgroundDiffFile { unBackgroundDiffFile :: TU.FilePath }+newtype OutputPath = OutputPath { unOutputPath :: TU.FilePath }+newtype GenomeFile = GenomeFile { unGenomeFile :: T.Text }+newtype MotifCommand = MotifCommand { unMotifCommand :: String }+newtype TopN = TopN { unTopN :: Int }+newtype Node = Node { unNode :: Int } deriving (Eq, Ord)
src/TooManyCells/Paths/Distance.hs view
@@ -62,7 +62,7 @@ removeN = filter (\(_, x) -> x /= n1) flipEdges _ (a1, n, l, a2) = (a2, n, l, a1) --- | Get the start of the diameter a subtree. Can specify which node.+-- | Get the start of the diameter of a subtree. Can specify which node. treeDiameterStart :: G.Graph gr => PathDistance -> G.Node -> gr (G.Node, b) Double -> G.Node@@ -77,13 +77,25 @@ . getGraphLeavesCycles [] gr $ n +-- | Get the shallowest leaf.+shallowLeaf :: (G.Graph gr, Show a)+ => PathDistance -> gr (G.Node, a) Double -> G.Node+shallowLeaf pd gr = fst+ . fromMaybe (error "No nodes in tree.")+ . headMay+ . sortBy (compare `on` snd)+ . fmap (getDistance pd gr 0 . fst)+ . F.toList+ . getGraphLeavesCycles [] gr+ $ 0+ -- | Get the distance of each item from a starting point. Make graph undirected. linearItemDistance :: (TreeItem a, Show a)- => FlipFlag -> PathDistance -> ClusterGraph a -> [(a, Double)]-linearItemDistance direction pd (ClusterGraph gr) =+ => ShallowFlag -> FlipFlag -> PathDistance -> ClusterGraph a -> [(a, Double)]+linearItemDistance shallow direction pd (ClusterGraph gr) = concatMap (uncurry assignItems)- . linearNodeDistance direction pd+ . linearNodeDistance shallow direction pd . G.undir . G.emap (fromMaybe 0 . L.view edgeDistance) $ gr@@ -107,17 +119,22 @@ -- | Get the distance of each leaf from a starting point. linearNodeDistance :: (G.Graph gr, Show a)- => FlipFlag -> PathDistance -> gr (G.Node, a) Double -> [(G.Node, Double)]-linearNodeDistance (FlipFlag direction) pd gr =+ => ShallowFlag+ -> FlipFlag+ -> PathDistance+ -> gr (G.Node, a) Double+ -> [(G.Node, Double)]+linearNodeDistance (ShallowFlag shallow) (FlipFlag direction) pd gr = flipDirection direction- . fmap (getDistance pd gr start . fst)+ . fmap (getDistance pd gr (start shallow) . fst) . F.toList . getGraphLeavesCycles [] gr $ 0 where flipDirection True = id flipDirection False = flipDistance- start = treeDiameterStart pd 0 gr+ start True = shallowLeaf pd gr+ start False = treeDiameterStart pd 0 gr -- | Get the distance of a base node n1 to another node n2. getDistance
src/TooManyCells/Paths/Plot.hs view
@@ -18,34 +18,49 @@ -- Remote import BirchBeer.Types+import Data.Char (toUpper)+import Data.Colour.SRGB (sRGB24show) import Language.R as R import Language.R.QQ (r)+import qualified Control.Lens as L import qualified Data.Graph.Inductive as G+import qualified Data.Map.Strict as Map import qualified Data.Text as T -- Local import TooManyCells.Paths.Types -- | Plot clusters.-plotPathDistanceR :: String -> Bandwidth -> [(Label, Double)] -> R s ()-plotPathDistanceR outputPlot (Bandwidth b) distances = do+plotPathDistanceR :: String+ -> LabelColorMap+ -> Bandwidth+ -> [(Label, Double)]+ -> R s ()+plotPathDistanceR outputPlot (LabelColorMap cm) (Bandwidth b) distances = do let xs = fmap snd distances ls = fmap (T.unpack . unLabel . fst) distances+ (cls, ccs) = unzip+ . fmap (L.over L._1 (T.unpack . unLabel) . L.over L._2 (fmap toUpper . sRGB24show))+ . Map.toAscList+ $ cm [r| suppressMessages(library(ggplot2)) suppressMessages(library(cowplot)) suppressMessages(library(RColorBrewer))+ suppressMessages(library(plyr)) df = data.frame(x = xs_hs, l = ls_hs)+ spikeDf = ddply(df, "l", here(summarise), groupSpike = density(x, adjust = b_hs)$x[which.max(density(x, adjust = b_hs)$y)]) p = ggplot(df, aes(x = x, color = l, fill = l)) + geom_density(adjust = b_hs, alpha = 0.1) ++ geom_vline(data = spikeDf, aes(xintercept = groupSpike, color = l), linetype="dashed") + xlab("Path distance") + ylab("Cell density") +- scale_color_discrete(guide = guide_legend(title = "")) +- guides(fill = FALSE) ++ scale_color_manual(guide = guide_legend(title = ""), aesthetics = c("color", "fill"), values = setNames(c(ccs_hs, "NA"), c(cls_hs, "#000000"))) ++ theme_cowplot() + theme(aspect.ratio = 1) - suppressMessages(ggsave(p, file = outputPlot_hs))+ suppressMessages(ggsave(p, file = outputPlot_hs, height = 8, width = 12)) |] return ()
src/TooManyCells/Paths/Types.hs view
@@ -15,6 +15,7 @@ -- Basic newtype FlipFlag = FlipFlag { unFlipFlag :: Bool }+newtype ShallowFlag = ShallowFlag { unShallowFlag :: Bool } newtype Bandwidth = Bandwidth { unBandwidth :: Double } -- Advanced
+ src/TooManyCells/Peaks/ClusterPeaks.hs view
@@ -0,0 +1,448 @@+{- TooManyCells.MakeTree.ClusterPeaks+Gregory W. Schwartz++Collects functions pertaining to printing cluster fragments and obtaining peaks+per cluster.+-}++{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TupleSections #-}++module TooManyCells.Peaks.ClusterPeaks+ ( saveClusterFragments+ , peakCallFiles+ ) where++-- Remote+import BirchBeer.Types (Cluster (..), Label, LabelMap (..), Id (..))+import Codec.Compression.GZip (compress, decompress)+import Control.Concurrent.Async.Pool (withTaskGroup, mapTasks)+import Control.Monad (unless, when, void)+import Control.Monad.Trans.Maybe (runMaybeT, MaybeT (..))+import Control.Monad.Trans.Resource (runResourceT, MonadResource)+import Data.Bool (bool)+import Data.List (intercalate)+import Data.List.Split (splitOn)+import Data.Maybe (mapMaybe, fromMaybe, isNothing, isJust, fromJust)+import Data.Streaming.Zlib (WindowBits (..))+import Data.Tuple (swap)+import GHC.Conc (getNumCapabilities)+import Safe (headMay, atMay)+import System.IO (hPutStrLn, stderr)+import System.IO.Temp (withSystemTempFile)+import System.Process (callCommand)+import Text.Read (readMaybe)+import TextShow (showt)+import Turtle+import Turtle.Line+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 Data.HashMap.Strict as HMap+import qualified Data.IntMap.Strict as IMap+import qualified Data.IntSet as ISet+import qualified Data.Map.Strict as Map+import qualified Data.Set as Set+import qualified Text.Printf as TP+import qualified Data.Text as T+import qualified Data.Text.Encoding as T+import qualified Streaming as S+import qualified Streaming.Prelude as S+import qualified System.Directory as FP+import qualified System.FilePath as FP+import qualified Filesystem.Path as FSP+import qualified System.IO as IO+import qualified Turtle.Bytes as TB++-- Local+import TooManyCells.File.Types+import TooManyCells.MakeTree.Types+import TooManyCells.Matrix.Types+import TooManyCells.Matrix.Utility+import TooManyCells.Peaks.Types++-- | External sort fragment file.+externalFragmentsSort :: IO.FilePath -> IO ()+externalFragmentsSort file =+ void . callCommand $ TP.printf "sort -k1,1 -k2,2n -o %s %s" file file++-- | External sort compressed fragment file.+externalCompressedFragmentsSort :: IO.FilePath -> IO ()+externalCompressedFragmentsSort file = do+ void . callCommand $ TP.printf "gzip -d %s" file+ void+ . callCommand+ $ TP.printf "cat %s | sort -k1,1 -k2,2n | gzip -c > %s" (FP.dropExtension file) file+ void . callCommand $ TP.printf "rm %s" (FP.dropExtension file)++-- | Get cluster fragment sizes from a split stream.+getClusterTotalMap :: (MonadResource m)+ => HMap.HashMap T.Text [Int]+ -> S.Stream (S.Of [T.Text]) m r+ -> m (IMap.IntMap Double)+getClusterTotalMap cellClusterMap =+ S.fold_ (\acc x -> IMap.insertWith (+) x 1 acc) mempty id+ . S.concat+ . S.mapMaybe getCluster+ where+ getCluster :: [T.Text] -> Maybe [Int]+ getCluster row = atMay row 3 >>= flip HMap.lookup cellClusterMap++-- | Get the genome coverage bedgraph.+toGenomecov :: GenomecovCommand+ -> GenomeFile+ -> IMap.IntMap Double+ -> IO.FilePath+ -> IO.FilePath+ -> IO (Maybe ())+toGenomecov (GenomecovCommand gcc) (GenomeFile gf) clusterTotalMap outputPath file = runMaybeT $ do+ output <- MaybeT+ . return+ . fmap ((FP.</>) outputPath . (<> "_cov.bdg"))+ . headMay+ . splitOn "."+ . FP.takeFileName+ $ file+ wigOut <- MaybeT+ . return+ . fmap ((FP.</>) outputPath . (<> "_cov.bw"))+ . headMay+ . splitOn "."+ . FP.takeFileName+ $ file++ cluster <- MaybeT+ . return+ $ readMaybe+ =<< (flip atMay 1 . splitOn "_" . FP.takeFileName $ file)++ scale <-+ MaybeT . return . fmap (1000000 /) $ IMap.lookup cluster clusterTotalMap++ _ <- MaybeT . fmap Just $ externalCompressedFragmentsSort file+ _ <- MaybeT . fmap Just . callCommand $ TP.printf gcc file gf scale output+ _ <- MaybeT . fmap Just $ externalFragmentsSort output+ _ <- MaybeT . fmap Just . callCommand $ TP.printf "bedGraphToBigWig %s %s %s" output gf wigOut+ return ()++-- | Write the line of the Fragment or BedGraph file, depending on+-- normalization.+writeLine :: (MonadResource m) => HMap.HashMap T.Text [Int]+ -> ClusterHandleMap+ -> [T.Text]+ -> m (Maybe ())+writeLine cellClusterMap (ClusterHandleMap chm) all@(chrom:start:end:barcode:duplicateCount:_) = S.liftIO . runMaybeT $ do+ clusters <- MaybeT . return $ HMap.lookup barcode cellClusterMap++ let outputRow = compress+ . flip B.append "\n"+ . B.intercalate "\t"+ $ fmap (B.fromStrict . T.encodeUtf8) all+ write :: Int -> MaybeT IO ()+ write c = do+ h <- MaybeT . return $ IMap.lookup c chm+ MaybeT . fmap Just $ B.hPutStr h outputRow++ -- Fragment file+ mapM_ write clusters+-- writeLine (Just ctm) cellClusterMap outputPath all@(chrom:start:end:barcode:_:_) = S.liftIO . runMaybeT $ do+-- cluster <- MaybeT . return $ HMap.lookup barcode cellClusterMap+-- clusterSize <- MaybeT . return $ IMap.lookup cluster ctm++-- let val = (1 / clusterSize) * 1000000 -- per million+-- file = outputPath FP.</> T.unpack ("cluster_" <> showt cluster <> "_fragments.bdg")++-- -- Write header if file does not exist.+-- outputPathExists <- MaybeT . fmap Just $ FP.doesFileExist file+-- MaybeT . fmap Just . unless outputPathExists $+-- B.appendFile file+-- . B.fromStrict+-- . T.encodeUtf8+-- $ "track type=bedGraph name=" <> "cluster_" <> showt cluster <> "\n"++-- -- BedGraph file+-- MaybeT+-- . fmap Just+-- . B.appendFile file+-- . B.intercalate "\t"+-- $ (fmap (B.fromStrict . T.encodeUtf8) [chrom, start, end, showt val])+-- <> ["\n"]+writeLine _ _ xs = error $ "writeLine: Unexpected number of columns, did you use the fragments.tsv.gz 10x format? Input: " <> show xs++-- | Get cluster handle map to have several handles open.+getClusterHandleMap :: IO.FilePath -> [Int] -> IO ClusterHandleMap+getClusterHandleMap outputPath =+ fmap (ClusterHandleMap . IMap.fromList)+ . mapM (\x -> IO.openBinaryFile (getFileName x) IO.AppendMode >>= return . (x,))+ where+ getFileName cluster =+ outputPath FP.</> T.unpack ("cluster_" <> showt cluster <> "_fragments.tsv.gz")++-- | Get the list of cluster output files.+getClusterFiles :: IO.FilePath -> [Int] -> [IO.FilePath]+getClusterFiles outputPath = fmap getFileName+ where+ getFileName cluster =+ outputPath FP.</> T.unpack ("cluster_" <> showt cluster <> "_fragments.tsv.gz")++-- | Save the fragments per cluster. Loads in the original fragments file to+-- stream out cluster specific fragments.+saveClusterFragments :: OutputDirectory+ -> BedGraphFlag+ -> GenomeFile+ -> GenomecovCommand+ -> Maybe LabelMap+ -> AllNodesFlag+ -> PeakNodes+ -> PeakNodesLabels+ -> ClusterResults+ -> Either MatrixFileType MatrixFileType+ -> IO ()+saveClusterFragments+ output+ (BedGraphFlag bgf)+ genome+ gcc+ lm+ (AllNodesFlag anf)+ (PeakNodes pkns)+ (PeakNodesLabels pknls)+ cr+ (Left (CompressedFragments (FragmentsFile file))) = do+ -- Where to place output files+ let outputPath = unOutputDirectory output FP.</> "cluster_fragments"++ restartDir outputPath++ -- Check for occurrence of label map.+ when (isNothing lm && (not $ IMap.null pknls)) $+ hPutStrLn IO.stderr "--labels-file argument required for --peak-node-labels. Ignoring filter ..."++ let checkLabel :: CellInfo -> Int -> Maybe Bool+ checkLabel cellInfo cluster = do+ validLabs <- IMap.lookup cluster pknls+ label <- join+ $ fmap+ ( Map.lookup (Id . unCell $ L.view barcode cellInfo)+ . unLabelMap+ )+ lm+ return $ Set.member label validLabs+ clusterAssocList =+ fmap+ (\ (!cellInfo, !c)+ -> ( unCell $ _barcode cellInfo+ , bool (filter (fromMaybe True . checkLabel cellInfo)) id (IMap.null pknls) -- Specific node label selection.+ . bool (filter (flip ISet.member pkns)) id (ISet.null pkns) -- Specific node selection.+ . fmap unCluster+ $ bool (take 1 c) c anf+ )+ )+ . _clusterList+ $ cr+ cellClusterMap :: HMap.HashMap T.Text [Int]+ cellClusterMap = HMap.fromList clusterAssocList+ processedStream = S.map (T.splitOn "\t" . T.decodeUtf8)+ . S.mapped BS.toStrict+ . BS.lines+ . decompressStreamAll (WindowBits 31) -- For gunzip+ . BS.readFile+ $ file+ clusters = Set.toList+ . Set.fromList+ . mconcat+ . HMap.elems+ $ cellClusterMap++ -- Open the cluster file handles.+ clusterHandleMap <- getClusterHandleMap outputPath clusters++ hPutStrLn IO.stderr "Streaming fragments to cluster fragment files ..."+ -- Stream in file for fragments.+ runResourceT+ . S.effects+ . S.mapMaybeM (writeLine cellClusterMap clusterHandleMap)+ . S.filter (maybe False (not . T.null) . headMay)+ $ processedStream++ hPutStrLn IO.stderr "Closing cluster fragment files ..."+ -- Close the cluster file handles.+ mapM_ IO.hClose . unClusterHandleMap $ clusterHandleMap++ -- Sort cluster files.+ hPutStrLn IO.stderr "Sorting cluster fragment files ..."+ mapM_ externalCompressedFragmentsSort . getClusterFiles outputPath $ clusters++ -- Stream in file for bedgraph.+ when bgf $ do+ hPutStrLn IO.stderr "Creating bedgraphs ..."+ let bedPath = unOutputDirectory output FP.</> "cluster_bedgraphs"++ restartDir bedPath++ clusterTotalMap <-+ runResourceT . getClusterTotalMap cellClusterMap $ processedStream++ files <- fmap (fmap (outputPath FP.</>)) . FP.listDirectory $ outputPath+ cores <- getNumCapabilities+ withTaskGroup cores $ \workers ->+ void . mapTasks workers . fmap (toGenomecov gcc genome clusterTotalMap bedPath) $ files -- mapTasks_ not working+ -- runResourceT+ -- . S.effects+ -- . S.mapMaybeM (writeLine (Just clusterTotalMap) cellClusterMap bedPath)+ -- $ processedStream+saveClusterFragments _ _ _ _ _ _ _ _ _ _ = error "saveClusterFragments: Does not supported this matrix type. See too-many-cells -h for each entry point for more information"++-- | Run the specified command on a file. Command should have \"%s\" where the+-- output file and then input file (order is important) should be.+-- For instance: `command --output %s --other --args --input %s`.+peakCallFile :: PeakCommand+ -> GenomeFile+ -> IO.FilePath+ -> IO.FilePath+ -> IO IO.FilePath+peakCallFile (PeakCommand c) (GenomeFile gf) outputDir file =+ withSystemTempFile "temp_fragments.tsv" $ \temp h -> do+ let output = (FP.</>) outputDir+ . intercalate ("_")+ . (<> ["peaks"])+ . take 1+ . splitOn "_fragments.tsv.gz"+ . FP.takeFileName+ $ file+ wigIn = output FP.</> (name <> "_treat_pileup.bdg")+ wigOut = output FP.</> (name <> "_treat_pileup.bw")+ name = fromMaybe (error "peakCallFile: Bad file name for fragments.")+ . headMay+ . splitOn "_fragments.tsv.gz"+ . FP.takeFileName+ $ file+ B.readFile file >>= B.hPut h . decompress >> IO.hClose h -- Decompress file in temporary directory+ _ <- callCommand $ TP.printf c temp name output+ _ <- externalFragmentsSort wigIn+ -- _ <- callCommand $ TP.printf "bedGraphToBigWig %s %s %s" wigIn gf wigOut+ return $ output FP.</> (name <> "_peaks.narrowPeak")++-- | Convert between filepaths.+toTurtleFile :: IO.FilePath -> Turtle.FilePath+toTurtleFile = Turtle.fromText . T.pack++-- | Convert narrowPeak to bdg format.+narrowPeakToBdg :: Turtle.FilePath -> IO Turtle.FilePath+narrowPeakToBdg file = do+ let outFile = FSP.replaceExtension file "bdg"++ sh $+ output outFile+ . fmap fromJust+ . mfilter isJust+ . fmap ( join+ . fmap (textToLine . T.intercalate "\t")+ . (\xs -> mapM (atMay xs) [0, 1, 2, 4])+ . T.splitOn "\t"+ . lineToText+ )+ . input+ $ file++ return outFile++-- | Return the union of peaks from all requested nodes.+unionPeaks :: OutputDirectory -> [IO.FilePath] -> IO IO.FilePath+unionPeaks (OutputDirectory path) files = do+ let outputPath = toTurtleFile $ path+ outputFile = outputPath Turtle.</> "union.bdg"+ union [x] _ = input x+ union xs stdin = inproc "bedtools" (["unionbedg", "-i"] <> fmap (format fp) xs) stdin++ bdgFiles <- mapM (narrowPeakToBdg . Turtle.fromText . T.pack) files++ sh $ do+ mktree outputPath++ output outputFile+ . inproc "bedtools" ["merge", "-i", "stdin"]+ . union bdgFiles+ $ mempty++ return . T.unpack . format fp $ outputFile++-- | Return the original fragments file with unioned peaks.+unionPeaksFragments :: OutputDirectory+ -> IO.FilePath+ -> [Either MatrixFileType MatrixFileType]+ -> IO ()+unionPeaksFragments (OutputDirectory path) unionFile files = sh $ do+ let outputFile = toTurtleFile $ path FP.</> "union_fragments.tsv.gz"+ getInputFile (Left (CompressedFragments (FragmentsFile file))) = file+ getInputFile _ = error "unionPeaksFragments: Does not supported this matrix type. See too-many-cells -h for each entry point for more information"++ TB.output outputFile+ . TB.compress 6 (TB.WindowBits 31)+ . TB.fromUTF8+ . fmap ( flip T.append "\n"+ . T.intercalate "\t"+ . fromMaybe (error "unionPeaksFragments: Not enough columns.")+ . (\xs -> mapM (atMay xs) [0, 1, 2, 6, 7])+ . T.splitOn "\t"+ . lineToText+ )+ . inproc "bedtools" ( ["intersect", "-wa", "-wb", "-a", T.pack unionFile, "-b"]+ <> fmap (T.pack . getInputFile) files+ )+ $ mempty++-- | Run a specified command on all cluster files asynchronously.+peakCallFiles :: PeakCommand+ -> GenomeFile+ -> OutputDirectory+ -> [Either MatrixFileType MatrixFileType]+ -> IO ()+peakCallFiles pc gf (OutputDirectory path) fragmentsFiles = do+ let output = path FP.</> "cluster_peaks"+ fragmentsPath = path FP.</> "cluster_fragments"++ restartDir output++ files <- fmap (fmap (fragmentsPath FP.</>)) . FP.listDirectory $ fragmentsPath++ cores <- getNumCapabilities+ peakFiles <- withTaskGroup cores $ \workers ->+ mapTasks workers . fmap (peakCallFile pc gf output) $ files -- mapTasks_ not working++ unionFile <- unionPeaks (OutputDirectory output) peakFiles+ unionPeaksFragments (OutputDirectory output) unionFile fragmentsFiles++-- | Delete a directory if it exists.+restartDir :: IO.FilePath -> IO ()+restartDir dir = do+ dirExists <- FP.doesDirectoryExist dir+ when dirExists $ FP.removeDirectoryRecursive dir+ FP.createDirectoryIfMissing True dir+ return ()++-- -- | Filter a BED file by labels.+-- labelFilterBed :: LabelMap+-- -> Set.Set Label+-- -> Either MatrixFileType MatrixFileType+-- -> FP.FilePath+-- -> IO FP.FilePath+-- labelFilterBed lm labels (Left (CompressedFragments (FragmentsFile file))) inFile = sh $ do+-- let validCells = Set.fromList+-- . fmap fst+-- . filter (flip Set.member labels . snd)+-- . Map.toAscList+-- . unLabelMap+-- $ lm++-- out <- mktempfile ".temp" "temp.bed"++-- output out+-- . inproc "bedtools" ["intersect", "-a", "stdin", "-b", inFile]+-- . mfilter (maybe False (flip Set.member validCells) . atMay 3 . T.splitOn "\t")+-- . toLines+-- . TB.decompress (WindowBits 31)+-- . TB.input+-- $ file
+ src/TooManyCells/Peaks/Types.hs view
@@ -0,0 +1,29 @@+{- TooManyCells.Peaks.Types+Gregory W. Schwartz++Collects the peak types used in the program.+-}++{-# LANGUAGE StrictData #-}++module TooManyCells.Peaks.Types where++-- Remote+import BirchBeer.Types (Label)+import System.IO (Handle (..))+import qualified Data.IntSet as ISet+import qualified Data.IntMap.Strict as IMap+import qualified Data.Set as Set++-- Local++-- Basic+newtype PeakCommand = PeakCommand { unPeakCommand :: String }+newtype GenomecovCommand = GenomecovCommand { unGenomecovCommand :: String }+newtype GenomeFile = GenomeFile { unGenomeFile :: String }+newtype BedGraphFlag = BedGraphFlag { unBedGraphFlag :: Bool }+newtype AllNodesFlag = AllNodesFlag { unAllNodesFlag :: Bool }+newtype ClusterHandleMap = ClusterHandleMap { unClusterHandleMap :: IMap.IntMap Handle }+newtype PeakNodes = PeakNodes { unPeakNodes :: ISet.IntSet }+newtype PeakNodesLabels =+ PeakNodesLabels { unPeakNodesLabels :: IMap.IntMap (Set.Set Label) }
+ src/TooManyCells/Program/Classify.hs view
@@ -0,0 +1,64 @@+{- TooManyCells.Program.Classify+Gregory W. Schwartz++Classify entry point for command line program.+-}++{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE OverloadedStrings #-}++module TooManyCells.Program.Classify where++-- Remote+import BirchBeer.Types (Label (..))+import Control.Monad (unless, when)+import Data.Maybe (fromMaybe)+import Options.Generic+import Safe (headMay)+import Text.Read (readMaybe)+import TextShow (showt)+import qualified Data.Text as T+import qualified Data.Text.IO as T+import qualified System.FilePath as FP++-- Local+import TooManyCells.Classify.Classify+import TooManyCells.Classify.Types+import TooManyCells.File.Types+import TooManyCells.Matrix.Types+import TooManyCells.Matrix.Utility+import TooManyCells.Program.LoadMatrix+import TooManyCells.Program.Options++-- | Classify path.+classifyMain :: Options -> IO ()+classifyMain opts = do+ let readOrErr err = fromMaybe (error err) . readMaybe+ refFiles' = unHelpful+ . referenceFile+ $ opts+ singleRefMatFlag' = unHelpful . singleReferenceMatrix $ opts+ getRef x = AggReferenceMat+ . aggSc+ . fst+ . fromMaybe (error $ "Could not load file in required field --matrix-path:" <> show x)+ <$> loadAllSSM (opts { matrixPath = Helpful [x] })+ 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]})+ (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+ refs <- getRefs refFiles'++ let classified = classifyCells sc refs+ outputRow (Cell !c, (_, 0)) = T.putStrLn $ c <> ",unknown,0.0"+ outputRow (Cell !c, (Label !l, !s)) = T.putStrLn $ c <> "," <> l <> "," <> showt s++ putStrLn "item,label,score"+ mapM_ outputRow classified
+ src/TooManyCells/Program/Differential.hs view
@@ -0,0 +1,201 @@+{- TooManyCells.Program.Differential+Gregory W. Schwartz++Differential entry point for command line program.+-}++{-# LANGUAGE QuasiQuotes #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE PackageImports #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TupleSections #-}++module TooManyCells.Program.Differential where++-- Remote+import BirchBeer.Load+import BirchBeer.Types+import BirchBeer.Utility+import Control.Monad (when, join)+import Control.Monad.Trans (liftIO)+import Data.Bool (bool)+import Data.Maybe (fromMaybe, isJust, isNothing)+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+import qualified Data.ByteString.Lazy.Char8 as B+import qualified Data.Set as Set+import qualified H.Prelude as H+import qualified System.Directory as FP+import qualified System.FilePath as FP++-- Local+import TooManyCells.Program.Options+import TooManyCells.Differential.Differential+import TooManyCells.Differential.Types+import TooManyCells.MakeTree.Types+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+ let readOrErr err = fromMaybe (error err) . readMaybe+ delimiter' =+ Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts+ labelsFile' =+ fmap LabelFile . unHelpful . labelsFile $ 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+ labels' = fmap ( DiffLabels+ . L.over L.both ( (\x -> bool (Just x) Nothing . Set.null $ x)+ . Set.fromList+ . fmap Label+ )+ . readOrErr "Cannot read --labels."+ )+ . unHelpful+ . labels+ $ opts+ (combined1, combined2) = combineNodesLabels nodes' labels'+ plotOutputR = fromMaybe "out.pdf" . unHelpful . plotOutput $ opts++ when (isNothing labelsFile' && isJust labels') $+ hPutStrLn stderr "Warning: labels requested with no label file, ignoring labels..."++ scRes <- loadAllSSM opts+ let processedSc = fmap fst scRes+ customLabelMap = join . fmap snd $ scRes++ labelMap <- if isJust labelsFile'+ then mapM (loadLabelData delimiter') $ labelsFile'+ else return customLabelMap++ let clInput = (FP.</> "cluster_list.json") . unPriorPath $ prior'+ treeInput = (FP.</> "cluster_tree.json") . unPriorPath $ prior'++ cr :: IO ClusterResults+ cr = loadClusterResultsFiles clInput treeInput++ gr <- treeToGraph+ . updateTreeRowBool updateTreeRows' processedSc+ . _clusterDend+ <$> cr++ H.withEmbeddedR defaultConfig $ H.runRegion $ do+ case features' of+ [] -> do+ case nodes' of+ (DiffNodes ([], [])) -> do+ res <- H.io+ . getAllDEGraphKruskalWallis+ seed'+ subsampleGroups'+ topN'+ labelMap+ (maybe (DiffLabels (Nothing, Nothing)) id labels')+ (extractSc processedSc)+ $ gr++ H.io . B.putStrLn . getAllDEStringKruskalWallis $ res+ (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'+ then do+ res <- H.io+ $ getDEGraphKruskalWallis+ seed'+ subsampleGroups'+ topN'+ labelMap+ (extractSc processedSc)+ combined1+ combined2+ gr++ H.io . B.putStrLn . getDEStringKruskalWallis $ res+ else do+ res <- getDEGraph+ seed'+ subsampleGroups'+ topN'+ labelMap+ (extractSc processedSc)+ combined1+ combined2+ gr++ H.io . B.putStrLn . getDEString $ res+ _ -> do+ let outputCsvR = FP.replaceExtension plotOutputR ".csv"++ diffPlot <- getSingleDiff+ seed'+ subsampleGroups'+ False+ violinFlag'+ noOutlierFlag'+ aggregate'+ separateNodes'+ separateLabels'+ labelMap+ (extractSc processedSc)+ combined1+ combined2+ features'+ gr+ [H.r| suppressMessages(write.csv(diffPlot_hs[[2]], file = outputCsvR_hs, row.names = FALSE, quote = FALSE)) |]+ [H.r| suppressMessages(ggsave(diffPlot_hs[[1]], file = plotOutputR_hs)) |]++ let normOutputR = FP.replaceBaseName+ plotOutputR+ (FP.takeBaseName plotOutputR <> "_scaled")+ normOutputCsvR = FP.replaceExtension normOutputR ".csv"++ diffNormPlot <- getSingleDiff+ seed'+ subsampleGroups'+ True+ violinFlag'+ noOutlierFlag'+ aggregate'+ separateNodes'+ separateLabels'+ labelMap+ (extractSc processedSc)+ combined1+ combined2+ features'+ gr+ [H.r| suppressMessages(write.csv(diffNormPlot_hs[[2]], file = normOutputCsvR_hs, row.names = FALSE, quote = FALSE)) |]+ [H.r| suppressMessages(ggsave(diffNormPlot_hs[[1]], file = normOutputR_hs)) |]++ return ()
+ src/TooManyCells/Program/Diversity.hs view
@@ -0,0 +1,129 @@+{- TooManyCells.Program.Diversity+Gregory W. Schwartz++Diversity entry point into program.+-}++{-# LANGUAGE QuasiQuotes #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE PackageImports #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TupleSections #-}++module TooManyCells.Program.Diversity where++-- Remote+import BirchBeer.Load+import BirchBeer.Types+import Control.Monad.Trans (liftIO)+import Data.Bool (bool)+import Data.Maybe (fromMaybe)+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+import qualified Data.Csv as CSV+import qualified Data.Text as T+import qualified Data.Text.Lazy.IO as T+import qualified Diagrams.Backend.Cairo as D+import qualified Diagrams.Prelude as D+import qualified H.Prelude as H+import qualified Plots as D+import qualified System.Directory as FP+import qualified System.FilePath as FP++-- Local+import TooManyCells.Program.Options+import TooManyCells.Diversity.Diversity+import TooManyCells.Diversity.Load+import TooManyCells.Diversity.Plot+import TooManyCells.Diversity.Types+import TooManyCells.File.Types+import TooManyCells.Program.Options++-- | Diversity path.+diversityMain :: Options -> IO ()+diversityMain opts = do+ let priors' =+ fmap PriorPath . unHelpful . priors $ opts+ delimiter' =+ Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts+ labelsFile' =+ fmap LabelFile . unHelpful . labelsFile $ 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++ -- Where to place output files.+ FP.createDirectoryIfMissing True . unOutputDirectory $ output'++ labelMap <- sequence . fmap (loadLabelData delimiter') $ labelsFile'++ pops <- fmap ( either+ (\err -> error $ err <> "\nEncountered error in population loading, aborting process.")+ id+ . sequence+ )+ . mapM (\x -> (fmap . fmap) (Label . T.pack . unPriorPath $ x,)+ . loadPopulation labelMap+ $ x+ )+ $ priors'++ popDiversities <-+ mapM+ (\ (l, pop) -> getPopulationDiversity+ l+ order'+ start'+ interval'+ endMay'+ pop+ )+ pops++ -- Output quantifications.+ B.writeFile (unOutputDirectory output' FP.</> "diversity.csv")+ . CSV.encodeDefaultOrderedByName+ $ popDiversities++ -- D.renderCairo (unOutputDirectory output' FP.</> "diversity.pdf") D.absolute+ -- . plotDiversity+ -- $ popDiversities++ -- D.renderCairo (unOutputDirectory output' FP.</> "chao1.pdf") D.absolute+ -- . plotChao1+ -- $ popDiversities++ -- D.renderCairo (unOutputDirectory output' FP.</> "rarefaction.pdf") D.absolute+ -- . plotRarefaction+ -- $ popDiversities++ -- Output plots.+ let colors = D.colorRamp (length pops) . D.brewerSet D.Set1 $ 9++ H.withEmbeddedR defaultConfig $ H.runRegion $ do+ let divFile = unOutputDirectory output' FP.</> "diversity.pdf"+ divPlot <- plotDiversityR colors popDiversities+ [H.r| suppressMessages(ggsave(divPlot_hs, file = divFile_hs)) |]++ -- let chao1File = unOutputDirectory output' FP.</> "chao_r.pdf"+ -- chao1Plot <- plotChao1R colors popDiversities+ -- [H.r| suppressMessages(ggsave(chao1Plot_hs, file = chao1File_hs)) |]++ let rarefactionFile = unOutputDirectory output' FP.</> "rarefaction.pdf"+ rarefactionPlot <- plotRarefactionR colors popDiversities+ [H.r| suppressMessages(ggsave(rarefactionPlot_hs, file = rarefactionFile_hs)) |]++ return ()++ return ()
+ src/TooManyCells/Program/Interactive.hs view
@@ -0,0 +1,88 @@+{- TooManyCells.Program.Interactive+Gregory W. Schwartz++Interactive entry point into program.+-}++{-# LANGUAGE QuasiQuotes #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE PackageImports #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TupleSections #-}++module TooManyCells.Program.Interactive where++-- Remote+import BirchBeer.Interactive+import BirchBeer.Load+import BirchBeer.Types+import Control.Monad (join)+import Control.Monad.Trans (liftIO)+import Data.Bool (bool)+import Data.Maybe (fromMaybe, isJust)+import Data.Tree (Tree (..))+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+import qualified Data.Vector as V+import qualified H.Prelude as H+import qualified System.Directory as FP+import qualified System.FilePath as FP++-- Local+import TooManyCells.File.Types+import TooManyCells.MakeTree.Types+import TooManyCells.MakeTree.Utility+import TooManyCells.Matrix.Types+import TooManyCells.Program.LoadMatrix+import TooManyCells.Program.Options+import TooManyCells.Program.Utility++-- | Interactive tree interface.+interactiveMain :: Options -> IO ()+interactiveMain opts = H.withEmbeddedR defaultConfig $ do+ let labelsFile' =+ fmap LabelFile . unHelpful . labelsFile $ opts+ prior' = maybe (error "\nRequires --prior") PriorPath+ . unHelpful+ . prior+ $ opts+ updateTreeRows' =+ UpdateTreeRowsFlag . not . unHelpful . noUpdateTreeRows $ opts+ delimiter' =+ Delimiter . fromMaybe ',' . unHelpful . delimiter $ opts++ scRes <- loadAllSSM opts+ let mat = fmap fst scRes+ customLabelMap = join . fmap snd $ scRes++ labelMap <- if isJust labelsFile'+ then mapM (loadLabelData delimiter') $ labelsFile'+ else return customLabelMap++ tree <- fmap ( updateTreeRowBool updateTreeRows' mat+ . either error id+ . A.eitherDecode+ )+ . B.readFile+ . (FP.</> "cluster_tree.json")+ . unPriorPath+ $ prior' :: IO (Tree (TreeNode (V.Vector CellInfo)))++ interactiveDiagram+ tree+ labelMap+ mat+ . fmap ( B2Matrix+ . L.over matrix (MatObsRow . unB2 . b1ToB2 . B1 . unMatObsRow)+ )+ $ mat++ return ()
+ src/TooManyCells/Program/LoadMatrix.hs view
@@ -0,0 +1,267 @@+{- TooManyCells.Program.LoadMatrix+Gregory W. Schwartz++Loading matrix for command line program.+-}++{-# LANGUAGE QuasiQuotes #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE PackageImports #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TupleSections #-}++module TooManyCells.Program.LoadMatrix where++-- Remote+import BirchBeer.Types+import Control.Concurrent.Async.Pool (withTaskGroup, mapReduce, wait)+import Control.Monad (when, join)+import Control.Monad.STM (atomically)+import Control.Monad.Trans (liftIO)+import Control.Monad.Trans.Maybe (MaybeT (..))+import Data.Bool (bool)+import Data.Either (isRight)+import Data.List (foldl')+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+import qualified Data.ByteString.Char8 as B+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 System.Directory as FP+import qualified System.FilePath as FP++-- Local+import TooManyCells.File.Types+import TooManyCells.Matrix.AtacSeq+import TooManyCells.Matrix.Types+import TooManyCells.Matrix.Preprocess+import TooManyCells.Matrix.Utility+import TooManyCells.Matrix.Load+import TooManyCells.Program.Options+import TooManyCells.Program.Utility+import qualified Data.Sparse.Common as S+import Control.Lens++-- | Load the single cell matrix, post-whitelist-filtered cells.+loadSSM :: Options -> FilePath -> IO SingleCells+loadSSM opts matrixPath' = do+ fileExist <- FP.doesFileExist matrixPath'+ directoryExist <- FP.doesDirectoryExist matrixPath'+ compressedFileExist <- FP.doesFileExist $ matrixPath' FP.</> "matrix.mtx.gz"++ matrixFile' <- getMatrixFileType matrixPath'++ let featuresFile' = FeatureFile+ $ matrixPath'+ FP.</> (bool "genes.tsv" "features.tsv.gz" compressedFileExist)+ cellsFile' = CellFile+ $ 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+ excludeFragments' = fmap (ExcludeFragments . T.pack)+ . unHelpful+ . excludeMatchFragments+ $ opts+ blacklistRegionsFile' = fmap BlacklistRegions+ . unHelpful+ . blacklistRegionsFile+ $ opts+ cellWhitelistFile' =+ fmap CellWhitelistFile . unHelpful . cellWhitelistFile $ opts++ cellWhitelist <- liftIO . sequence $ fmap getCellWhitelist cellWhitelistFile'++ let unFilteredSc =+ case matrixFile' of+ (Left (DecompressedMatrix file)) ->+ loadSparseMatrixDataStream delimiter' file+ (Left (CompressedFragments file)) -> do+ liftIO $ when (isNothing binWidth') $+ hPutStrLn stderr "\nWarning: No binwidth specified for fragments file\+ \ input. This will make the feature list extremely large\+ \ and may result in many outliers. Please see --binwidth.\+ \ Ignore this message if using peaks.\+ \ Continuing..."+ liftIO $ when (isNothing cellWhitelist) $+ hPutStrLn stderr "\nWarning: No cell whitelist specified for fragments file\+ \ input. This will use all barcodes in the file. Most times\+ \ this file contains barcodes that are not cells. Please see\+ \ --cell-whitelist-file. Continuing..."+ fmap (bool binarizeSc id . unNoBinarizeFlag $ noBinarizeFlag')+ . loadFragments cellWhitelist blacklistRegionsFile' excludeFragments' binWidth'+ $ file+ (Left file@(BedGraph _)) -> do+ liftIO $ when (isNothing binWidth') $+ hPutStrLn stderr "\nWarning: No binwidth specified for BedGraph file\+ \ input. This will make the feature list extremely large\+ \ and may result in many outliers. Please see --binwidth.\+ \ Ignore this message if using peaks.\+ \ Continuing..."+ fmap (bool binarizeSc id . unNoBinarizeFlag $ noBinarizeFlag')+ . loadBdgBW blacklistRegionsFile' excludeFragments' binWidth'+ $ file+ (Left file@(BigWig _)) -> do+ liftIO $ when (isNothing binWidth') $+ hPutStrLn stderr "\nWarning: No binwidth specified for bigWig file\+ \ input. This will make the feature list extremely large\+ \ and may result in many outliers. Please see --binwidth.\+ \ Ignore this message if using peaks.\+ \ Continuing..."+ fmap (bool binarizeSc id . unNoBinarizeFlag $ noBinarizeFlag')+ . loadBdgBW blacklistRegionsFile' excludeFragments' binWidth'+ $ file+ (Right (DecompressedMatrix file)) ->+ loadCellrangerData featureColumn' featuresFile' cellsFile' file+ (Right (CompressedMatrix file)) ->+ loadCellrangerDataFeatures featureColumn' featuresFile' cellsFile' file+ _ -> error "Does not supported this matrix type. See too-many-cells -h for each entry point for more information"++ let whiteListFilter Nothing = id+ whiteListFilter (Just wl) = filterWhitelistSparseMat wl+ -- Convert non-read data (all matrix formats) to bins. This is done+ -- separately from the read data as it's faster to do this process in the+ -- fragment reading before the matrix creation.+ windowSc Nothing sc = sc+ windowSc (Just bw) sc =+ case matrixFile' of+ (Left file@(CompressedFragments _)) -> sc+ (Left file@(BigWig _)) -> sc+ otherwise -> bool sc (rangeToBinSc bw sc) . isChrRegionMat $ sc+ transposeFunc = bool id transposeSc $ unTransposeFlag transpose' -- Whether to transpose matrix++ fmap (windowSc binWidth' . whiteListFilter cellWhitelist . transposeFunc) unFilteredSc++-- | 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+ shiftPositiveFlag' =+ ShiftPositiveFlag . unHelpful . 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++ liftIO $ when ((isJust pca' || isJust lsa' || isJust svd') && (elem TfIdfNorm normalizations')) $+ hPutStrLn stderr "\nWarning: Dimensionality reduction (creating negative numbers) with tf-idf\+ \ normalization may lead to NaNs or 0s before spectral\+ \ clustering (leading to svdlibc to hang or dense SVD\+ \ to error out)! Continuing..."++ cores <- MaybeT . fmap Just $ getNumCapabilities+ (whitelistSc, whitelistLM) <- MaybeT+ $ if null matrixPaths'+ then return Nothing+ else+ withTaskGroup cores $ \workers ->+ fmap ( Just+ . (L.over L._1 (fastBinJoinCols binWidth' False)) -- Possibly undo feature change for fast bin joining+ . unLabelMat+ )+ . join+ . atomically+ . 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+ $ matrixPaths'++ let sc = maybe+ id+ (\x -> filterNumSparseMat x)+ filterThresholds'+ $ whitelistSc+ normMat TfIdfNorm = tfidfScaleSparseMat+ normMat UQNorm = uqScaleSparseMat+ normMat MedNorm = medScaleSparseMat+ normMat TotalMedNorm = scaleSparseMat+ normMat TotalNorm = totalScaleSparseMat+ normMat (LogCPMNorm b) = logCPMSparseMat b+ normMat QuantileNorm = quantileScaleSparseMat+ normMat NoneNorm = id+ applyNorms mat = foldl' (\acc norm -> L.over matrix (normMat norm) acc) mat normalizations'+ processSc = L.over matrix (MatObsRow . S.sparsifySM . unMatObsRow)+ . L.over matrix ( bool id shiftPositiveMat+ $ unShiftPositiveFlag shiftPositiveFlag'+ )+ . (\m -> maybe m (flip (pcaSparseSc dropDimensionFlag') m) pca')+ . (\m -> maybe m (flip (lsaSparseSc dropDimensionFlag') m) lsa')+ . (\m -> maybe m (flip (svdSparseSc dropDimensionFlag') m) svd')+ . bool (transformChrRegions customRegions') id (null $ unCustomRegions customRegions')+ . applyNorms+ . bool id binarizeSc (unBinarizeFlag binarizeFlag')+ processedSc = processSc sc+ -- Filter label map if necessary.+ labelMap = (\ valid -> fmap ( LabelMap+ . Map.filterWithKey (\k _ -> Set.member k valid)+ . unLabelMap+ )+ $ whitelistLM+ )+ . Set.fromList+ . fmap (Id . unCell)+ . V.toList+ . L.view rowNames+ $ processedSc+++ -- Check for empty matrix.+ when (V.null . getRowNames $ processedSc) . error $ emptyMatErr "cells"++ liftIO $ when ( length matrixPaths' > 1+ && ( fromMaybe False+ . fmap (isRight . parseChrRegion)+ . flip (V.!?) 0+ $ getColNames processedSc+ )+ )+ $+ hPutStrLn stderr "\nNote: Detected chromosome region features.\+ \ Matrices were combined by overlapping features.\+ \ To disable this feature, make sure the feature names\+ \ are NOT in the form `chrN:START-END`. For instance,\+ \ just add any character to the beginning of the feature.\+ \ Continuing..."++ liftIO . mapM_ print . matrixValidity $ processedSc++ return (processedSc, labelMap)
+ src/TooManyCells/Program/MakeTree.hs view
@@ -0,0 +1,501 @@+{- TooManyCells.Program.MakeTree+Gregory W. Schwartz++MakeTree entrypoint into the program.+-}++{-# LANGUAGE QuasiQuotes #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE PackageImports #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TupleSections #-}++module TooManyCells.Program.MakeTree where++-- Remote+import Control.Concurrent+import BirchBeer.ColorMap+import BirchBeer.Interactive+import BirchBeer.Load+import BirchBeer.MainDiagram+import BirchBeer.Plot+import BirchBeer.Types+import BirchBeer.Utility+import Control.Monad (when, unless, join)+import Control.Monad.Trans (liftIO)+import Control.Monad.Trans.Maybe (MaybeT (..))+import Data.Bool (bool)+import Data.Colour.SRGB (sRGB24read)+import Data.Maybe (fromMaybe, isJust)+import Data.Monoid ((<>))+import Language.R as R+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)+import qualified "find-clumpiness" Types as Clump+import qualified Control.Lens as L+import qualified Data.Aeson as A+import qualified Data.ByteString.Lazy.Char8 as B+import qualified Data.Colour.Palette.BrewerSet as D+import qualified Data.Colour.Palette.Harmony as D+import qualified Data.Csv as CSV+import qualified Data.GraphViz as G+import qualified Data.Set as Set+import qualified Data.Text as T+import qualified Data.Text.Lazy.IO as T+import qualified Diagrams.Backend.Cairo as D+import qualified Diagrams.Prelude as D+import qualified H.Prelude as H+import qualified Plots as D+import qualified System.Directory as FP+import qualified System.FilePath as FP+import qualified System.ProgressBar as Progress++-- Local+import TooManyCells.File.Types+import TooManyCells.MakeTree.Clumpiness+import TooManyCells.MakeTree.Cluster+import TooManyCells.MakeTree.Load+import TooManyCells.MakeTree.Plot+import TooManyCells.MakeTree.Print+import TooManyCells.MakeTree.Types+import TooManyCells.MakeTree.Utility+import TooManyCells.Matrix.Load+import TooManyCells.Matrix.Types+import TooManyCells.Matrix.Utility+import TooManyCells.Program.LoadMatrix+import TooManyCells.Program.Options+import TooManyCells.Program.Utility++makeTreeMain :: Options -> IO ()+makeTreeMain opts = H.withEmbeddedR defaultConfig $ do+ let readOrErr err = fromMaybe (error err) . readMaybe+ matrixPaths' = unHelpful . matrixPath $ opts+ labelsFile' =+ fmap LabelFile . unHelpful . labelsFile $ opts+ prior' =+ fmap PriorPath . unHelpful . prior $ 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+ 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+ elbowCutoff' =+ fmap ( ElbowCutoff+ . readOrErr "Cannot read --elbow-cutoff."+ )+ . unHelpful+ . elbowCutoff+ $ opts+ customCut' = CustomCut . Set.fromList . unHelpful . customCut $ opts+ rootCut' = fmap RootCut . unHelpful . rootCut $ opts+ dendrogramOutput' = DendrogramFile+ . fromMaybe "dendrogram.svg"+ . unHelpful+ . dendrogramOutput+ $ opts+ matrixOutput' = fmap (getMatrixOutputType . (unOutputDirectory output' FP.</>))+ . unHelpful+ . matrixOutput+ $ opts+ labelMapOutputFlag' =+ LabelMapOutputFlag . unHelpful . labelsOutput $ opts+ fragmentsOutputFlag' =+ FragmentsOutputFlag . unHelpful . 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+ 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+ drawColors' = fmap ( CustomColors+ . fmap sRGB24read+ . (\x -> readOrErr "Cannot read --draw-colors." x :: [String])+ )+ . unHelpful+ . drawColors+ $ opts+ drawDiscretize' = (=<<) (\x -> either error Just+ . either+ (\ err -> either+ (\y -> Left $ finalError err y)+ (Right . SegmentColorMap)+ (readEither x :: Either String Int)+ )+ (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+ 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++ pb <- Progress.newProgressBar Progress.defStyle 10 (Progress.Progress 0 11 ())+ -- Increment progress bar.+ Progress.incProgress pb 1++ -- Load matrix once.+ scRes <- loadAllSSM opts+ let processedSc = fmap fst scRes+ customLabelMap = join . fmap snd $ scRes++ -- Increment progress bar.+ Progress.incProgress pb 1++ -- 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++ --R.withEmbeddedR R.defaultConfig $ R.runRegion $ do+ -- For r clustering.+ -- mat <- scToRMat processedSc+ -- clusterRes <- hdbscan mat+ -- clusterList <- clustersToClusterList sc clusterRes++ -- For agglomerative clustering.+ --let clusterResults = fmap hClust processedSc++ -- Load previous results or calculate results if first run.+ originalClusterResults <- case prior' of+ Nothing -> do+ (fullCr, _) <-+ hSpecClust dense' eigenGroup' numEigen' minModularity' numRuns'+ . extractSc+ $ processedSc++ return fullCr :: IO ClusterResults+ (Just x) -> do+ let clInput = unPriorPath x FP.</> "cluster_list.json"+ treeInput = unPriorPath x FP.</> "cluster_tree.json"++ -- Strict loading in order to avoid locked file.+ fmap (L.over clusterDend (updateTreeRowBool updateTreeRows' processedSc))+ $ loadClusterResultsFiles clInput treeInput++ -- Increment progress bar.+ Progress.incProgress pb 1++ let birchMat = processedSc+ birchSimMat =+ case (not . null . unHelpful . matrixPath $ opts, drawCollection') of+ (True, CollectionGraph{} ) ->+ Just+ . B2Matrix+ . L.over matrix (MatObsRow . unB2 . b1ToB2 . B1 . unMatObsRow)+ . extractSc+ $ processedSc+ _ -> Nothing++ -- Increment progress bar.+ Progress.incProgress pb 1++ let config :: BirchBeer.Types.Config CellInfo SingleCells+ config = BirchBeer.Types.Config+ { _birchLabelMap = labelMap+ , _birchMinSize = minSize'+ , _birchMaxStep = maxStep'+ , _birchMaxProportion = maxProportion'+ , _birchMinDistance = minDistance'+ , _birchMinDistanceSearch = minDistanceSearch'+ , _birchSmartCutoff = smartCutoff'+ , _birchElbowCutoff = elbowCutoff'+ , _birchCustomCut = customCut'+ , _birchRootCut = rootCut'+ , _birchOrder = Just order'+ , _birchDrawLeaf = drawLeaf'+ , _birchDrawCollection = drawCollection'+ , _birchDrawMark = drawMark'+ , _birchDrawNodeNumber = drawNodeNumber'+ , _birchDrawMaxNodeSize = drawMaxNodeSize'+ , _birchDrawMaxLeafNodeSize = drawMaxLeafNodeSize'+ , _birchDrawNoScaleNodes = drawNoScaleNodes'+ , _birchDrawLegendSep = drawLegendSep'+ , _birchDrawLegendAllLabels = drawLegendAllLabels'+ , _birchDrawPalette = drawPalette'+ , _birchDrawColors = drawColors'+ , _birchDrawDiscretize = drawDiscretize'+ , _birchDrawScaleSaturation = drawScaleSaturation'+ , _birchDrawFont = drawFont'+ , _birchDrawItemLineWeight = drawItemLineWeight'+ , _birchDrawBarBounds = drawBarBounds'+ , _birchTree = _clusterDend originalClusterResults+ , _birchMat = birchMat+ , _birchSimMat = birchSimMat+ }++ (plot, labelColorMap, itemColorMap, markColorMap, tree', gr') <- mainDiagram config++ -- Increment progress bar.+ Progress.incProgress pb 1++ -- Write results.+ clusterResults <- case prior' of+ Nothing -> do+ let clusterList' = treeToClusterList tree'+ cr' = ClusterResults clusterList' tree'++ return cr'++ (Just x) -> do+ let clusterList' = treeToClusterList tree'+ cr' = ClusterResults clusterList' tree'++ return cr'++ B.writeFile+ (unOutputDirectory output' FP.</> "cluster_list.json")+ . A.encode+ . _clusterList+ $ clusterResults+ B.writeFile+ (unOutputDirectory output' FP.</> "cluster_tree.json")+ . A.encode+ . _clusterDend+ $ clusterResults+ T.writeFile+ (unOutputDirectory output' FP.</> "graph.dot")+ . G.printDotGraph+ . G.graphToDot G.nonClusteredParams+ . unClusterGraph+ $ gr'+ B.writeFile+ (unOutputDirectory output' FP.</> "cluster_info.csv")+ . printClusterInfo+ $ gr'+ -- Write matrix+ mapM_ (\x -> writeMatrixLike (MatrixTranspose False) x . extractSc $ processedSc) matrixOutput'+ -- Write label map+ mapM_+ ( bool+ (const (return ()))+ ( B.writeFile (unOutputDirectory output' FP.</> "labels.csv")+ . printLabelMap+ )+ . unLabelMapOutputFlag+ $ labelMapOutputFlag'+ )+ labelMap+ -- Write fragments.tsv.gz+ when (unFragmentsOutputFlag fragmentsOutputFlag')+ $ mapM_+ (\x -> saveFragments output' x =<< mapM getMatrixFileType matrixPaths')+ processedSc+ -- Write node info+ B.writeFile+ (unOutputDirectory output' FP.</> "node_info.csv")+ . printNodeInfo labelMap+ $ gr'+ case labelMap of+ Nothing -> return ()+ (Just lm) ->+ -- Write cluster diversity+ case clusterDiversity order' lm clusterResults of+ (Left err) -> hPutStrLn stderr+ $ err+ <> "\nWarning: Problem in diversity, skipping cluster_diversity.csv output ..."+ (Right result) ->+ B.writeFile+ ( unOutputDirectory output'+ FP.</> "cluster_diversity.csv"+ )+ . printClusterDiversity+ $ result++ -- Increment progress bar.+ Progress.incProgress pb 1++ -- Header+ B.putStrLn $ "cell,cluster,path"++ -- Body+ B.putStrLn+ . CSV.encode+ . fmap (\ (!ci, !(c:cs))+ -> ( unCell . _barcode $ ci+ , showt $ unCluster c+ , T.intercalate "/" . fmap (showt . unCluster) $ c:cs+ )+ )+ . _clusterList+ $ clusterResults++ -- Increment progress bar.+ Progress.incProgress pb 1++ -- Plot only if needed and ignore non-tree analyses if dendrogram is+ -- supplied.+ H.runRegion $ do+ -- Calculations with plotting the label map (clumpiness).+ case labelMap of+ Nothing ->+ H.io $ hPutStrLn stderr "\nClumpiness requires labels for cells, skipping ..."+ (Just lm) -> do+ -- Get clumpiness.+ case treeToClumpList clumpinessMethod' lm . _clusterDend $ clusterResults of+ (Left err) -> H.io+ . hPutStrLn stderr+ $ err+ <> "\nWarning: Problem in clumpiness, skipping clumpiness output ..."+ (Right clumpList) -> do+ -- Save clumpiness to a file.+ H.io+ . B.writeFile (unOutputDirectory output' FP.</> "clumpiness.csv")+ . clumpToCsv+ $ clumpList++ -- Plot clumpiness.+ either (H.io . hPutStrLn stderr) id+ $ plotClumpinessHeatmapR+ (unOutputDirectory output' FP.</> "clumpiness.pdf")+ clumpList++ -- View cutting location for modularity.+ case minDistanceSearch' of+ Nothing -> return ()+ (Just _) -> plotRankedModularityR+ (unOutputDirectory output' FP.</> "modularity_rank.pdf")+ . L.view clusterDend+ $ clusterResults++ -- Increment progress bar.+ H.io $ Progress.incProgress pb 1++ -- Plot.+ H.io $ do+ -- cr <- clusterResults+ -- gr <- graph+ -- cm <- itemColorMap++ -- plot <- if drawDendrogram'+ -- then return . plotDendrogram legend drawLeaf' cm . _clusterDend $ cr+ -- else do+ -- plotGraph legend drawConfig cm markColorMap gr++ D.renderCairo+ ( unOutputDirectory output'+ FP.</> unDendrogramFile dendrogramOutput'+ )+ (D.mkHeight 1000)+ plot++ -- Increment progress bar.+ H.io $ Progress.incProgress pb 1++ -- Plot clustering of the projections.+ case projectionFile' of+ Nothing -> return ()+ (Just f) -> do+ -- Load projection map if provided.+ projectionMap <- H.io $ loadProjectionMap f++ plotClustersR+ (unOutputDirectory output' FP.</> "projection.pdf")+ projectionMap+ . _clusterList+ $ clusterResults++ -- Plot clustering with labels.+ case (labelMap, itemColorMap) of+ (Just lm, Just icm) ->+ plotLabelClustersR+ (unOutputDirectory output' FP.</> "label_projection.pdf")+ projectionMap+ lm+ icm+ (_clusterList clusterResults)+ _ -> return ()++ -- Increment progress bar.+ H.io $ Progress.incProgress pb 1++ return ()
+ src/TooManyCells/Program/MatrixOutput.hs view
@@ -0,0 +1,31 @@+{- TooManyCells.Program.MatrixOutput+Gregory W. Schwartz++MatrixOutput entrypoint into the program.+-}++{-# LANGUAGE OverloadedStrings #-}++module TooManyCells.Program.MatrixOutput where++-- Remote+import Options.Generic+import qualified System.Directory as FP++-- Local+import TooManyCells.Program.Options+import TooManyCells.Program.LoadMatrix+import TooManyCells.File.Types+import TooManyCells.Matrix.Utility+import TooManyCells.Matrix.Types++matrixOutputMain :: Options -> IO ()+matrixOutputMain opts = do+ let matOutput' = getMatrixOutputType . unHelpful . matOutput $ opts++ -- Load matrix once.+ scRes <- loadAllSSM opts+ let processedSc = fmap fst scRes++ -- Write matrix+ writeMatrixLike (MatrixTranspose False) matOutput' . extractSc $ processedSc
+ src/TooManyCells/Program/Motifs.hs view
@@ -0,0 +1,64 @@+{- TooManyCells.Program.Motifs+Gregory W. Schwartz++Motif entry point for command line program.+-}++{-# LANGUAGE OverloadedStrings #-}++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++-- Local+import TooManyCells.File.Types+import TooManyCells.Motifs.FindMotif+import TooManyCells.Motifs.Types+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+ diffFile' = DiffFile . TU.fromText . unHelpful . diffFile $ opts+ backgroundDiffFile' = fmap (BackgroundDiffFile . TU.fromText)+ . unHelpful+ . backgroundDiffFile+ $ opts+ outDir = OutputDirectory+ . fromMaybe "out"+ . unHelpful+ . output+ $ opts++ TU.mktree . TU.fromText . T.pack . unOutputDirectory $ outDir++ TU.sh $ do+ nodes <- TU.liftIO $ getNodes diffFile'+ node <- if null nodes then pure Nothing else TU.select . fmap Just $ nodes++ let outPath = OutputPath+ $ (TU.fromText . T.pack $ unOutputDirectory outDir)+ TU.</> (maybe mempty (TU.fromText . ("node_" <>) . showt . unNode) node)++ TU.liftIO+ $ getMotif+ diffFile'+ backgroundDiffFile'+ outPath+ motifCommand'+ genome'+ topN'+ node
+ src/TooManyCells/Program/Options.hs view
@@ -0,0 +1,305 @@+{- TooManyCells.Program.Options+Gregory W. Schwartz++Options for the command line program.+-}++{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TypeOperators #-}++module TooManyCells.Program.Options where++-- Remote+import Options.Generic+import qualified Data.Text as T++-- Local++-- | 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 argument. An example of homer: `/path/to/findMotifs.pl %s fasta %s`. Uses meme by default."+ , 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)++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++instance ParseRecord Options where+ parseRecord = parseRecordWithModifiers modifiers
+ src/TooManyCells/Program/Paths.hs view
@@ -0,0 +1,102 @@+{- TooManyCells.Program.Paths+Gregory W. Schwartz++Paths entry point into program.+-}++{-# LANGUAGE QuasiQuotes #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE PackageImports #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TupleSections #-}++module TooManyCells.Program.Paths where++-- Remote+import BirchBeer.ColorMap (getLabelColorMap)+import BirchBeer.Load+import BirchBeer.Types+import BirchBeer.Utility+import Data.Bool (bool)+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+import qualified H.Prelude as H+import qualified System.Directory as FP+import qualified System.FilePath as FP++-- Local+import TooManyCells.File.Types+import TooManyCells.MakeTree.Plot+import TooManyCells.Matrix.Types+import TooManyCells.Paths.Distance+import TooManyCells.Paths.Plot+import TooManyCells.Paths.Types+import TooManyCells.Program.Options+import TooManyCells.File.Types+import TooManyCells.Matrix.Types++-- | Paths path.+pathsMain :: Options -> IO ()+pathsMain opts = do+ let readOrErr err = fromMaybe (error err) . readMaybe+ labelsFile' =+ maybe (error "\nNeed a label file.") LabelFile+ . unHelpful+ . labelsFile+ $ opts+ prior' =+ maybe (error "\nNeed a prior path containing tree.") PriorPath+ . unHelpful+ . prior+ $ 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+ output' =+ OutputDirectory . fromMaybe "out" . unHelpful . output $ opts++ -- Where to place output files.+ FP.createDirectoryIfMissing True . unOutputDirectory $ output'++ -- Get the label map from a file.+ labelMap <- loadLabelData delimiter' labelsFile'++ -- Load previous results or calculate results if first run.+ tree <- fmap (either error id . A.eitherDecode)+ . B.readFile+ . (FP.</> "cluster_tree.json")+ . unPriorPath+ $ prior'++ let gr = treeToGraph tree+ pathDistances :: [(CellInfo, Double)]+ pathDistances = linearItemDistance shallow' direction' pathDistance' gr+ labeledPathDistances =+ labelItemDistance labelMap pathDistances++ H.withEmbeddedR defaultConfig $ H.runRegion $ do+ let labelColorMap = getLabelColorMap palette' labelMap+ plotPathDistanceR+ (unOutputDirectory output' FP.</> "path_distances.pdf")+ labelColorMap+ bandwidth'+ labeledPathDistances+ return ()++ return ()
+ src/TooManyCells/Program/Peaks.hs view
@@ -0,0 +1,99 @@+{- TooManyCells.Program.Peaks+Gregory W. Schwartz++Peaks entry point for command line program.+-}++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)+import Text.Read (readMaybe)+import qualified Control.Lens as L+import qualified Data.IntSet as ISet+import qualified Data.IntMap.Strict as IMap+import qualified Data.Map.Strict as Map+import qualified Data.Set as Set+import qualified System.Directory as FP+import qualified System.FilePath as FP++-- Local+import TooManyCells.File.Types+import TooManyCells.MakeTree.Load+import TooManyCells.MakeTree.Types+import TooManyCells.Matrix.Types+import TooManyCells.Peaks.ClusterPeaks+import TooManyCells.Peaks.Types+import TooManyCells.Program.Options+import TooManyCells.Program.LoadMatrix+import TooManyCells.Program.Utility++-- | Peaks path.+peaksMain :: Options -> IO ()+peaksMain opts = do+ let readOrErr err = fromMaybe (error err) . readMaybe+ fragmentPaths' = unHelpful . fragmentsPath $ opts+ prior' = PriorPath+ . fromMaybe (error "\nRequires a previous run to get the clusters.")+ . unHelpful+ . prior+ $ 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+ 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+ 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 ..."++ cr <- loadClusterResultsFiles clInput treeInput :: IO ClusterResults+ labelMap <- mapM (loadLabelData delimiter') labelsFile'++ unless (unHelpful $ skipFragments opts) $ do+ hPutStrLn stderr "Splitting fragment file by cluster ..."+ mapM_+ ( (=<<) ( saveClusterFragments+ output'+ bedgraphFlag'+ genome'+ genomecovCommand'+ labelMap+ allNodesFlag'+ peakNodes'+ peakNodesLabels'+ cr+ )+ . getMatrixFileType+ )+ fragmentPaths'++ hPutStrLn stderr "Calling peaks ..."+ peakCallFiles peakCommand' genome' output' =<< mapM getMatrixFileType fragmentPaths'
+ src/TooManyCells/Program/Utility.hs view
@@ -0,0 +1,70 @@+{- TooManyCells.Program.Utility+Gregory W. Schwartz++Utility functions for the command line program.+-}++module TooManyCells.Program.Utility where++-- Remote+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++-- Local+import TooManyCells.File.Types+import TooManyCells.Matrix.Types+import TooManyCells.Program.Options++-- | Read or return an error.+readOrErr :: (Read a) => String -> String -> a+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++-- | Get the file type of an input matrix. Returns either a left file (i.e. CSV)+-- or a right matrix market.+getMatrixFileType :: FilePath -> IO (Either MatrixFileType MatrixFileType)+getMatrixFileType path = do+ fileExist <- FP.doesFileExist path+ directoryExist <- FP.doesDirectoryExist path+ compressedFileExist <- FP.doesFileExist $ path FP.</> "matrix.mtx.gz"++ let fragmentsFile = (\ (x, y) -> isInfixOf "fragments" x && y == ".tsv.gz")+ . FP.splitExtensions+ . FP.takeFileName+ $ path+ matrixFile+ | fileExist && (FP.takeExtension path == ".bdg") = Left . BedGraph $ path+ | fileExist && (FP.takeExtension path == ".bw") = Left . BigWig $ path+ | fileExist && not fragmentsFile = Left . DecompressedMatrix . MatrixFile $ path+ | fileExist && fragmentsFile = Left . CompressedFragments . FragmentsFile $ path+ | directoryExist && not compressedFileExist = Right . DecompressedMatrix . MatrixFile $ path FP.</> "matrix.mtx"+ | directoryExist && compressedFileExist = Right . CompressedMatrix . MatrixFile $ path FP.</> "matrix.mtx.gz"+ | directoryExist = error "Cannot determine matrix pointed to, are there too many matrices here?"+ | otherwise = error "\nMatrix path does not exist."++ return matrixFile
too-many-cells.cabal view
@@ -1,134 +1,140 @@-cabal-version: >=1.10-name: too-many-cells-version: 0.2.2.2-license: GPL-3-license-file: LICENSE-copyright: 2019 Gregory W. Schwartz-maintainer: gsch@pennmedicine.upenn.edu-author: Gregory W. Schwartz-homepage: http://github.com/GregorySchwartz/too-many-cells#readme-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.-category: Bioinformatics-build-type: Simple--source-repository head- type: git- location: https://github.com/GregorySchwartz/too-many-cells+name: too-many-cells+version: 2.1.0.1+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+category: Bioinformatics+build-type: Simple+-- extra-source-files:+cabal-version: >=1.10 library- exposed-modules:- TooManyCells.Diversity.Diversity- TooManyCells.Differential.Differential- TooManyCells.Differential.Types- TooManyCells.Diversity.Load- TooManyCells.Diversity.Plot- TooManyCells.Diversity.Types- TooManyCells.File.Types- TooManyCells.MakeTree.Adjacency- TooManyCells.MakeTree.Clumpiness- TooManyCells.MakeTree.Cluster- TooManyCells.MakeTree.Load- TooManyCells.MakeTree.Plot- TooManyCells.MakeTree.Print- TooManyCells.MakeTree.Types- TooManyCells.Matrix.Load- TooManyCells.Matrix.Preprocess- TooManyCells.Matrix.Types- TooManyCells.Matrix.Utility- TooManyCells.Paths.Distance- TooManyCells.Paths.Plot- TooManyCells.Paths.Types- hs-source-dirs: src- default-language: Haskell2010- ghc-options: -O2- build-depends:- base >=4.7 && <5,- aeson >=1.4.0.0,- birch-beer >=0.2.2.0,- bytestring >=0.10.8.2,- cassava >=0.5.1.0,- colour >=2.3.4,- containers >=0.5.11.0,- deepseq >=1.4.3.0,- diagrams >=1.4,- diagrams-lib >=1.4.2.3,- diagrams-cairo >=1.4.1,- diagrams-graphviz >=1.4.1,- differential >=0.1.2.0,- directory >=1.3.1.5,- diversity >=0.8.1.0,- find-clumpiness >=0.2.3.1,- filepath >=1.4.2,- fgl >=5.6.0.0,- foldl >=1.4.2,- graphviz >=2999.20.0.2,- hierarchical-clustering >=0.4.6,- hierarchical-spectral-clustering >=0.5.0.1,- hmatrix >=0.19.0.0,- inline-r >=0.9.2,- lens >=4.16.1,- managed >=1.0.6,- matrix-market-attoparsec >=0.1.0.8,- modularity >=0.2.1.1,- mtl >=2.2.2,- palette >=0.3.0.1,- parallel >=3.2.1.1,- plots >=0.1.1.0,- safe >=0.3.17,- scientific >=0.3.6.2,- sparse-linear-algebra >=0.3.1,- split >=0.2.3.3,- statistics >=0.14.0.2,- streaming >=0.2.1.0,- streaming-bytestring >=0.1.6,- streaming-cassava >=0.1.0.1,- streaming-with >=0.2.2.1,- SVGFonts >=1.6.0.3,- temporary >=1.3,- text >=1.2.3.0,- text-show >=3.7.4,- vector >=0.12.0.1,- vector-algorithms >=0.7.0.1,- zlib >=0.6.2+ hs-source-dirs: src+ exposed-modules: TooManyCells.Classify.Classify+ , TooManyCells.Classify.Types+ , TooManyCells.Differential.Differential+ , TooManyCells.Differential.Types+ , TooManyCells.Diversity.Diversity+ , TooManyCells.Diversity.Load+ , TooManyCells.Diversity.Plot+ , TooManyCells.Diversity.Types+ , TooManyCells.File.Types+ , TooManyCells.MakeTree.Adjacency+ , TooManyCells.MakeTree.Clumpiness+ , TooManyCells.MakeTree.Cluster+ , TooManyCells.MakeTree.Load+ , TooManyCells.MakeTree.Plot+ , TooManyCells.MakeTree.Print+ , TooManyCells.MakeTree.Types+ , TooManyCells.MakeTree.Utility+ , TooManyCells.Matrix.AtacSeq+ , TooManyCells.Matrix.Load+ , TooManyCells.Matrix.Preprocess+ , TooManyCells.Matrix.Types+ , TooManyCells.Matrix.Utility+ , TooManyCells.Motifs.FindMotif+ , TooManyCells.Motifs.Types+ , TooManyCells.Paths.Distance+ , TooManyCells.Paths.Plot+ , TooManyCells.Paths.Types+ , TooManyCells.Peaks.ClusterPeaks+ , TooManyCells.Peaks.Types+ , TooManyCells.Program.Classify+ , TooManyCells.Program.Differential+ , TooManyCells.Program.Diversity+ , TooManyCells.Program.Interactive+ , TooManyCells.Program.LoadMatrix+ , TooManyCells.Program.MakeTree+ , TooManyCells.Program.MatrixOutput+ , TooManyCells.Program.Motifs+ , TooManyCells.Program.Options+ , TooManyCells.Program.Paths+ , TooManyCells.Program.Peaks+ , TooManyCells.Program.Utility+ build-depends: base >= 4.7 && < 5+ , IntervalMap+ , SVGFonts+ , aeson+ , async+ , async-pool+ , attoparsec+ , birch-beer+ , bytestring+ , cassava+ , colour+ , containers+ , deepseq+ , diagrams+ , diagrams-cairo+ , diagrams-graphviz+ , diagrams-lib+ , differential+ , directory+ , diversity+ , fgl+ , filepath+ , find-clumpiness+ , foldl+ , graphviz+ , hashable+ , hierarchical-clustering+ , hierarchical-spectral-clustering+ , hmatrix+ , hmatrix-svdlibc+ , inline-r+ , lens+ , managed+ , matrix-market-attoparsec+ , modularity+ , mtl+ , mwc-random+ , optparse-generic+ , palette+ , parallel+ , plots+ , process+ , resourcet+ , safe+ , scientific+ , sparse-linear-algebra+ , spectral-clustering >= 0.3.0.2+ , split+ , statistics+ , stm+ , streaming+ , streaming-bytestring+ , streaming-cassava+ , streaming-commons+ , streaming-utils+ , streaming-with+ , system-filepath+ , temporary+ , terminal-progress-bar+ , text+ , text-show+ , transformers+ , turtle >= 1.5.18+ , unordered-containers+ , vector+ , vector-algorithms+ , zlib+ ghc-options: -O2+ default-language: Haskell2010 executable too-many-cells- main-is: Main.hs- hs-source-dirs: app- default-language: Haskell2010- ghc-options: -threaded -rtsopts -O2- build-depends:- base >=4.11.1.0,- too-many-cells -any,- aeson >=1.4.0.0,- birch-beer >=0.2.2.0,- bytestring >=0.10.8.2,- cassava >=0.5.1.0,- colour >=2.3.4,- containers >=0.5.11.0,- diagrams-lib >=1.4.2.3,- diagrams-cairo >=1.4.1,- directory >=1.3.1.5,- fgl >=5.6.0.0,- filepath >=1.4.2,- find-clumpiness >=0.2.3.1,- graphviz >=2999.20.0.2,- hierarchical-spectral-clustering >=0.5.0.1,- inline-r >=0.9.2,- lens >=4.16.1,- matrix-market-attoparsec >=0.1.0.8,- modularity >=0.2.1.1,- mtl >=2.2.2,- optparse-generic >=1.3.0,- palette >=0.3.0.1,- plots >=0.1.1.0,- spectral-clustering >=0.3.0.2,- streaming >=0.2.1.0,- streaming-bytestring >=0.1.6,- streaming-utils >=0.1.4.7,- terminal-progress-bar >=0.2,- text >=1.2.3.0,- text-show >=3.7.4,- transformers >=0.5.5.0,- vector >=0.12.0.1+ hs-source-dirs: app+ main-is: Main.hs+ ghc-options: -threaded -rtsopts -O2+ build-depends: base+ , too-many-cells+ , optparse-generic+ default-language: Haskell2010++source-repository head+ type: git+ location: https://github.com/GregorySchwartz/too-many-cells