diff --git a/app/Main.hs b/app/Main.hs --- a/app/Main.hs +++ b/app/Main.hs @@ -92,9 +92,10 @@ , 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 | 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). 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." + , 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." @@ -108,7 +109,7 @@ , 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) 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." + , 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." @@ -133,7 +134,7 @@ , 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 | 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). 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." + , 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." @@ -148,7 +149,7 @@ , 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 | 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). 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." + , 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." @@ -211,6 +212,7 @@ 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 @@ -336,19 +338,19 @@ normMat MedNorm = medScaleSparseMat normMat TotalMedNorm = scaleSparseMat normMat BothNorm = scaleSparseMat + normMat LogCPMNorm = logCPMSparseMat normMat NoneNorm = id - processMat = ( bool id shiftPositiveMat + processMat = normMat normalization' . _matrix + processedSc = ( bool id shiftPositiveSc $ unShiftPositiveFlag shiftPositiveFlag' ) - . (\m -> maybe m (flip pcaDenseMat m) pca') - . normMat normalization' - . _matrix - processedSc = sc { _matrix = processMat sc } + . (\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_ print . matrixValidity $ processedSc + liftIO . mapM_ (hPutStrLn stderr) . matrixValidity $ processedSc return processedSc @@ -370,6 +372,7 @@ 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' = @@ -519,8 +522,8 @@ -- Load previous results or calculate results if first run. originalClusterResults <- case prior' of Nothing -> do - let (fullCr, _) = - hSpecClust dense' eigenGroup' normalization' numEigen' minModularity' + (fullCr, _) <- + hSpecClust dense' eigenGroup' normalization' numEigen' minModularity' numRuns' . extractSc $ processedSc diff --git a/src/TooManyCells/MakeTree/Cluster.hs b/src/TooManyCells/MakeTree/Cluster.hs --- a/src/TooManyCells/MakeTree/Cluster.hs +++ b/src/TooManyCells/MakeTree/Cluster.hs @@ -21,7 +21,7 @@ -- Remote import BirchBeer.Types -import BirchBeer.Utility (getGraphLeaves, getGraphLeavesWithParents, dendrogramToGraph, dendToTree, treeToGraph) +import BirchBeer.Utility (getGraphLeaves, getGraphLeavesWithParents, dendrogramToGraph, dendToTree, clusteringTreeToTree, treeToGraph) import Control.Monad (join) import Data.Function (on) import Data.List (sortBy, groupBy, zip4, genericLength) @@ -137,17 +137,82 @@ -> NormType -> Maybe NumEigen -> Maybe Q + -> Maybe NumRuns -> SingleCells - -> (ClusterResults, ClusterGraph CellInfo) -hSpecClust (DenseFlag isDense) eigenGroup norm numEigen minModMay sc = - ( ClusterResults { _clusterList = clustering - , _clusterDend = dendToTree dend - } - , gr - ) - where - clustering :: [(CellInfo, [Cluster])] - clustering = + -> IO (ClusterResults, ClusterGraph CellInfo) +hSpecClust (DenseFlag isDense) eigenGroup norm 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 = + hierarchicalSpectralCluster + eigenGroup + False + (fmap unNumEigen numEigen) + Nothing + minModMay + (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 = + HSD.hierarchicalSpectralCluster + eigenGroup + False + (fmap unNumEigen numEigen) + Nothing + minModMay + (fmap unNumRuns runsMay) + items + . Left + . sparseToHMat + + tree <- hSpecCommand norm isDense . unMatObsRow . _matrix $ sc + + let clustering :: [(CellInfo, [Cluster])] + clustering = concatMap (\ (!ns, (_, !xs)) -> zip (maybe [] F.toList xs) . repeat . fmap Cluster $ ns ) @@ -155,73 +220,14 @@ . flip getGraphLeavesWithParents 0 . unClusterGraph $ gr - gr = dendrogramToGraph dend - dend = clusteringTreeToDendrogram tree - tree = hSpecCommand norm isDense - . unMatObsRow - . _matrix - $ sc - 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 - items - . Left - hSpecCommand BothNorm False = - hierarchicalSpectralCluster - eigenGroup - True - (fmap unNumEigen numEigen) - Nothing - minModMay - items - . Left - hSpecCommand _ False = - hierarchicalSpectralCluster - eigenGroup - False - (fmap unNumEigen numEigen) - Nothing - minModMay - items - . Left - hSpecCommand TfIdfNorm True = - HSD.hierarchicalSpectralCluster - eigenGroup - True - (fmap unNumEigen numEigen) - Nothing - minModMay - items - . Left - . sparseToHMat - hSpecCommand BothNorm True = - HSD.hierarchicalSpectralCluster - eigenGroup - True - (fmap unNumEigen numEigen) - Nothing - minModMay - items - . Left - . sparseToHMat - hSpecCommand _ True = - HSD.hierarchicalSpectralCluster - eigenGroup - False - (fmap unNumEigen numEigen) - Nothing - minModMay - items - . Left - . sparseToHMat + dend = clusteringTreeToTree tree + gr = treeToGraph dend + + return ( ClusterResults { _clusterList = clustering + , _clusterDend = dend + } + , gr + ) dendrogramToClusterList :: HC.Dendrogram (V.Vector CellInfo) -> [(CellInfo, [Cluster])] diff --git a/src/TooManyCells/MakeTree/Print.hs b/src/TooManyCells/MakeTree/Print.hs --- a/src/TooManyCells/MakeTree/Print.hs +++ b/src/TooManyCells/MakeTree/Print.hs @@ -17,6 +17,7 @@ -- 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.Monoid ((<>)) @@ -62,7 +63,7 @@ getQs :: ([G.Node], a) -> [Double] getQs = mapMaybe getQ . fst getQ :: G.Node -> Maybe Double - getQ = fmap snd . headMay . G.lsuc gr + getQ = join . fmap (L.view (L._2. edgeDistance)) . headMay . G.lsuc gr -- | Get the information of each leaf cluster path. Modularity -- starts from the parent of the cluster to the root for modularity. @@ -102,12 +103,16 @@ (getSize x) (getProportion x) (getQ x) + (getSignificance x) (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 getQ :: G.Node -> Maybe Double - getQ = fmap snd . headMay . G.lsuc gr + getQ = join . fmap (L.view edgeDistance . snd) . headMay . G.lsuc gr + getSignificance :: G.Node -> Maybe Double + getSignificance = + join . fmap (L.view edgeSignificance . snd) . headMay . G.lsuc gr getProportion :: G.Node -> Maybe Double getProportion n = case fmap fst $ G.lsuc gr n of @@ -123,13 +128,14 @@ printNodeInfo :: (TreeItem a) => Maybe LabelMap -> ClusterGraph a -> B.ByteString printNodeInfo lm = - (<>) "node,size,proportion,modularity,composition,subtree\n" + (<>) "node,size,proportion,modularity,significance,composition,subtree\n" . CSV.encode - . fmap (\ (NodeInfo n s p m c ch) + . fmap (\ (NodeInfo n s p m sig c ch) -> ( n , s , maybe "" show p , maybe "" show m + , maybe "" show sig , maybe "" show c , intercalate "/" . fmap show $ ch ) diff --git a/src/TooManyCells/MakeTree/Types.hs b/src/TooManyCells/MakeTree/Types.hs --- a/src/TooManyCells/MakeTree/Types.hs +++ b/src/TooManyCells/MakeTree/Types.hs @@ -62,6 +62,9 @@ newtype NumEigen = NumEigen { unNumEigen :: Int } deriving (Read, Show) +newtype NumRuns = NumRuns + { unNumRuns :: Int + } deriving (Read, Show) newtype L = L Double newtype C = C Double newtype H = H Double @@ -90,6 +93,7 @@ , _nodeSize :: Int , _nodeSplitProportion :: Maybe Double , _nodeModularity :: Maybe Double + , _nodeSignificance :: Maybe Double , _nodeLabelComposition :: Maybe LabelCompositions , _nodeChildren :: [G.Node] } diff --git a/src/TooManyCells/Matrix/Preprocess.hs b/src/TooManyCells/Matrix/Preprocess.hs --- a/src/TooManyCells/Matrix/Preprocess.hs +++ b/src/TooManyCells/Matrix/Preprocess.hs @@ -6,11 +6,13 @@ {-# LANGUAGE BangPatterns #-} {-# LANGUAGE QuasiQuotes #-} +{-# LANGUAGE OverloadedStrings #-} module TooManyCells.Matrix.Preprocess ( scaleRMat , scaleDenseMat , scaleSparseMat + , logCPMSparseMat , uqScaleSparseMat , medScaleSparseMat , filterRMat @@ -21,8 +23,8 @@ , featureSelectionRandomForest , removeCorrelated , pcaRMat - , pcaDenseMat - , shiftPositiveMat + , pcaDenseSc + , shiftPositiveSc ) where -- Remote @@ -34,6 +36,7 @@ import Language.R as R import Language.R.QQ (r) import Statistics.Quantile (continuousBy, s) +import TextShow (showt) import qualified Control.Lens as L import qualified Data.Set as Set import qualified Data.Sparse.Common as S @@ -94,6 +97,16 @@ . 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 + -- | Scale a matrix based on the median. medScaleSparseMat :: MatObsRow -> MatObsRow medScaleSparseMat (MatObsRow mat) = MatObsRow @@ -116,6 +129,13 @@ where total = sum xs +-- | log(CPM + 1) normalization for cell. +logCPMSparseCell :: S.SpVector Double -> S.SpVector Double +logCPMSparseCell xs = fmap (logBase 2 . (+ 1) . cpm) xs + where + cpm x = x / tpm + tpm = sum xs / 1000000 + -- | Upper quartile scale cells. uqScaleSparseCell :: S.SpVector Double -> S.SpVector Double uqScaleSparseCell xs = fmap (/ uq) xs @@ -312,7 +332,7 @@ [r| mat = prcomp(t(mat_hs), tol = 0.95)$x |] --- | Conduct PCA on a matrix, retaining a percentage of variance. +-- | Conduct PCA on a matrix, taking the first principal components. pcaDenseMat :: PCADim -> MatObsRow -> MatObsRow pcaDenseMat (PCADim pcaDim) (MatObsRow matObs) = MatObsRow @@ -328,6 +348,12 @@ where mat = sparseToHMat matObs +-- | 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) + -- | Shift features to positive values. shiftPositiveMat :: MatObsRow -> MatObsRow shiftPositiveMat = MatObsRow @@ -341,3 +367,7 @@ . unMatObsRow where shift xs = S.sparsifySV . S.vr . fmap (\x -> x + minimum xs + 1) $ xs + +-- | Shift features to positive values for SingleCells. +shiftPositiveSc :: SingleCells -> SingleCells +shiftPositiveSc = L.over matrix shiftPositiveMat diff --git a/src/TooManyCells/Matrix/Types.hs b/src/TooManyCells/Matrix/Types.hs --- a/src/TooManyCells/Matrix/Types.hs +++ b/src/TooManyCells/Matrix/Types.hs @@ -102,7 +102,8 @@ } deriving (Eq,Ord,Read,Show,Generic,A.ToJSON,A.FromJSON) L.makeLenses ''CellInfo -data NormType = TfIdfNorm | UQNorm | MedNorm | TotalMedNorm | BothNorm | NoneNorm deriving (Read, Show, Eq) +data NormType = TfIdfNorm | UQNorm | MedNorm | TotalMedNorm | BothNorm | LogCPMNorm | NoneNorm + deriving (Read, Show, Eq) instance (Generic a) => Generic (Vector a) diff --git a/src/TooManyCells/Paths/Distance.hs b/src/TooManyCells/Paths/Distance.hs --- a/src/TooManyCells/Paths/Distance.hs +++ b/src/TooManyCells/Paths/Distance.hs @@ -85,6 +85,7 @@ concatMap (uncurry assignItems) . linearNodeDistance direction pd . G.undir + . G.emap (fromMaybe 0 . L.view edgeDistance) $ gr where assignItems n v = diff --git a/too-many-cells.cabal b/too-many-cells.cabal --- a/too-many-cells.cabal +++ b/too-many-cells.cabal @@ -1,6 +1,6 @@ cabal-version: >=1.10 name: too-many-cells -version: 0.1.12.4 +version: 0.2.0.0 license: GPL-3 license-file: LICENSE copyright: 2019 Gregory W. Schwartz @@ -46,7 +46,7 @@ build-depends: base >=4.7 && <5, aeson >=1.4.0.0, - birch-beer >=0.1.4.4, + birch-beer >=0.2.0.0, bytestring >=0.10.8.2, cassava >=0.5.1.0, colour >=2.3.4, @@ -65,7 +65,7 @@ foldl >=1.4.2, graphviz >=2999.20.0.2, hierarchical-clustering >=0.4.6, - hierarchical-spectral-clustering >=0.4.1.2, + hierarchical-spectral-clustering >=0.5.0.1, hmatrix >=0.19.0.0, inline-r >=0.9.2, lens >=4.16.1, @@ -102,7 +102,7 @@ base >=4.11.1.0, too-many-cells -any, aeson >=1.4.0.0, - birch-beer >=0.1.4.4, + birch-beer >=0.2.0.0, bytestring >=0.10.8.2, cassava >=0.5.1.0, colour >=2.3.4, @@ -114,7 +114,7 @@ filepath >=1.4.2, find-clumpiness >=0.2.3.1, graphviz >=2999.20.0.2, - hierarchical-spectral-clustering >=0.4.1.2, + hierarchical-spectral-clustering >=0.5.0.1, inline-r >=0.9.2, lens >=4.16.1, matrix-market-attoparsec >=0.1.0.8,