differential 0.1.2.1 → 0.2.0.0

raw patch · 5 files changed

+35/−12 lines, 5 files

Files

differential.cabal +1/−1
src/Differential.hs +3/−3
src/Plot.hs +21/−4
src/Types.hs +4/−0
src/Utility.hs +6/−4

differential.cabal view

@@ -1,5 +1,5 @@ name:                differential-version:             0.1.2.1+version:             0.2.0.0 synopsis:            Finds out whether an entity comes from different distributions (statuses). description:         Uses statistical tests to find whether an entity comes from one or two distributions. homepage:            http://github.com/GregorySchwartz/differential#readme

src/Differential.hs view

@@ -119,7 +119,7 @@ differentialMatrixObsRow :: [Int] -- ^ as                          -> [Int] -- ^ bs                          -> S.SpMatrix Double-                         -> [(Log2Diff, Maybe PValue, Maybe FDR)]+                         -> [(Log2Diff, Maybe PValue, Maybe FDR, Maybe QValue)] differentialMatrixObsRow as bs = differentialMatrixFeatRow as bs . S.transpose  -- | Get differentials between columns (observations) of select rows (features)@@ -127,10 +127,10 @@ differentialMatrixFeatRow :: [Int] -- ^ as                           -> [Int] -- ^ bs                           -> S.SpMatrix Double-                          -> [(Log2Diff, Maybe PValue, Maybe FDR)]+                          -> [(Log2Diff, Maybe PValue, Maybe FDR, Maybe QValue)] differentialMatrixFeatRow as bs = withFDR . fmap obsToDiff . S.toRowsL   where-    withFDR xs = zipWith (\(!l, !p) fdr -> (l, p, fdr)) xs+    withFDR xs = zipWith (\(!l, !p) (!fdr, !qval) -> (l, p, fdr, qval)) xs                . getFDR 0.05                . fmap snd                $ xs

src/Plot.hs view

@@ -54,9 +54,11 @@ -- to the feature while Status refers to the differential group. Choose whether -- to normalize by the maximum value for a name. Returns a list where the first -- entry is the plot and the second entry is a data frame of the values.-plotSingleDiff :: Bool -> [Entity] -> R.R s (R.SomeSEXP s)-plotSingleDiff normalizeBool vals = do+plotSingleDiff :: Bool -> Bool -> Bool -> [Entity] -> R.R s (R.SomeSEXP s)+plotSingleDiff normalizeBool violin noOutlier vals = do     let normalize = if normalizeBool then 1 else 0 :: Double+        violinFlag = if violin then 1 else 0 :: Double+        noOutlierFlag = if noOutlier then 1 else 0 :: Double         names = fmap (T.unpack . unName . _name) vals         statuses = fmap (T.unpack . unStatus . _status) vals         values = fmap _value vals@@ -73,9 +75,14 @@           df = ddply(df, "name", transform, value = value / max(value))         } +        colors = brewer.pal(9, "Set1")+        pal = colorRampPalette(colors)+        statuses = sort(unique(df$status))+        myColors = pal(length(statuses))+        names(myColors) = statuses+         p = ggplot(df, aes(x = name, y = value, fill = status)) +-                geom_violin(alpha = 0.5, draw_quantiles = c(0.25, 0.5, 0.75), scale = "width") +-                scale_fill_brewer(palette = "Set1") ++                scale_fill_manual(values=myColors) +                 xlab("Feature") +                 theme_classic() +                 theme( axis.text = element_text(color = "black")@@ -84,6 +91,16 @@                     , axis.ticks.y = element_line(color = "black")                     , axis.text.x = element_text(angle=315, hjust=0)                     )++        if(violinFlag_hs) {+          p = p + geom_violin(draw_quantiles = c(0.25, 0.5, 0.75), scale = "width")+        } else {+          if(noOutlierFlag_hs) {+            p = p + geom_boxplot(outlier.shape = NA)+          } else {+            p = p + geom_boxplot()+          }+        }          if(normalize_hs) {           p = p + ylab("Normalized abundance")

src/Types.hs view

@@ -40,6 +40,9 @@ newtype FDR = FDR     { unFDR :: Double     } deriving (Eq,Ord,Read,Show,Generic)+newtype QValue = QValue+    { unQValue :: Double+    } deriving (Eq,Ord,Read,Show,Generic) newtype Name       = Name { unName :: T.Text } deriving (Eq, Ord, Read, Show, A.ToJSON,Generic) newtype Id         = Id { unId :: T.Text } deriving (Eq, Ord, Read, Show, A.ToJSON,Generic) newtype Status = Status@@ -51,6 +54,7 @@  instance NFData Log2Diff instance NFData PValue+instance NFData QValue instance NFData FDR instance NFData Name instance NFData Id

src/Utility.hs view

@@ -65,11 +65,11 @@ toInt32 x = fromIntegral x :: Int32  -- | Convert p-values to FDR using the Benjamini-Hochberg procedure.-getFDR :: Double -> [Maybe PValue] -> [Maybe FDR]+getFDR :: Double -> [Maybe PValue] -> [(Maybe FDR, Maybe QValue)] getFDR alpha xs =   fmap snd     . sortBy (compare `on` fst)-    . fmap (\(!r, (!o, !p)) -> (o, getFDR r p))+    . fmap (\(!r, (!o, !p)) -> (o, (getFDRCritical r p, getQValue r p)))     . drop 1 -- Get rid of starting value     . scanl scanFunc (0, (-1, Nothing))     . sortBy (compare `on` snd)@@ -80,5 +80,7 @@     scanFunc (!r, (_, !prev)) x@(_, Nothing) = (r, x)     scanFunc (!r, (_, !prev)) !x =       if snd x == prev then (r, x) else (r + 1, x) -- Handle ties-    getFDR _ Nothing     = Nothing-    getFDR rank (Just _) = Just . FDR $ alpha * (rank / m)+    getFDRCritical _ Nothing     = Nothing+    getFDRCritical rank (Just _) = Just . FDR $ alpha * (rank / m)+    getQValue _ Nothing     = Nothing+    getQValue rank (Just (PValue p)) = Just . QValue $ p * (m / rank)