too-many-cells-2.1.0.1: src/TooManyCells/Program/MakeTree.hs
{- 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 ()