diff --git a/bamstats.cabal b/bamstats.cabal
--- a/bamstats.cabal
+++ b/bamstats.cabal
@@ -1,5 +1,5 @@
 Name:          bamstats
-Version:       0.2
+Version:       0.3
 License:       GPL
 Cabal-Version: >= 1.6
 Build-Type:    Simple
diff --git a/src/Bam.hs b/src/Bam.hs
--- a/src/Bam.hs
+++ b/src/Bam.hs
@@ -8,7 +8,7 @@
 import System.Console.CmdArgs
 
 data Options = Stats { numrd :: Maybe Int, inputs :: [FilePath] }
-             | Hist { numrd :: Maybe Int, bins, maxdist :: Int, inputs :: [FilePath] }
+             | Hist { numrd :: Maybe Int, bins, maxdist :: Int, inputs :: [FilePath], plot :: Maybe String }
              | Quants { numrd :: Maybe Int, mindist :: Int, delta :: Double, inputs :: [FilePath] }
              | Dump { numrd :: Maybe Int, inputs :: [FilePath] }
              deriving (Data,Typeable,Read,Show)
@@ -21,6 +21,7 @@
             , bins  = 20      &= help "number of bins" &= typ "INT"
             , maxdist = 1000  &= help "max insert size to count" &= typ "INT"
             , inputs = []     &= args &= typ "BAM file(s)"
+            , plot = Nothing  &= help "gnuplot output, with optional gnuplot commands" &= opt ""
             } &= help "Collect a histogram of insert sizes"
 quants = Quants {
               numrd = Nothing &= help "max number of reads to include"
@@ -44,10 +45,11 @@
 main = do
   o <- cmdArgs $ modes [clfy,hst,quants,dump] 
        &= help "Extract information from BAM files" 
-       &= program "bam" &= summary "bam v0.0, ©2012 Ketil Malde"
+       &= program "bam" &= summary "bam v0.3, ©2012 Ketil Malde"
   let geninp f = (case numrd o of Just x -> take x; Nothing -> id) `fmap` readBams f
       genout = putStrLn . case o of Stats {} -> display . classify (CS 0 0 0 0 0)
-                                    Hist {}  -> display . classify (histgen o)
+                                    Hist {}  -> (case plot o of Nothing -> display
+                                                                Just str -> genplot str) . classify (histgen o)
                                     Quants {} -> showQuants . classify (quantgen o)
                                     Dump {}  -> display . classify (Bams [])
   mapM_ (\f -> genout =<< geninp f) $ inputs o
diff --git a/src/Bio/SamTools/Classify.hs b/src/Bio/SamTools/Classify.hs
--- a/src/Bio/SamTools/Classify.hs
+++ b/src/Bio/SamTools/Classify.hs
@@ -4,6 +4,7 @@
 
 module Bio.SamTools.Classify where
 
+import Prelude hiding (sum)
 import Bio.SamTools.Bam
 import Data.List (foldl', intercalate)
 import Data.Maybe (isNothing, fromJust)
@@ -27,14 +28,14 @@
   where t = total c
 
 showQuants :: Stats Hist -> String
-showQuants c = unlines $ map (intercalate "\t")
-  [  "Alignment" : header (innies c)
-  ,  "innies   " : quants t (innies c)
-  ,  "outies   " : quants t (outies c)
-  ,  "lefties  " : quants t (lefties c)
-  ,  "righties " : quants t (righties c)
+showQuants cs = unlines $ map (intercalate "\t")
+  [  "Alignment" : header (innies cs)
+  ,  "innies   " : quants t (innies cs)
+  ,  "outies   " : quants t (outies cs)
+  ,  "lefties  " : quants t (lefties cs)
+  ,  "righties " : quants t (righties cs)
   , ["Total reads: ", show t]]
-  where t = total c
+  where t = total cs
         percentiles = [0.05,0.25,0.5,0.75,0.95] :: [Double]
         quants tot h = printf "%14d" (hcount h)
               : printf "%5.2f%%" (200*fromIntegral (hcount h)/fromIntegral tot::Double)
@@ -46,12 +47,28 @@
         go (f:fs) sum prev ((b,c):rest) = 
           let next = sum+c in 
           if next >= f
-          then let b' = b - round (fromIntegral (next-f)/fromIntegral c*fromIntegral (b-fst prev))
+          then let b' = b - round (fromIntegral (next-f)/fromIntegral c*fromIntegral (b-fst prev)::Double)
                in printf "%7d" b' : (go fs sum prev ((b,c):rest))
           else     go (f:fs) next (b,c) rest
         go [] _ _ _ = []
         go _ _ _ [] = error "ran out of reads!?"
 
+genplot :: String -> Stats Hist -> String
+genplot cmds h = preamble ++ cmds ++ "\n" ++ plot
+  where preamble = unlines ["set style data boxes"
+                           ,"set style fill solid border 0"
+                           ,"set xlabel 'insert length'"
+                           ,"set ylabel 'read count'"]
+        plot = "plot '-' ti 'innies', '-' ti 'outies', '-' ti 'lefties', '-' ti 'righties'\n"
+                  ++concat [ plot1 $ buckets $ innies h
+                           , plot1 $ map (\(x,y)->(-x,-y)) $ buckets $ outies h
+                           , plot1 $ map (\(x,y)->(-x, y)) $ buckets $ lefties h
+                           , plot1 $ map (\(x,y)->( x,-y)) $ buckets $ righties h]
+        plot1 = unlines . go 0
+        go prev ((b,c):rest@(_:_)) = (show ((prev+b) `div` 2) ++ " " ++ show c) : go (b+1) rest
+        go _ [_] = ["e"] -- ignore last bucket, which is a catch-all
+        go _ [] = error "no data"
+        
 class Insertable x where
   insert :: Bam1 -> x -> x
   disp1  :: Int  -> x -> [String]
@@ -135,7 +152,8 @@
   insert b h = H (hcount h + 1) (go (fromIntegral $ fromJust $ insertSize b) (buckets h))
     where go x ((b1,v):bs@(_:_)) =
             if x > b1 then (b1,v):go x bs else inc (b1,v):bs
-          go x [(b1,v)] = [inc (b1,v)]
+          go _ [(b1,v)] = [inc (b1,v)]
+          go _ [] = error "this never happens"
           inc (x,y) = let y' = y+1 in y' `seq` (x,y')
   disp1 tot h = printf "%14d" (hcount h)
               : printf "%5.2f%%" (200*fromIntegral (hcount h)/fromIntegral tot::Double)