diff --git a/differential.cabal b/differential.cabal
--- a/differential.cabal
+++ b/differential.cabal
@@ -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
diff --git a/src/Differential.hs b/src/Differential.hs
--- a/src/Differential.hs
+++ b/src/Differential.hs
@@ -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
diff --git a/src/Plot.hs b/src/Plot.hs
--- a/src/Plot.hs
+++ b/src/Plot.hs
@@ -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")
diff --git a/src/Types.hs b/src/Types.hs
--- a/src/Types.hs
+++ b/src/Types.hs
@@ -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
diff --git a/src/Utility.hs b/src/Utility.hs
--- a/src/Utility.hs
+++ b/src/Utility.hs
@@ -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)