diff --git a/src/TooManyCells/Motifs/FindMotif.hs b/src/TooManyCells/Motifs/FindMotif.hs --- a/src/TooManyCells/Motifs/FindMotif.hs +++ b/src/TooManyCells/Motifs/FindMotif.hs @@ -10,11 +10,13 @@ module TooManyCells.Motifs.FindMotif ( getMotif + , getMotifGenome , getNodes ) where -- Remote import Control.Monad (mfilter) +import Data.List (foldl') import Data.Maybe (catMaybes, fromMaybe, isJust) import System.Directory (getTemporaryDirectory) import TextShow (showt) @@ -118,6 +120,32 @@ ] $ mempty +-- | Make a temporary bed file for input into the genome motif command. +mkTmpBed :: TU.FilePath + -> DiffFile + -> TopN + -> Maybe Node + -> TU.Shell () +mkTmpBed tmpBed diffFile (TopN topN) node = TU.sh $ do + tmpDir <- TU.liftIO getTmpDir + + let formatRow xs = -- Format recommended by findMotifsGenome.pl + TU.unsafeTextToLine + $ foldl' (\acc x -> T.replace x "\t" acc) (TU.lineToText xs) [":", "-"] + <> "\t" + <> TU.lineToText xs + <> "\t\t+" + + (=<<) (TU.output tmpBed . TU.select) + . TU.reduce (Fold.lastN topN) + . fmap (formatRow . fromMaybe (error "feature column not found.")) + . csvCut "feature" + . csvNumSort "log2FC" + . csvMatch (ColMatch "node") (Match $ maybe "0" (showt . unNode) node) + . readCsv + . unDiffFile + $ diffFile + getMotif :: DiffFile -> Maybe BackgroundDiffFile -> OutputPath @@ -150,6 +178,44 @@ TP.printf (unMotifCommand mc) (TU.format TU.fp tmpFasta) + (TU.format TU.fp . unOutputPath $ outPath) + + TU.stdout . TU.inshell cmd $ mempty + +getMotifGenome :: DiffFile + -> Maybe BackgroundDiffFile + -> OutputPath + -> MotifGenomeCommand + -> GenomeFile + -> TopN + -> Maybe Node + -> IO () +getMotifGenome diffFile bgDiffFile outPath mc gf topN node = TU.sh $ do + tmpDir <- TU.liftIO getTmpDir + tmpBed <- TU.mktempfile tmpDir "motif_input.bed" + mkTmpBed tmpBed diffFile topN node + + tmpBgBed <- TU.mktempfile tmpDir "motif_bg_input.bed" + maybe + (return ()) + (\x -> mkTmpBed tmpBgBed (DiffFile . unBackgroundDiffFile $ x) topN node) + bgDiffFile + + let cmd = if isJust bgDiffFile + then + T.pack $ + TP.printf + (unMotifGenomeCommand mc) + (TU.format TU.fp tmpBed) + (unGenomeFile gf) + (TU.format TU.fp . unOutputPath $ outPath) + (TU.format TU.fp tmpBgBed) + else + T.pack $ + TP.printf + (unMotifGenomeCommand mc) + (TU.format TU.fp tmpBed) + (unGenomeFile gf) (TU.format TU.fp . unOutputPath $ outPath) TU.stdout . TU.inshell cmd $ mempty diff --git a/src/TooManyCells/Motifs/Types.hs b/src/TooManyCells/Motifs/Types.hs --- a/src/TooManyCells/Motifs/Types.hs +++ b/src/TooManyCells/Motifs/Types.hs @@ -22,5 +22,6 @@ newtype OutputPath = OutputPath { unOutputPath :: TU.FilePath } newtype GenomeFile = GenomeFile { unGenomeFile :: T.Text } newtype MotifCommand = MotifCommand { unMotifCommand :: String } +newtype MotifGenomeCommand = MotifGenomeCommand { unMotifGenomeCommand :: String } newtype TopN = TopN { unTopN :: Int } newtype Node = Node { unNode :: Int } deriving (Eq, Ord) diff --git a/src/TooManyCells/Program/Motifs.hs b/src/TooManyCells/Program/Motifs.hs --- a/src/TooManyCells/Program/Motifs.hs +++ b/src/TooManyCells/Program/Motifs.hs @@ -32,6 +32,10 @@ . unHelpful . motifCommand $ opts + motifGenomeCommand' = fmap MotifGenomeCommand + . unHelpful + . motifGenomeCommand + $ opts diffFile' = DiffFile . TU.fromText . unHelpful . diffFile $ opts backgroundDiffFile' = fmap (BackgroundDiffFile . TU.fromText) . unHelpful @@ -53,12 +57,22 @@ $ (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 + case motifGenomeCommand' of + Nothing -> TU.liftIO + $ getMotif + diffFile' + backgroundDiffFile' + outPath + motifCommand' + genome' + topN' + node + (Just x) -> TU.liftIO + $ getMotifGenome + diffFile' + backgroundDiffFile' + outPath + x + genome' + topN' + node diff --git a/src/TooManyCells/Program/Options.hs b/src/TooManyCells/Program/Options.hs --- a/src/TooManyCells/Program/Options.hs +++ b/src/TooManyCells/Program/Options.hs @@ -194,7 +194,8 @@ | 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." + , motifCommand :: Maybe String "([meme %s -nmotifs 50 -oc %s] | STRING) The command to find motifs in a fasta file. Can be any command that will be run on each fasta file converted from the bed optionally per node, but the first \"%s\" must be the input file, the second \"%s\" is the output. An example of homer: `/path/to/findMotifs.pl %s fasta %s`. Uses meme by default." + , motifGenomeCommand :: Maybe String "([Nothing] | STRING) The command to find motifs from a bed file instead of --motif-command (replaces that argument), as in homer's findMotifsGenome.pl. Can be any command that will be run on each bed file optionally per node, but the first \"%s\" must be the input file, the second \"%s\" is the genome file, and the last is the output. An example of homer: `/path/to/findMotifsGenome.pl %s %s %s`." , topN :: Maybe Int "([100] | INT ) The top INT differentially expressed features." , output :: Maybe String "([out] | STRING) The folder containing output." } | MatrixOutput { matrixPath :: [String] "(PATH) The path to the input directory containing the matrix output of cellranger (cellranger < 3 (matrix.mtx, genes.tsv, and barcodes.tsv) or cellranger >= 3 (matrix.mtx.gz, features.tsv.gz, and barcodes.tsv.gz) or an input csv file containing feature row names and cell column names. scATAC-seq is supported if input file contains \"fragments\", ends with \".tsv.gz\" (such as \"fragments.tsv.gz\" or \"sample1_fragments.tsv.gz\"), and is in the SORTED (sort -k1,1 -k2,2n) 10x fragments format (see also --binwidth, --no-binarize). If given as a list (--matrix-path input1 --matrix-path input2 etc.) then will join all matrices together. Assumes the same number and order of features in each matrix, so only cells are added." diff --git a/too-many-cells.cabal b/too-many-cells.cabal --- a/too-many-cells.cabal +++ b/too-many-cells.cabal @@ -1,5 +1,5 @@ name: too-many-cells -version: 2.1.0.1 +version: 2.1.1.0 synopsis: Cluster single cells and analyze cell clade relationships. description: Different methods to cluster and analyze single cell data with diversity indices and differential expression. homepage: http://github.com/GregorySchwartz/too-many-cells#readme