diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +The MIT License (MIT) + +Copyright (c) 2014 Kai Zhang + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,2 @@ +the-bioinformatician-s-toolkit +============================== diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/benchmarks/bench.hs b/benchmarks/bench.hs new file mode 100644 --- /dev/null +++ b/benchmarks/bench.hs @@ -0,0 +1,53 @@ +{-# LANGUAGE OverloadedStrings #-} + +import Criterion.Main +import System.Random +import qualified Data.ByteString.Char8 as B +import Data.Default.Class +import qualified Data.Conduit.List as CL +import Data.Conduit +import Control.Monad.Identity +import System.IO.Unsafe +import AI.Clustering.Hierarchical + +import Bio.Data.Fasta +import Bio.Motif +import Bio.Motif.Search +import Bio.Motif.Alignment +import Bio.Seq + +dna :: DNA Basic +dna = fromBS $ B.pack $ map f $ take 5000 $ randomRs (0, 3) (mkStdGen 2) + where + f :: Int -> Char + f x = case x of + 0 -> 'A' + 1 -> 'C' + 2 -> 'G' + 3 -> 'T' + _ -> undefined + +pwm :: PWM +pwm = toPWM [ "0.3 0.3 0.3 0.1" + , "0 0.5 0 0.5" + , "0.1 0.2 0.5 0.3" + , "0.1 0.1 0.1 0.7" + , "0 0 0 1" + , "0.25 0.25 0.25 0.25" + , "0.1 0.1 0.3 0.5" + , "0.25 0.25 0 0.5" + , "0.1 0.1 0.7 0.1" + , "0 0 0 1" + ] + +motifs :: [Motif] +motifs = unsafePerformIO $ readFasta' "data/motifs.fasta" + +main :: IO () +main = defaultMain + [ bench "motif score" $ nf (scores def pwm) dna + , bgroup "TFBS scanning" [ bench "Naive" $ nf (\x -> runIdentity $ findTFBSSlow def pwm x (0.6 * optimalScore def pwm) $$ CL.consume) dna + , bench "look ahead" $ nf (\x -> runIdentity $ findTFBS def pwm x (0.6 * optimalScore def pwm) $$ CL.consume) dna + ] + , bench "motif merge" $ nf (\x -> fmap show $ flatten $ buildTree x) motifs + ] diff --git a/bioinformatics-toolkit.cabal b/bioinformatics-toolkit.cabal new file mode 100644 --- /dev/null +++ b/bioinformatics-toolkit.cabal @@ -0,0 +1,163 @@ +-- Initial bioinformatics-toolkit.cabal generated by cabal init. For +-- further documentation, see http://haskell.org/cabal/users-guide/ + +name: bioinformatics-toolkit +version: 0.1.0 +synopsis: A collection of bioinformatics tools +description: A collection of bioinformatics tools +license: MIT +license-file: LICENSE +author: Kai Zhang +maintainer: kai@kzhang.org +copyright: (c) 2014 Kai Zhang +category: Bio +build-type: Simple +extra-source-files: README.md +cabal-version: >=1.10 +data-files: + data/*.fasta + +library + hs-source-dirs: src + ghc-options: -Wall + + exposed-modules: + Bio.ChIPSeq + Bio.ChIPSeq.FragLen + Bio.Data.Bam + Bio.Data.Bed + Bio.Data.BigWig + Bio.Data.Fasta + Bio.GO + Bio.GO.Parser + Bio.GO.GREAT + Bio.Motif + Bio.Motif.Alignment + Bio.Motif.Search + Bio.RealWorld.BioGRID + Bio.RealWorld.ENCODE + Bio.RealWorld.Ensembl + Bio.RealWorld.ID + Bio.RealWorld.UCSC + Bio.Seq + Bio.Seq.IO + Bio.Utils.Functions + Bio.Utils.Misc + Bio.Utils.Overlap + Bio.Utils.Types + + -- other-modules: + -- other-extensions: + build-depends: + base >=4.0 && <5.0 + , aeson + , aeson-pretty + , bbi + , binary + , bytestring >=0.10 + , bytestring-lexing >=0.5 + , colour + , conduit + , containers >=0.5 + , data-default-class + , double-conversion + , clustering + , http-conduit + , hexpat + , matrices >=0.4.3 + , mtl >=2.1.3.1 + , parallel >=3.2 + , primitive + , palette + , samtools + , split + , statistics >=0.13.2.1 + , text >=0.11 + , transformers >=0.3.0.0 + , unordered-containers >=0.2 + , vector + , vector-algorithms + , word8 + , IntervalMap >=0.3 + + default-language: Haskell2010 + +executable mkindex + hs-source-dirs: exe + main-is: mkindex.hs + build-depends: + base >=4.6 && <5.0 + , bioinformatics-toolkit + , shelly + , text + default-language: Haskell2010 + +executable viewSeq + hs-source-dirs: exe + main-is: viewSeq.hs + build-depends: + base >=4.6 && <5.0 + , bioinformatics-toolkit + , bytestring + default-language: Haskell2010 + +--executable mergeMotifs +-- hs-source-dirs: exe +-- main-is: mergeMotifs.hs +-- build-depends: +-- base >=4.6 && <5.0 +-- , bioinformatics-toolkit +-- , bytestring +-- , clustering +-- , data-default-class +-- , split +-- , optparse-applicative +-- , haskell-plot >=0.2.0.0 +-- , diagrams-cairo >=1.3 +-- , diagrams-lib >=1.3 +-- default-language: Haskell2010 + +benchmark bench + type: exitcode-stdio-1.0 + main-is: benchmarks/bench.hs + default-language: Haskell2010 + build-depends: + base >=4.6 && <5.0 + , bioinformatics-toolkit + , random + , criterion + , clustering + , bytestring + , data-default-class + , conduit + , mtl + default-language: Haskell2010 + +test-suite tests + type: exitcode-stdio-1.0 + hs-source-dirs: tests + main-is: test.hs + other-modules: + Tests.Bed + , Tests.ChIPSeq + , Tests.Motif + , Tests.Seq + + default-language: Haskell2010 + build-depends: + base + , bytestring + , random + , vector + , data-default-class + , tasty + , tasty-golden + , tasty-hunit + , bioinformatics-toolkit + , conduit + , unordered-containers + , mtl + +source-repository head + type: git + location: https://github.com/kaizhang/bioinformatics-toolkit.git diff --git a/data/motifs.fasta b/data/motifs.fasta new file mode 100644 --- /dev/null +++ b/data/motifs.fasta @@ -0,0 +1,1000 @@ +>H3k27me3_1 +0.606774 0.016294 0.13288 0.244052 +0.036584 0.220797 0.734391 0.008228 +0.007188 0.022713 0.970099 0 +0.880368 0.00655 0 0.113082 +0.022933 0.370323 0.588998 0.017746 +0.03691 0.064294 0.898796 0 +0.036733 0.847295 0.05167 0.064303 +0.035226 0.936624 0.009752 0.018398 +0.038192 0 0.961808 0 +0.86875 0.034126 0 0.097123 + +>H3k27me3_2 +0.027569 0.836023 0.098072 0.038336 +0.066084 0.796257 0.084066 0.053593 +0.053722 0.101794 0.788767 0.055717 +0.027115 0.851343 0.096392 0.02515 +0.057722 0.783978 0.10818 0.05012 +0.036805 0.14198 0.733304 0.087911 +0.015704 0.891019 0.070288 0.022988 +0.053043 0.784949 0.107395 0.054612 +0.050031 0.135567 0.737444 0.076958 + +>H3k27me3_3 +0.057438 0.817772 0.062863 0.061926 +0.106073 0.731281 0.120311 0.042334 +0.030535 0.018387 0.88613 0.064948 +0.030741 0.879362 0.049584 0.040314 +0.015252 0.788258 0.137458 0.059032 +0.10851 0.817518 0.033884 0.040088 +0.036468 0.01862 0.816303 0.128608 +0.009933 0.730805 0.207425 0.051837 +0.03466 0.800567 0.091489 0.073284 + +>H3k27me3_4 +0.282 0.275 0.302 0.141 +0.278 0.207 0.301 0.215 +0.646 0.04 0.091 0.223 +0.053 0.707 0.062 0.178 +0.09 0.09 0.693 0.127 +0.224 0.462 0.078 0.236 +0.178 0.053 0.57 0.199 +0.128 0.646 0.086 0.14 +0.04 0.186 0.608 0.166 +0.165 0.16 0.535 0.14 +0.228 0.24 0.279 0.253 +0.242 0.293 0.287 0.179 + +>H3k27me3_5 +0 0 0.951071 0.048929 +0 0 0 1 +0.627359 0.095216 0.049576 0.227849 +0 0.024042 0.580115 0.395842 +0.028143 0.971857 0 0 +0 0.014661 0 0.985339 +0.005613 0.019322 0.975064 0 + +>H3k27me3_6 +0.057628 0.768463 0.028808 0.145102 +0.004678 0.995322 0 0 +0.84383 0.054519 0.070009 0.031642 +0.013305 0.045173 0.836304 0.105218 +0.380657 0.550477 0.068865 0 +0.020839 0.91739 0.061771 0 +0.112091 0.040031 0.815495 0.032383 +0 0.682259 0.274559 0.043182 +0 0.058162 0.941838 0 + +>H3k27me3_7 +0.065186 0.381721 0.450664 0.102429 +0.177422 0.633336 0.1695 0.019743 +0.022156 0.928003 0.020696 0.029145 +0.00933 0.028055 0.914386 0.048229 +0.851568 0.036666 0.030642 0.081124 +0.195895 0.028873 0.743155 0.032077 +0.030455 0.668012 0.089362 0.212171 +0.05471 0.033974 0.90287 0.008447 +0 0.97092 0.018842 0.010238 +0.00204 0.035057 0.82321 0.139693 +0.189764 0.229928 0.536579 0.043729 + +>H3k27me3_8 +0.221868 0.141908 0 0.636224 +0.20352 0.751002 0.020317 0.025161 +0.201842 0.316531 0.384956 0.096671 +0.048677 0.037684 0.872111 0.041527 +0.165474 0.038256 0.665057 0.131213 +0.072201 0.842767 0.019817 0.065215 +0.146093 0.028365 0.825541 0 +0.058821 0.802287 0.041842 0.09705 +0.144684 0.067184 0.750231 0.037901 +0.08984 0.805021 0.040715 0.064423 +0.839481 0.054112 0.059255 0.047152 +0 0.934819 0.03832 0.026862 + +>H3k27me3_9 +0.167882 0.520832 0.275958 0.035328 +0.186393 0.813607 0 0 +0.314604 0 0.670072 0.015324 +0 0.926258 0.069624 0.004118 +0 0.920505 0.079495 0 +0 0.770725 0.040069 0.189206 +0.009384 0.020025 0.970591 0 +0.861156 0.05519 0 0.083654 +0 0.035565 0.964435 0 + +>H3k27me3_10 +0 0 1 0 +0.887436 0.050346 0.010708 0.05151 +0.045451 0.004784 0.935866 0.013899 +0.026369 0.089881 0.88375 0 +0.983318 0.016682 0 0 +0 0.25089 0 0.74911 +0 0.675374 0.324626 0 +0 0 0.932479 0.067521 +0 0.932663 0 0.067337 + +>H3k27me3_11 +0.224457 0.725389 0.036157 0.013997 +0.046447 0 0.953553 0 +0.063569 0.918074 0.016276 0.002081 +0.812151 0 0.15863 0.029219 +0 0.842258 0 0.157742 +0.022002 0.020292 0.957706 0 +0.081601 0.82946 0.049243 0.039696 +0.044825 0.029444 0.819662 0.10607 +0 0.64916 0.03712 0.31372 + +>H3k27me3_12 +0.001305 0.985088 0.006533 0.007073 +0.030098 0.137468 0.031795 0.800639 +0.000641 0.947329 0.050635 0.001395 +0.708868 0.138382 0.048309 0.104441 +0.01663 0.177837 0.764505 0.041027 +0.019051 0.285605 0.03219 0.663153 +0.018078 0.053027 0.222331 0.706565 +0.009817 0.183174 0.02938 0.77763 +0.003817 0.963035 0.008029 0.025119 +0.008196 0.944821 0.003097 0.043885 +0.190656 0.284665 0.035054 0.489625 +0.007074 0.235282 0.708158 0.049486 + +>H3k27me3_13 +0.201553 0.786786 0.011662 0 +0.056087 0.822486 0.109843 0.011584 +0.025109 0.017029 0.957862 0 +0.048154 0.648583 0.029026 0.274237 +0 0.034838 0.955107 0.010055 +0.887569 0 0 0.112431 +0.00911 0.869571 0.116133 0.005186 +0.224225 0.524249 0.027113 0.224412 +0.100124 0.825883 0.056214 0.017778 + +>H3k27me3_14 +0.041 0.778533 0.042471 0.137996 +0.123524 0.756722 0.072554 0.0472 +0.000288 0.917455 0.082257 0 +0.775207 0.042638 0.03771 0.144445 +0.006508 0.054321 0.897 0.042171 +0.070115 0.848451 0.062052 0.019382 +0.013144 0.870986 0.107971 0.007899 +0.147341 0.025824 0.057405 0.76943 +0.008468 0.74943 0.196232 0.04587 +0.078526 0.383205 0.380263 0.158006 +0.210026 0.546259 0.16247 0.081245 +0.13618 0.285236 0.265129 0.313454 + +>H3k27me3_15 +0 0.923859 0.049081 0.027059 +0.030773 0.054177 0.89204 0.02301 +0.03561 0.710155 0 0.254236 +0.028979 0.845767 0.059673 0.06558 +0.01451 0.083137 0.902353 0 +0.016964 0.039033 0.495903 0.448099 +0.153442 0 0.821859 0.024698 +0.003027 0.016753 0.895212 0.085008 +0.827733 0.10104 0.033935 0.037292 + +>H3k27me3_16 +0.102 0.672 0.131 0.095 +0.093 0.096 0.679 0.132 +0.077 0.744 0.095 0.084 +0.097 0.13 0.118 0.655 +0.061 0.728 0.117 0.094 +0.097 0.696 0.095 0.112 +0.096 0.117 0.694 0.093 +0.111 0.091 0.713 0.085 +0.109 0.087 0.725 0.079 +0.113 0.091 0.694 0.102 + +>H3k27me3_17 +0.001 0.898 0.1 0.001 +0.001 0.001 0.997 0.001 +0.219 0.001 0.543 0.237 +0.164 0.001 0.674 0.161 +0.098 0.001 0.9 0.001 +0.116 0.001 0.882 0.001 +0.001 0.097 0.901 0.001 +0.001 0.001 0.997 0.001 + +>H3k27me3_18 +0.030606 0.828389 0.056442 0.084563 +0.111882 0.773152 0.089181 0.025784 +0.094575 0.059243 0.746942 0.09924 +0 0.954034 0.025967 0.019999 +0.043145 0.011024 0.74285 0.202981 +0.042381 0.028283 0.791169 0.138166 +0.05721 0.033757 0.878238 0.030795 +0.727749 0.098964 0 0.173287 +0.125497 0.038022 0.795772 0.040709 + +>H3k27me3_19 +0.102 0.084 0.712 0.102 +0.104 0.675 0.115 0.106 +0.095 0.097 0.705 0.103 +0.107 0.111 0.672 0.11 +0.105 0.09 0.102 0.703 +0.089 0.105 0.111 0.695 +0.11 0.108 0.105 0.677 +0.711 0.088 0.105 0.096 +0.713 0.108 0.107 0.072 +0.738 0.082 0.094 0.086 + +>H3k27me3_20 +0.051681 0.823102 0.058038 0.06718 +0.04667 0.806892 0.04487 0.101568 +0.086563 0.079983 0.716515 0.116939 +0.003699 0.809887 0.142735 0.04368 +0.038984 0.031869 0.817261 0.111886 +0.049845 0.892505 0.051869 0.005781 +0.11927 0.072643 0.0486 0.759487 +0.010058 0.048014 0.901099 0.040829 +0.064414 0.721648 0.116811 0.097128 + +>H3k27me3_21 +0.012404 0.145428 0.788559 0.05361 +0 0.987277 0 0.012723 +0.020786 0.03082 0.877875 0.070519 +0.004679 0.016532 0.812235 0.166554 +0.052582 0.039503 0.671561 0.236354 +0.163953 0.656069 0.01867 0.161309 +0.059608 0.148713 0.084866 0.706813 +0 0.985746 0.010521 0.003733 +0 0.859543 0.092123 0.048335 + +>H3k27me3_22 +0 0 0 1 +0.012876 0 0.980217 0.006906 +0.03792 0.004296 0.957784 0 +0.053699 0 0.946301 0 +0.94609 0.013517 0 0.040393 +0.128065 0.022685 0.84925 0 +0.111168 0.053925 0.043403 0.791504 +0.473428 0.497362 0.019718 0.009492 + +>H3k27me3_23 +0.035442 0.762288 0.041676 0.160594 +0.109849 0.043451 0.732698 0.114002 +0.077134 0.893037 0.008684 0.021145 +0.185442 0.038754 0.747105 0.028699 +0.06961 0.026831 0.817424 0.086134 +0.019668 0.030807 0.911916 0.037608 +0.097454 0.082235 0.705757 0.114554 +0.8089 0.044458 0.086125 0.060517 +0.13045 0.857966 0.011584 0 + +>H3k27me3_24 +0.037785 0.189216 0 0.772999 +0 1 0 0 +0 1 0 0 +0.036078 0 0.886905 0.077017 +0 0.124555 0.83069 0.044755 +0.015533 0 0.984467 0 +0.121058 0 0.067534 0.811408 +0 0.058872 0.879147 0.06198 +0.009905 0.87634 0 0.113755 + +>H3k27me3_25 +0.761215 0.084224 0 0.154562 +0 0.095532 0.904468 0 +0.091023 0.047337 0.83473 0.02691 +0.136245 0.016155 0.726058 0.121543 +0.02227 0.026057 0.859822 0.091851 +0.206104 0.269121 0.521486 0.003288 +0.029464 0.841482 0.031084 0.097971 +0.026948 0 0.952093 0.020958 +0.016535 0.907432 0.049712 0.026321 + +>H3k27me3_26 +0.049624 0.873643 0.028015 0.048718 +0.252582 0.665923 0.025361 0.056133 +0.045616 0.891054 0.018 0.04533 +0.030306 0.009068 0.941617 0.019009 +0.008131 0.031025 0.950739 0.010106 +0.849258 0.045227 0.042674 0.062842 +0.03959 0.891805 0.053673 0.014932 +0.323839 0.048021 0.609607 0.018533 +0.006592 0.972602 0.005527 0.015279 + +>H3k27me3_27 +0.088867 0.8221 0.044076 0.044957 +0.042556 0.797011 0.143696 0.016737 +0.009207 0.88788 0.087999 0.014914 +0.781947 0.038703 0.126388 0.052963 +0.012945 0.030843 0.855373 0.100839 +0.034445 0.812728 0.144246 0.008581 +0.131431 0.047159 0.066305 0.755105 +0.00705 0.819166 0.161102 0.012682 +0.094177 0.720151 0.081899 0.103773 +0.184854 0.39469 0.329478 0.090978 +0.058616 0.370292 0.535453 0.035639 +0.006249 0.262678 0.726294 0.004778 + +>H3k27me3_28 +0.237059 0.085644 0 0.677297 +0.036429 0 0.74756 0.216011 +0 0.685895 0 0.314105 +0 0.990671 0.009329 0 +0.822357 0.086003 0.09164 0 +0 0.938432 0.061568 0 +0 0.557102 0 0.442898 +0 0.961506 0.038494 0 +0.916383 0 0.083617 0 + +>H3k36me3_1 +0.412511 0.030259 0.427783 0.129446 +0.916679 0.062226 0.01331 0.007785 +0.004271 0.008974 0.001146 0.985609 +0.00015 0.88472 0.003421 0.111709 +0.642004 0.002774 0.3484 0.006822 +0 0.983702 0 0.016298 +0.002246 0.232527 0.025492 0.739734 +0.001817 0.047181 0.00822 0.942782 +0.018288 0.012159 0.969052 0.000501 +0.972232 0.005142 0.012633 0.009994 + +>H3k36me3_2 +0.003553 0.962184 0.021214 0.013049 +0.010587 0.007471 0.007149 0.974793 +0.056903 0.002065 0.940403 0.00063 +0.029341 0.127911 0.012368 0.83038 +0.785189 0.080043 0.121928 0.01284 +0.714592 0.118011 0.145336 0.022061 +0.004274 0.118391 0.008479 0.868856 +0.008939 0.958376 0.01189 0.020794 +0.008619 0.934744 0.005097 0.05154 +0.164029 0.677445 0.093925 0.064601 + +>H3k36me3_3 +0.06728 0.003463 0.928563 0.000695 +0.042745 0.00829 0.002542 0.946423 +0.7652 0.083152 0.151537 0.000111 +0.752591 0.056491 0.174836 0.016082 +0.037203 0.041632 0.061134 0.860031 +0.009265 0.953082 0.018104 0.019549 +0.012137 0.928617 0.003188 0.056058 +0 0.953697 0.002099 0.044204 +0.988637 0.00156 0.008536 0.001267 + +>H3k36me3_4 +0.002921 0.975182 0.019283 0.002614 +0.389395 0.000673 0.019729 0.590203 +0.009176 0.003798 0.9865 0.000525 +0.0083 0.017494 0.966122 0.008083 +0.860109 0.077394 0.019347 0.04315 +0.024896 0.006935 0.96517 0.002999 +0.02811 0.022893 0.022859 0.926139 +0.031456 0.018898 0.528436 0.42121 +0.023575 0.914853 0.018795 0.042778 +0.456293 0 0.543707 0 + +>H3k36me3_5 +0.040974 0.015834 0.806165 0.137026 +0.000482 0.951344 0.005066 0.043108 +0.003276 0.952792 0.001466 0.042467 +0.015066 0.000847 0.008931 0.975156 +0.001057 0.962622 0.011316 0.025005 +0.008398 0.962897 0.007132 0.021573 +0.012183 0.890617 0.008478 0.088722 +0.732122 0.079602 0.16567 0.022606 +0.720023 0.13387 0.126859 0.019248 + +>H3k36me3_6 +0.000676 0.978498 0.018077 0.002749 +0.007447 0.021657 0.001438 0.969458 +0.02526 0.005255 0.963432 0.006053 +0.07272 0.003095 0.908564 0.015622 +0.024316 0.245887 0.724174 0.005622 +0.023618 0.808526 0.050988 0.116868 +0.780085 0.03254 0.176005 0.01137 +0.838104 0.073281 0.055882 0.032733 +0.018699 0.895298 0.064705 0.021298 +0.882149 0.01584 0.017592 0.084419 + +>H3k36me3_7 +0.932231 0.007543 0.057018 0.003208 +0.033304 0.038088 0.898609 0.029999 +0.020811 0.019055 0.953415 0.00672 +0.972527 0.020669 0.004139 0.002664 +0.020958 0.000775 0.978098 0.000169 +0.720805 0.014303 0.264037 0.000855 +0.974348 0.006919 0.010349 0.008383 +0.03354 0.013051 0.120669 0.832741 +0.006406 0.879915 0.00736 0.10632 +0.309651 0.049391 0.629367 0.01159 +0.038361 0.93892 0 0.022719 +0 0.004613 0.003955 0.991433 + +>H3k36me3_8 +0.008194 0.049311 0.015898 0.926597 +0.090861 0.564763 0.035339 0.309037 +0.178644 0.025419 0.794867 0.001069 +0.952999 0.030655 0.002213 0.014134 +0.006241 0.038922 0.951929 0.002908 +0.870211 0.060451 0.062974 0.006364 +0.01117 0.937917 0.005864 0.045049 +0.124574 0.853348 0.009747 0.012331 +0.925735 0.044778 0.020247 0.00924 + +>H3k36me3_9 +0.904352 0.033657 0.020888 0.041103 +0.025522 0.316953 0.654973 0.002552 +0.60872 0.051391 0.110131 0.229758 +0.02877 0.868779 0.089199 0.013253 +0.011526 0.964667 0.011639 0.012168 +0.933367 0.006451 0.006633 0.053549 +0.008692 0.020278 0.96526 0.00577 +0.008984 0.884634 0.094607 0.011775 +0.018634 0.916658 0.007392 0.057316 +0.172494 0.135223 0.23282 0.459462 + +>H3k36me3_10 +0.919685 0.012377 0.043492 0.024446 +0.07339 0.006966 0.901569 0.018075 +0.02117 0.041673 0.927641 0.009516 +0.039129 0.872396 0.042334 0.046141 +0.106186 0.120244 0.071998 0.701573 +0.038702 0.053091 0.907428 0.00078 +0.686984 0.061773 0.191688 0.059555 +0.212205 0.012475 0.759124 0.016196 +0.052953 0.044408 0.887356 0.015283 + +>H3k36me3_11 +0.919696 0.003455 0.049962 0.026887 +0.694257 0.264183 0.018576 0.022983 +0.009563 0.960945 0.002416 0.027076 +0.024854 0.004572 0.001935 0.968639 +0.01113 0.925409 0.014847 0.048614 +0.000416 0.958359 0.033243 0.007983 +0.046887 0.036202 0.008047 0.908863 +0.025582 0.00419 0.968692 0.001536 +0.813344 0.020974 0.142688 0.022994 + +>H3k36me3_12 +0.000105 0.991402 0.007214 0.001279 +0.00254 0.001259 0.000272 0.995929 +0.001064 0.992408 0.004293 0.002235 +0 0.010284 0.002637 0.987079 +0.949042 0.024663 0.017051 0.009244 +0.006653 0.987799 0.000557 0.00499 +0.230691 0.031043 0.008711 0.729555 +0.944702 0.008495 0.037851 0.008953 +0.90633 0.089644 0.001487 0.002539 +0.854696 0.130088 0.009795 0.005422 + +>H3k36me3_13 +0.00307 0.670061 0.000684 0.326185 +0.248509 0.001502 0.728177 0.021812 +0 0.985736 0 0.014264 +0 0 0.050282 0.949718 +0.000754 0.003429 0.011552 0.984265 +0.019913 0 0.979321 0.000766 +0.989083 0.004283 0.004944 0.00169 +0.992549 0.000812 0.005428 0.001212 +0 0.989639 0 0.010361 +0.004247 0.932986 0.048163 0.014604 + +>H3k36me3_14 +0.81174 0.019665 0.140907 0.027687 +0.080253 0.047986 0.097927 0.773834 +0.019734 0.914965 0.03224 0.033061 +0.045194 0.878126 0.026402 0.050277 +0.008454 0.956647 0.005975 0.028924 +0.928903 0.015369 0.02804 0.027687 +0.031443 0.073964 0.883553 0.01104 +0.006821 0.844347 0.128921 0.019911 +0.685602 0.075362 0.010203 0.228833 +0.162236 0.589332 0.228093 0.020339 + +>H3k36me3_15 +0.882402 0.033396 0.068885 0.015316 +0.830804 0.119149 0.04082 0.009227 +0.831424 0.037807 0.111478 0.019291 +0.048734 0.195949 0.754178 0.001138 +0.035274 0.016171 0.019563 0.928992 +0.037379 0.136887 0.820595 0.005138 +0.01301 0.923482 0.038584 0.024925 +0.046365 0.003483 0.001277 0.948875 +0.031845 0.001009 0.959658 0.007488 + +>H3k36me3_16 +0.22142 0.661973 0.052298 0.064309 +0.065963 0.723622 0.168499 0.041916 +0.672261 0.085026 0.135066 0.107646 +0.047246 0.166603 0.767594 0.018557 +0.016197 0.846993 0.012173 0.124637 +0.058715 0.772271 0.146818 0.022196 +0.062856 0.034337 0.136185 0.766622 +0.079586 0.143703 0.67515 0.101561 +0.047924 0.130594 0.775836 0.045647 + +>H3k36me3_17 +0.05465 0.033561 0.893461 0.018328 +0.084497 0.018699 0.88612 0.010684 +0.006624 0.016727 0.958797 0.017852 +0.007727 0.018333 0.002741 0.971199 +0.02039 0.170582 0.042949 0.76608 +0.003044 0.009809 0.00288 0.984266 +0.000433 0.920931 0.000922 0.077714 +0.673304 0.014996 0.159279 0.152421 +0.001107 0.990042 0.004468 0.004383 +0.019357 0.577936 0.006314 0.396393 + +>H3k36me3_18 +0 0.973926 0.007321 0.018753 +0.950074 0.005704 0.007072 0.03715 +0.010533 0.040159 0.928987 0.020321 +0.028786 0.880034 0.034524 0.056655 +0.017522 0.033617 0.016692 0.932169 +0.796171 0.092793 0.062744 0.048293 +0.069759 0.915717 0.001579 0.012946 +0.064312 0.01818 0.006552 0.910956 +0.013023 0.635436 0.094048 0.257493 +0.506109 0.152271 0.317617 0.024004 + +>H3k36me3_19 +0.009099 0.768966 0.174896 0.04704 +0.005736 0.949914 0.009358 0.034992 +0.008481 0.699832 0.004802 0.286885 +0.637203 0.048626 0.304202 0.00997 +0.700166 0.085486 0.129955 0.084392 +0.005993 0 0.994007 0 +0.022416 0.026348 0.030887 0.920349 +0.899538 0.006936 0.040875 0.052651 +0.034571 0.027553 0.922674 0.015202 +0.012281 0.974302 0.003413 0.010005 + +>H3k36me3_20 +0 0.954634 0.020821 0.024545 +0.959916 0.010502 0.01948 0.010102 +0 0.964407 0.004213 0.031381 +0.006762 0.006861 0.034083 0.952295 +0.001316 0.022362 0.188212 0.78811 +0.004588 0.039112 0.052845 0.903455 +0.097028 0.008342 0.887728 0.006902 +0.035378 0.005466 0.920855 0.038301 +0.062433 0.039808 0.892798 0.004961 + +>H3k36me3_21 +0.859214 0.02391 0.093131 0.023745 +0.011283 0.67469 0.030872 0.283154 +0.222336 0.011084 0.646888 0.119693 +0.018063 0.914944 0.018984 0.048009 +0.005282 0.966931 0.012938 0.01485 +0.018928 0.00973 0.005204 0.966138 +0.042841 0.025086 0.931121 0.000953 +0.078391 0.025731 0.053763 0.842115 +0.868978 0.057046 0.069564 0.004412 + +>H3k36me3_22 +0.60841 0.005091 0.2906 0.095899 +0.003456 0.99488 0.001023 0.00064 +0.004149 0.698563 0.081454 0.215835 +0.873008 0.0357 0.090623 0.00067 +0.138171 0.03005 0.009922 0.821856 +0.001478 0.068544 0.92731 0.002668 +0.025352 0 0.005404 0.969244 +0.16004 0.005446 0.031354 0.803159 +0.002049 0.037332 0.957441 0.003178 + +>H3k36me3_23 +0.016215 0.899832 0.033387 0.050566 +0.894176 0.061417 0.013931 0.030477 +0.000727 0.981877 0.014869 0.002526 +0.093005 0.057132 0.09865 0.751213 +0.026234 0.644096 0.256871 0.072798 +0.007867 0.906725 0.066711 0.018698 +0.755779 0.044818 0.031628 0.167774 +0.11542 0.054551 0.724506 0.105524 +0.014303 0.901416 0.055745 0.028536 + +>H3k36me3_24 +0.101078 0.638138 0.230415 0.03037 +0.04837 0.189042 0.164568 0.598021 +0.003459 0.929466 0.046241 0.020833 +0.909886 0.040184 0.017452 0.032478 +0 0.985452 0.010825 0.003723 +0.02308 0.035562 0.017328 0.924031 +0.026869 0.195428 0.681141 0.096562 +0.013393 0.791357 0.110883 0.084367 +0.852359 0.066776 0.030189 0.050677 + +>H3k36me3_25 +1 0 0 0 +0.997517 0 0.001911 0.000572 +0.989396 0.010604 0 0 +0.997988 0 0.002012 0 +0 0.004599 0.009508 0.985893 +0.002489 0.022044 0.007222 0.968244 +1 0 0 0 +0 0 1 0 + +>H3k36me3_26 +0 0.050941 0.9481 0.000958 +0.917975 0.004396 0.04399 0.033638 +0.080107 0.01739 0.009995 0.892508 +0.003122 0.910747 0.071006 0.015124 +0.856541 0.004123 0.092481 0.046855 +0 0.954192 0.021677 0.024131 +0.324062 0.003366 0.672572 0 +0.007931 0.891422 0.022519 0.078128 +0.019644 0.947779 0.008891 0.023686 +0.801951 0.055033 0.13526 0.007756 +0.064093 0.803326 0.002803 0.129778 + +>H3k36me3_27 +0.977359 0.010318 0.003964 0.008359 +0.003647 0.0648 0.922837 0.008716 +0.02984 0.134662 0.678923 0.156574 +0.057661 0.09916 0.035379 0.8078 +0.068183 0.827109 0.101843 0.002865 +0.944543 0.002283 0.03896 0.014213 +0.2072 0.0287 0.760046 0.004054 +0.072374 0.025213 0.887926 0.014487 +0.871529 0.014581 0.070648 0.043242 + +>H3k36me3_28 +0.001057 0.012283 0.98666 0 +0.080845 0.003927 0.914187 0.001042 +0.053964 0.89972 0.009015 0.037301 +0.584604 0.017961 0.397435 0 +0.000917 0.011319 0.987138 0.000625 +0.968632 0.013738 0.008291 0.009339 +0.006691 0 0.111601 0.881708 +0.001018 0.961972 0.025565 0.011445 +0.944747 0 0.055253 0 + +>H3k36me3_29 +0.919554 0.010718 0.050146 0.019582 +0.026146 0.922507 0.035952 0.015394 +0.038863 0.669532 0.038846 0.252759 +0.089452 0.798815 0.071686 0.040047 +0.003881 0.209842 0.009636 0.776641 +0.163106 0.038113 0.792181 0.0066 +0.008861 0.069439 0.010575 0.911125 +0.008895 0.891756 0.073583 0.025765 +0.042975 0.048174 0.02657 0.882281 + +>H3k36me3_30 +0.005893 0.014907 0.9792 0 +0.70184 0.008031 0 0.290129 +0.045003 0.012658 0.940072 0.002267 +0.015291 0.001339 0.983371 0 +0.023043 0.009685 0.004654 0.962617 +0.001328 0.011501 0.584003 0.403169 +0.009777 0.004999 0.978982 0.006243 +0.014365 0.929085 0.045343 0.011207 +0.965775 0.006529 0.009204 0.018492 +0 0.480552 0.457615 0.061833 + +>H3k36me3_31 +0.00587 0.212446 0.771873 0.009811 +0.013448 0.024481 0 0.962071 +0.018008 0.085861 0.223536 0.672595 +0.001198 0.066174 0.029182 0.903447 +0 0.928497 0.015111 0.056392 +0.685426 0.032575 0.123895 0.158104 +0.00504 0.949661 0.010899 0.0344 +0.000889 0.93481 0.004957 0.059345 +0.938088 0.016945 0.018351 0.026615 +0.002431 0.058992 0.241834 0.696743 + +>H3k36me3_32 +0.027825 0.835602 0.030067 0.106506 +0.147398 0.016565 0.015743 0.820294 +0.01865 0.017184 0.96346 0.000705 +0.811514 0.017567 0.103109 0.06781 +0.058709 0.010857 0.923336 0.007097 +0.22332 0.067862 0.702356 0.006463 +0.077166 0.812135 0.020988 0.08971 +0.794925 0.063402 0.106321 0.035352 +0.076781 0.017228 0.896096 0.009896 + +>H3k36me3_33 +0.28239 0.00308 0.66605 0.04848 +0.006395 0.980825 0.005253 0.007527 +0.943723 0.019739 0.007425 0.029113 +0.023457 0.026908 0.949635 0 +0.024511 0.045941 0.040224 0.889324 +0.087171 0.001414 0.90956 0.001856 +0.032953 0.086727 0.858562 0.021758 +0.043944 0.756823 0.004067 0.195166 +0.530697 0.017996 0.373231 0.078076 + +>H3k36me3_34 +1 0 0 0 +0.899924 0.012201 0.077813 0.010061 +0.983382 0.003567 0.010553 0.002498 +0.054735 0.942875 0.00239 0 +0.020183 0.708769 0.085405 0.185643 +0.010158 0.836989 0.004265 0.148587 +0.015848 0.933207 0.000393 0.050552 +0.834088 0.002349 0.160635 0.002928 +0 0.130753 0.006978 0.862269 + +>H3k36me3_35 +0.010712 0.847385 0.035623 0.10628 +0.003611 0.856391 0.005876 0.134122 +0.705357 0.008514 0.215873 0.070257 +0 0.076906 0.003708 0.919386 +0.010812 0.940153 0.010889 0.038146 +0.019558 0.011158 0.009881 0.959402 +0 0.912582 0.069699 0.017719 +0.205869 0.167232 0.034953 0.591946 +0.809279 0.042703 0.116952 0.031065 + +>H3k36me3_36 +0.774845 0.020892 0.188444 0.015819 +0.046869 0.046224 0.892087 0.014819 +0.038936 0.093175 0.711246 0.156643 +0.024784 0.926375 0.019708 0.029133 +0.845001 0.013071 0.106818 0.03511 +0.02403 0.044879 0.742407 0.188685 +0.021144 0.070082 0.903262 0.005513 +0.789184 0.064255 0.013028 0.133533 +0.057027 0.015332 0.894381 0.03326 + +>H3k36me3_37 +0.892006 0.050529 0.004518 0.052947 +0.929454 0.004419 0.02214 0.043986 +0.026154 0.889842 0.073229 0.010774 +0.916297 0.001224 0.062965 0.019515 +0.006583 0.293303 0.153373 0.546741 +0.794208 0.068959 0.115636 0.021197 +0.002683 0.179151 0.809093 0.009073 +0.053853 0.165226 0.049169 0.731753 +0.025423 0.002331 0.961049 0.011197 + +>H3k36me3_38 +0.611015 0.381214 0.004009 0.003763 +0.952794 0.025639 0.009368 0.012199 +0.153521 0.82543 0.012908 0.008141 +0.069069 0.848914 0.062295 0.019722 +0.021454 0.009984 0.010235 0.958328 +0.003052 0.943686 0.009798 0.043464 +0.015209 0.766826 0.152517 0.065448 +0.029932 0.046145 0.858404 0.065519 +0.006437 0.773855 0.007769 0.21194 + +>H3k36me3_39 +0 0.997159 0.001235 0.001606 +0.070598 0.018906 0 0.910496 +0.080071 0.670981 0.016896 0.232052 +0.847996 0.086118 0.004285 0.061601 +0.941722 0 0.058278 0 +0.945067 0 0.054022 0.000911 +0.947439 0.042437 0.006122 0.004001 +0.950897 0.042896 0.004833 0.001374 +0.972425 0.027575 0 0 +0.808207 0.189053 0.00274 0 + +>H3k36me3_40 +0 0.488302 0 0.511698 +0.326084 0.045974 0.555576 0.072367 +0 0.902679 0.003992 0.093329 +0.048088 0.945195 0.002591 0.004127 +0 0.948647 0.00164 0.049713 +0.930135 0.00143 0.056366 0.012069 +0 0.014649 0.981362 0.003989 +0.006494 0.989833 0.003673 0 +0 0.104576 0.001271 0.894152 +0.96452 0.0203 0.015181 0 + +>H3k36me3_41 +0.001894 0.950087 0.025304 0.022715 +0.921449 0.018594 0.03558 0.024377 +0.005648 0.941668 0.008425 0.044259 +0.00218 0.914476 0.023949 0.059394 +0.846917 0.039877 0.093493 0.019713 +0.02074 0.313675 0.06233 0.603255 +0.090929 0.028919 0.87049 0.009662 +0.020382 0.76811 0.015853 0.195655 +0.003603 0.882984 0.0047 0.108713 +0.11662 0.479464 0.102319 0.301596 + +>H3k36me3_42 +0.810977 0.004871 0.159193 0.024959 +0.011096 0.955056 0.019831 0.014016 +0.015746 0.961809 0.005249 0.017197 +0.076484 0.013566 0.007299 0.902651 +0 0.989289 0.010711 0 +0.886521 0.031408 0.079126 0.002945 +0.255718 0.06 0.597657 0.086624 +0.097263 0.042865 0.77098 0.088892 +0.026254 0.006217 0.021618 0.945911 + +>H3k36me3_43 +0.006143 0.003397 0.963095 0.027365 +0.064462 0.848995 0 0.086543 +0.274996 0.004109 0.202708 0.518186 +0.00106 0.965502 0.015939 0.017499 +0.834139 0.036895 0.126006 0.00296 +0.01903 0.222489 0.011339 0.747143 +0.196842 0.013587 0.781826 0.007745 +0.001173 0.95275 0.011071 0.035006 +0.006254 0.985028 0.004546 0.004172 + +>H3k36me3_44 +0.10144 0.56942 0.095451 0.233689 +0.11475 0.008777 0.699136 0.177337 +0.009882 0.005122 0.015821 0.969174 +0.013055 0.004043 0.972658 0.010244 +0.88118 0.040719 0.010247 0.067853 +0.063349 0.016347 0.798433 0.121871 +0.016712 0.94582 0.019492 0.017976 +0.009032 0.91158 0.004368 0.07502 +0.850769 0.129613 0.005472 0.014146 + +>H3k36me3_45 +0.025097 0.879699 0.024053 0.071151 +0.108123 0.03167 0.04076 0.819447 +0.01298 0.878373 0.090374 0.018274 +0.027724 0.038849 0.007367 0.92606 +0.013234 0.007289 0.965278 0.014199 +0.021564 0.046201 0.01581 0.916424 +0.042964 0.639999 0.120265 0.196773 +0.595976 0.200519 0.160376 0.043129 +0.008453 0.940496 0.012248 0.038804 +0.061462 0.542457 0.137743 0.258339 + +>H3k36me3_46 +0.114206 0.108627 0.775148 0.00202 +0.040201 0.849107 0.075996 0.034697 +0.891444 0.020683 0.048754 0.039119 +0.053378 0.047139 0.897656 0.001827 +0.091474 0.040686 0.866214 0.001626 +0.674014 0.107546 0.008192 0.210248 +0.166723 0.022014 0.806289 0.004974 +0.067126 0.021483 0.904298 0.007093 +0.804812 0.136065 0.020322 0.0388 + +>H3k36me3_47 +0.024889 0.105094 0.767292 0.102726 +0.912153 0.006795 0.081053 0 +0.013073 0.009563 0.940711 0.036653 +0.010903 0.935927 0.05317 0 +0 0.159601 0.006981 0.833418 +0.293496 0.120075 0.586429 0 +0.007238 0.023724 0.002171 0.966867 +0 0.021055 0.971398 0.007547 +0.589138 0.04026 0.035987 0.334615 +0.006491 0.013307 0 0.980203 + +>H3k36me3_48 +0.934702 0 0.062228 0.00307 +0.003373 0.970668 0.001714 0.024245 +0.356315 0.39856 0.245125 0 +0 0.003698 0 0.996302 +0.000898 0.003489 0.995613 0 +0.003142 0.029837 0.937759 0.029262 +0.014567 0.971985 0.007442 0.006006 +0.147584 0.026321 0.163884 0.66221 +0.655472 0.018438 0.083972 0.242118 + +>H3k36me3_49 +0.759851 0.087607 0.048312 0.10423 +0.024825 0.934356 0.030984 0.009836 +0.003299 0.980731 0.010745 0.005226 +0.842888 0.080928 0.066673 0.009511 +0.00237 0.978292 0.008187 0.011152 +0.85562 0.068684 0.020269 0.055427 +0.038686 0.766025 0.064353 0.130936 +0.136724 0.652591 0.163121 0.047564 +0.020079 0.82756 0.106334 0.046027 +0.339958 0.044072 0.449717 0.166253 +0.002456 0.023074 0.87749 0.096979 +0.003583 0.779471 0.203144 0.013801 + +>H3k36me3_50 +0.128161 0.129664 0.713792 0.028383 +0.004565 0.9862 0.001003 0.008232 +0.798998 0.081084 0.042152 0.077766 +0.004875 0.778633 0.205189 0.011303 +0.019553 0.939324 0.020115 0.021008 +0.868394 0.018873 0.066862 0.045871 +0.003594 0.874075 0.063006 0.059324 +0.008774 0.86022 0.035407 0.095598 +0.811466 0.007365 0.141121 0.040049 + +>H3k36me3_51 +0.737976 0.018171 0.130776 0.113077 +0.011037 0.085834 0.898712 0.004417 +0.028514 0.031026 0.073368 0.867093 +0.011541 0.027051 0.951167 0.010241 +0.842861 0.009248 0.085301 0.06259 +0.027145 0.063794 0.689289 0.219772 +0.020781 0.880483 0.088308 0.010428 +0.12984 0.012688 0.024012 0.83346 +0 0.012803 0.971691 0.015507 + +>H3k36me3_52 +0.015139 0.34128 0.00703 0.636552 +0.705418 0.180077 0.07348 0.041025 +0 0.978434 0.014328 0.007238 +0.176047 0.203605 0.000578 0.619769 +0.003157 0.895704 0.097982 0.003157 +0.879347 0.005983 0.037696 0.076973 +0.054011 0.024876 0.919087 0.002026 +0.082555 0.068701 0.827464 0.02128 +0.782412 0.02585 0.008528 0.18321 + +>H3k36me3_53 +0.101 0.587 0.141 0.171 +0.073 0.1 0.678 0.149 +0.348 0.2 0.102 0.35 +0.178 0.476 0.055 0.291 +0.042 0.226 0.705 0.027 +0.85 0.034 0.065 0.051 +0.025 0.65 0.115 0.21 +0.184 0.081 0.705 0.03 + +>H3k36me3_54 +0.135136 0.697312 0.074937 0.092615 +0 0.947702 0.044472 0.007827 +0.029985 0.931444 0.029358 0.009212 +0.628134 0.001826 0.01032 0.359719 +0 0.938682 0.030387 0.030932 +0.078092 0.814343 0.022242 0.085322 +0.014109 0.045762 0.005665 0.934464 +0 0.943892 0.056108 0 +0.010588 0.148248 0.784909 0.056254 +0.320235 0.016901 0.521988 0.140877 + +>H3k36me3_55 +0.706908 0.042083 0.075769 0.17524 +0.013156 0.001979 0.950923 0.033942 +0.00974 0.851683 0.135991 0.002586 +0.035698 0.955895 0.008407 0 +0.726227 0.007418 0.245694 0.020661 +0 0.067194 0.922738 0.010068 +0.081936 0.063635 0.834154 0.020275 +0.014378 0.785319 0.011176 0.189127 +0.944722 0.01517 0.010148 0.02996 + +>H3k36me3_56 +0.923371 0.040613 0.035004 0.001011 +0.834497 0.090343 0.04551 0.029651 +0.007527 0.030823 0.959929 0.001721 +0.125878 0 0.870903 0.003219 +0.007889 0.832338 0 0.159773 +0.025808 0.084089 0.857356 0.032747 +0.029025 0.001524 0.962293 0.007158 +0.014345 0.238691 0.738412 0.008552 +0.018055 0.714245 0.101775 0.165925 + +>H3k4me1_1 +0.021 0.018 0.002 0.959 +0.041 0.136 0.015 0.808 +0.245 0.001 0.014 0.74 +0.243 0.078 0.666 0.013 +0.142 0.835 0.002 0.021 +0.878 0.04 0.001 0.081 +0.256 0.026 0.021 0.697 +0.34 0.116 0.217 0.327 +0.452 0.062 0.114 0.373 +0.085 0.464 0.321 0.13 +0.863 0.072 0.002 0.063 +0.805 0.071 0.091 0.033 + +>H3k4me1_2 +0.126221 0.096725 0.652362 0.124692 +0.157713 0.76706 0.074905 0.000322 +0.123264 0.011594 0.004509 0.860634 +0.006538 0.043494 0.947226 0.002742 +0.850792 0.026025 0.048173 0.07501 +0.029324 0.225233 0.739268 0.006175 +0.048641 0.103353 0.059443 0.788562 +0.091813 0.767385 0.096649 0.044154 +0.823017 0.070473 0.006845 0.099666 + +>H3k4me1_3 +0.142656 0.046092 0.628207 0.183044 +0.077532 0.721399 0.068271 0.132798 +0.862808 0.039974 0.036854 0.060365 +0.01964 0.009924 0.018916 0.95152 +0.053006 0 0.025394 0.9216 +0.033778 0.023827 0.012824 0.929571 +0.123322 0.106461 0.047877 0.72234 +0.764536 0.032298 0.029472 0.173695 +0.754843 0.118183 0.056065 0.070909 diff --git a/exe/mkindex.hs b/exe/mkindex.hs new file mode 100644 --- /dev/null +++ b/exe/mkindex.hs @@ -0,0 +1,11 @@ +module Main where + +import Bio.Seq.IO (mkIndex) +import Shelly +import qualified Data.Text as T +import System.Environment + +main :: IO () +main = do [inDir, outF] <- getArgs + fls <- shelly . lsT . fromText . T.pack $ inDir + mkIndex (map T.unpack fls) outF diff --git a/exe/viewSeq.hs b/exe/viewSeq.hs new file mode 100644 --- /dev/null +++ b/exe/viewSeq.hs @@ -0,0 +1,15 @@ +module Main where + +import qualified Data.ByteString.Char8 as B +import Bio.Seq +import Bio.Seq.IO +import System.Environment + +main :: IO () +main = do + [fl, chr, start, end] <- getArgs + g <- pack fl + s <- getSeqs g [(B.pack chr, read start, read end)] :: IO [Either String (DNA IUPAC)] + case head s of + Left err -> error err + Right x -> B.putStrLn . toBS $ x diff --git a/src/Bio/ChIPSeq.hs b/src/Bio/ChIPSeq.hs new file mode 100644 --- /dev/null +++ b/src/Bio/ChIPSeq.hs @@ -0,0 +1,251 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE BangPatterns #-} +module Bio.ChIPSeq + ( rpkmBed + , rpkmSortedBed + , monoColonalize + , profiling + , profilingCoverage + , rpkmBam + , tagCountDistr + , peakCluster + ) where + +import Bio.SamTools.Bam +import qualified Bio.SamTools.BamIndex as BI +import Control.Monad (liftM, forM_, forM) +import Control.Monad.Primitive (PrimMonad) +import Control.Monad.Trans.Class (lift) +import Data.Conduit +import qualified Data.Conduit.List as CL +import Data.Function (on) +import qualified Data.Foldable as F +import qualified Data.HashMap.Strict as M +import qualified Data.IntervalMap as IM +import Data.Maybe (fromJust) +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Algorithms.Intro as I +import qualified Data.Vector.Generic as G +import qualified Data.Vector.Generic.Mutable as GM + +import Bio.Data.Bam +import Bio.Data.Bed + +-- | process a sorted BED stream, keep only mono-colonal tags +monoColonalize :: Monad m => Conduit BED m BED +monoColonalize = do + x <- await + F.forM_ x loop + where + loop prev = do + x <- await + case x of + Nothing -> yield prev + Just current -> + case () of + _ | compareBed prev current == GT -> + error $ "Bio.ChIPSeq.monoColonalize: Input is not sorted: " ++ show prev ++ " > " ++ show current + | chromStart prev == chromStart current && + chromEnd prev == chromEnd current && + chrom prev == chrom current && + _strand prev == _strand current -> loop prev + | otherwise -> yield prev >> loop current +{-# INLINE monoColonalize #-} + +-- | calculate RPKM on a set of unique regions. Regions (in bed format) would be kept in +-- memory but not tag file. +-- RPKM: Readcounts per kilobase per million reads. Only counts the starts of tags +rpkmBed :: (PrimMonad m, BEDLike b, G.Vector v Double) + => [b] -> Sink BED m (v Double) +rpkmBed regions = do + v <- lift $ do v' <- V.unsafeThaw . V.fromList . zip [0..] $ regions + I.sortBy (compareBed `on` snd) v' + V.unsafeFreeze v' + let (idx, sortedRegions) = V.unzip v + n = G.length idx + rc <- rpkmSortedBed $ Sorted sortedRegions + + lift $ do + result <- GM.new n + G.sequence_ . G.imap (\x i -> GM.unsafeWrite result i (rc U.! x)) $ idx + G.unsafeFreeze result +{-# INLINE rpkmBed #-} + +-- | calculate RPKM on a set of regions. Regions must be sorted. The Sorted data +-- type is used to remind users to sort their data. +rpkmSortedBed :: (PrimMonad m, BEDLike b, G.Vector v Double) + => Sorted (V.Vector b) -> Sink BED m (v Double) +rpkmSortedBed (Sorted regions) = do + vec <- lift $ GM.replicate l 0 + n <- CL.foldM (count vec) (0 :: Int) + let factor = fromIntegral n / 1e9 + lift $ liftM (G.imap (\i x -> x / factor / (fromIntegral . size) (regions V.! i))) + $ G.unsafeFreeze vec + where + count v nTags tag = do + let chr = chrom tag + p | _strand tag == Just True = chromStart tag + | _strand tag == Just False = chromEnd tag - 1 + | otherwise = error "Unkown strand" + xs = concat . snd . unzip $ + IM.containing (M.lookupDefault IM.empty chr intervalMap) p + addOne v xs + return $ succ nTags + + intervalMap = sortedBedToTree (++) . Sorted . G.toList . G.zip regions . + G.map return . G.enumFromN 0 $ l + addOne v' = mapM_ $ \x -> GM.unsafeRead v' x >>= GM.unsafeWrite v' x . (+1) + l = G.length regions +{-# INLINE rpkmSortedBed #-} + +-- | divide each region into consecutive bins, and count tags for each bin. The +-- total number of tags is also returned +profiling :: (PrimMonad m, G.Vector v Int, BEDLike b) + => Int -- ^ bin size + -> [b] -- ^ regions + -> Sink BED m ([v Int], Int) +profiling k beds = do + initRC <- lift $ forM beds $ \bed -> do + let start = chromStart bed + end = chromEnd bed + num = (end - start) `div` k + index i = (i - start) `div` k + v <- GM.replicate num 0 + return (v, index) + + sink 0 $ V.fromList initRC + where + sink !nTags vs = do + tag <- await + case tag of + Just (BED chr start end _ _ strand) -> do + let p | strand == Just True = start + | strand == Just False = end - 1 + | otherwise = error "profiling: unkown strand" + overlaps = concat . snd . unzip $ + IM.containing (M.lookupDefault IM.empty chr intervalMap) p + lift $ forM_ overlaps $ \x -> do + let (v, idxFn) = vs `G.unsafeIndex` x + i = idxFn p + GM.unsafeRead v i >>= GM.unsafeWrite v i . (+1) + sink (nTags+1) vs + + _ -> do rc <- lift $ mapM (G.unsafeFreeze . fst) $ G.toList vs + return (rc, nTags) + + intervalMap = bedToTree (++) $ zip beds $ map return [0..] +{-# INLINE profiling #-} + +-- | divide each region into consecutive bins, and count tags for each bin. The +-- total number of tags is also returned +profilingCoverage :: (PrimMonad m, G.Vector v Int, BEDLike b1, BEDLike b2) + => Int -- ^ bin size + -> [b1] -- ^ regions + -> Sink b2 m ([v Int], Int) +profilingCoverage k beds = do + initRC <- lift $ forM beds $ \bed -> do + let start = chromStart bed + end = chromEnd bed + num = (end - start) `div` k + index i = (i - start) `div` k + v <- GM.replicate num 0 + return (v, index) + + sink 0 $ V.fromList initRC + where + sink !nTags vs = do + tag <- await + case tag of + Just bed -> do + let chr = chrom bed + start = chromStart bed + end = chromEnd bed + overlaps = concat . snd . unzip $ IM.intersecting + (M.lookupDefault IM.empty chr intervalMap) $ IM.IntervalCO start end + lift $ forM_ overlaps $ \x -> do + let (v, idxFn) = vs `G.unsafeIndex` x + lo = idxFn start + hi = idxFn end + forM_ [lo..hi] $ \i -> + GM.unsafeRead v i >>= GM.unsafeWrite v i . (+1) + sink (nTags+1) vs + + _ -> do rc <- lift $ mapM (G.unsafeFreeze . fst) $ G.toList vs + return (rc, nTags) + + intervalMap = bedToTree (++) $ zip beds $ map return [0..] +{-# INLINE profilingCoverage #-} + + +-- | calculate RPKM using BAM file (*.bam) and its index file (*.bam.bai), using +-- constant space +rpkmBam :: BEDLike b => FilePath -> Conduit b IO Double +rpkmBam fl = do + nTags <- lift $ readBam fl $$ CL.foldM (\acc bam -> return $ + if isUnmap bam then acc else acc + 1) 0.0 + handle <- lift $ BI.open fl + conduit nTags handle + where + conduit n h = do + x <- await + case x of + Nothing -> lift $ BI.close h + Just bed -> do let chr = chrom bed + s = chromStart bed + e = chromEnd bed + rc <- lift $ viewBam h (chr, s, e) $$ readCount s e + yield $ rc * 1e9 / n / fromIntegral (e-s) + conduit n h + readCount l u = CL.foldM f 0.0 + where + f acc bam = do let p1 = fromIntegral . fromJust . position $ bam + rl = fromIntegral . fromJust . queryLength $ bam + p2 = p1 + rl - 1 + return $ if isReverse bam + then if l <= p2 && p2 < u then acc + 1 + else acc + else if l <= p1 && p1 < u then acc + 1 + else acc +{-# INLINE rpkmBam #-} + +tagCountDistr :: PrimMonad m => G.Vector v Int => Sink BED m (v Int) +tagCountDistr = loop M.empty + where + loop m = do + x <- await + case x of + Just (BED chr s e _ _ (Just str)) -> do + let p | str = s + | otherwise = 1 - e + case M.lookup chr m of + Just table -> loop $ M.insert chr (M.insertWith (+) p 1 table) m + _ -> loop $ M.insert chr (M.fromList [(p,1)]) m + _ -> lift $ do + vec <- GM.replicate 100 0 + F.forM_ m $ \table -> + F.forM_ table $ \v -> do + let i = min 99 v + GM.unsafeRead vec i >>= GM.unsafeWrite vec i . (+1) + G.unsafeFreeze vec +{-# INLINE tagCountDistr #-} + +-- | cluster peaks +peakCluster :: (BEDLike b, Monad m) + => [b] -- ^ peaks + -> Int -- ^ radius + -> Int -- ^ cutoff + -> Source m BED +peakCluster peaks r th = mergeBedWith mergeFn peaks' $= CL.filter g + where + peaks' = map f peaks + f b = let chr = chrom b + c = (chromStart b + chromEnd b) `div` 2 + in asBed chr (c-r) (c+r) + mergeFn xs = BED (chrom $ head xs) lo hi Nothing (Just $ fromIntegral $ length xs) Nothing + where + lo = minimum . map chromStart $ xs + hi = maximum . map chromEnd $ xs + g b = fromJust (_score b) >= fromIntegral th +{-# INLINE peakCluster #-} diff --git a/src/Bio/ChIPSeq/FragLen.hs b/src/Bio/ChIPSeq/FragLen.hs new file mode 100644 --- /dev/null +++ b/src/Bio/ChIPSeq/FragLen.hs @@ -0,0 +1,52 @@ +{-# LANGUAGE BangPatterns #-} + +module Bio.ChIPSeq.FragLen + ( fragLength + , naiveCCWithSmooth + ) where + +import qualified Data.ByteString.Char8 as B +import qualified Data.HashSet as S +import qualified Data.HashMap.Strict as M +import Data.List (foldl', maximumBy) +import Data.Ord (comparing) +import Bio.Data.Bed +import Control.Parallel.Strategies (parMap, rpar) + +-- | estimate fragment length for a ChIP-seq experiment +fragLength :: (Int, Int) -> [BED] -> Int +fragLength (start, end) beds = fst $ maximumBy (comparing snd) $ + naiveCCWithSmooth 4 beds [start, start+2 .. end] +{-# INLINE fragLength #-} + +-- sizeDistribution :: (Int, Int) -> [BED] -> (Int, Int) + +fromBED :: [BED] -> [(B.ByteString, (S.HashSet Int, S.HashSet Int))] +fromBED = map toSet . M.toList . M.fromListWith f . map parseLine + where + parseLine :: BED -> (B.ByteString, ([Int], [Int])) + parseLine x = case _strand x of + Just True -> (_chrom x, ([_chromStart x], [])) + Just False -> (_chrom x, ([], [_chromEnd x])) + _ -> error "Unknown Strand!" + f (a,b) (a',b') = (a++a', b++b') + toSet (chr, (forwd, rev)) = (chr, (S.fromList forwd, S.fromList rev)) +{-# INLINE fromBED #-} + +-- | fast relative cross-correlation with smoothing +apprxCorr :: S.HashSet Int -> S.HashSet Int -> Int -> Int -> Int +apprxCorr forwd rev smooth d = S.foldl' f 0 rev + where + f :: Int -> Int -> Int + f !acc r | any (`S.member` forwd) [r-d-smooth..r-d+smooth] = acc + 1 + | otherwise = acc +{-# INLINE apprxCorr #-} + +-- | calcuate cross corrlation with different shifts +naiveCCWithSmooth :: Int -> [BED] -> [Int] -> [(Int, Int)] +naiveCCWithSmooth smooth input range = f . fromBED $ input + where + cc :: [(B.ByteString, (S.HashSet Int, S.HashSet Int))] -> Int -> Int + cc xs d = foldl' (+) 0 . map ((\(forwd, rev) -> apprxCorr forwd rev smooth d).snd) $ xs + f xs = zip range $ parMap rpar (cc xs) range +{-# INLINE naiveCCWithSmooth #-} diff --git a/src/Bio/Data/Bam.hs b/src/Bio/Data/Bam.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Data/Bam.hs @@ -0,0 +1,63 @@ +-------------------------------------------------------------------------------- +-- | +-- Module : $Header$ +-- Copyright : (c) 2014 Kai Zhang +-- License : MIT + +-- Maintainer : kai@kzhang.org +-- Stability : experimental +-- Portability : portable + +-- functions for processing BED files +-------------------------------------------------------------------------------- + +module Bio.Data.Bam + ( readBam + , bamToBed + , viewBam + ) where + +import Bio.SamTools.Bam +import Bio.SamTools.BamIndex +import Control.Monad.Trans.Class (lift) +import qualified Data.ByteString.Char8 as B +import Data.Conduit (Source, Conduit, yield) +import qualified Data.Conduit.List as CL +import Data.Maybe (fromJust) +import Bio.Data.Bed + +readBam :: FilePath -> Source IO Bam1 +readBam fl = do handle <- lift $ openBamInFile fl + go handle + where + go h = do x <- lift $ get1 h + case x of + Nothing -> lift $ closeInHandle h + Just bam -> yield bam >> go h +{-# INLINE readBam #-} + +bamToBed :: Monad m => Conduit Bam1 m BED +bamToBed = CL.mapMaybe f + where + f bam =case targetName bam of + Just chr -> + let start = fromIntegral . fromJust . position $ bam + end = start + (fromIntegral . fromJust . queryLength) bam + nm = Just . queryName $ bam + strand = Just . not . isReverse $ bam + in Just $ BED chr start end nm Nothing strand + _ -> Nothing +{-# INLINE bamToBed #-} + +viewBam :: IdxHandle -> (B.ByteString, Int, Int) -> Source IO Bam1 +viewBam handle (chr, s, e) = case lookupTarget (idxHeader handle) chr of + Nothing -> return () + Just chrId -> do + q <- lift $ query handle chrId (fromIntegral s,fromIntegral e) + go q + where + go q' = do r <- lift $ next q' + case r of + Nothing -> return () + Just bam -> yield bam >> go q' +{-# INLINE viewBam #-} diff --git a/src/Bio/Data/Bed.hs b/src/Bio/Data/Bed.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Data/Bed.hs @@ -0,0 +1,418 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE FlexibleContexts #-} +-------------------------------------------------------------------------------- +-- | +-- Module : $Header$ +-- Copyright : (c) 2014 Kai Zhang +-- License : MIT + +-- Maintainer : kai@kzhang.org +-- Stability : experimental +-- Portability : portable + +-- functions for processing BED files +-------------------------------------------------------------------------------- + +module Bio.Data.Bed + ( BEDLike(..) + , BEDTree + , bedToTree + , sortedBedToTree + , splitBed + , splitBedBySize + , splitBedBySizeLeft + , Sorted(..) + , sortBed + , intersectBed + , intersectBedWith + , intersectSortedBed + , intersectSortedBedWith + , isOverlapped + , mergeBed + , mergeBedWith + , mergeSortedBed + , mergeSortedBedWith + -- * BED6 format + , BED(..) + -- * BED3 format + , BED3(..) + , fetchSeq + , fetchSeq' + , compareBed + ) where + +import Control.Arrow ((***)) +import Control.Monad.State.Strict +import qualified Data.ByteString.Char8 as B +import Data.Conduit +import qualified Data.Conduit.List as CL +import Data.Default.Class (Default(..)) +import Data.Function (on) +import qualified Data.Foldable as F +import qualified Data.HashMap.Strict as M +import qualified Data.IntervalMap.Strict as IM +import Data.List (groupBy) +import Data.Maybe (fromMaybe) +import qualified Data.Vector as V +import qualified Data.Vector.Algorithms.Intro as I +import System.IO + +import Bio.Seq +import Bio.Seq.IO +import Bio.Utils.Misc (binBySizeLeft, binBySize, bins, readInt, readDouble) + +-- | a class representing BED-like data, e.g., BED3, BED6 and BED12. BED format +-- uses 0-based index (see documentation). +class BEDLike b where + -- | construct bed record from chromsomoe, start location and end location + asBed :: B.ByteString -> Int -> Int -> b + + -- | convert bytestring to bed format + fromLine :: B.ByteString -> b + + -- | convert bed to bytestring + toLine :: b -> B.ByteString + + -- | field accessor + chrom :: b -> B.ByteString + chromStart :: b -> Int + chromEnd :: b -> Int + + toBed3 :: b -> BED3 + toBed3 bed = BED3 (chrom bed) (chromStart bed) (chromEnd bed) + + -- | return the size of a bed region + size :: b -> Int + size bed = chromEnd bed - chromStart bed + + hReadBed :: MonadIO m => Handle -> Source m b + hReadBed h = do eof <- liftIO $ hIsEOF h + unless eof $ do + line <- liftIO $ B.hGetLine h + yield $ fromLine line + hReadBed h + {-# INLINE hReadBed #-} + + -- | non-streaming version + hReadBed' :: MonadIO m => Handle -> m [b] + hReadBed' h = hReadBed h $$ CL.consume + {-# INLINE hReadBed' #-} + + readBed :: MonadIO m => FilePath -> Source m b + readBed fl = do handle <- liftIO $ openFile fl ReadMode + hReadBed handle + liftIO $ hClose handle + {-# INLINE readBed #-} + + -- | non-streaming version + readBed' :: MonadIO m => FilePath -> m [b] + readBed' fl = readBed fl $$ CL.consume + {-# INLINE readBed' #-} + + hWriteBed :: MonadIO m => Handle -> Sink b m () + hWriteBed handle = do + x <- await + case x of + Nothing -> return () + Just bed -> (liftIO . B.hPutStrLn handle . toLine) bed >> hWriteBed handle + {-# INLINE hWriteBed #-} + + hWriteBed' :: MonadIO m => Handle -> [b] -> m () + hWriteBed' handle beds = CL.sourceList beds $$ hWriteBed handle + {-# INLINE hWriteBed' #-} + + writeBed :: MonadIO m => FilePath -> Sink b m () + writeBed fl = do handle <- liftIO $ openFile fl WriteMode + hWriteBed handle + liftIO $ hClose handle + {-# INLINE writeBed #-} + + writeBed' :: MonadIO m => FilePath -> [b] -> m () + writeBed' fl beds = CL.sourceList beds $$ writeBed fl + {-# INLINE writeBed' #-} + + {-# MINIMAL asBed, fromLine, toLine, chrom, chromStart, chromEnd #-} + +type BEDTree a = M.HashMap B.ByteString (IM.IntervalMap Int a) + +-- | convert a set of bed records to interval tree, with combining function for +-- equal keys +sortedBedToTree :: (BEDLike b, F.Foldable f) + => (a -> a -> a) + -> Sorted (f (b, a)) + -> BEDTree a +sortedBedToTree f (Sorted xs) = + M.fromList + . map ((head *** IM.fromAscListWith f) . unzip) + . groupBy ((==) `on` fst) + . map (\(bed, x) -> (chrom bed, (IM.IntervalCO (chromStart bed) (chromEnd bed), x))) + . F.toList + $ xs +{-# INLINE sortedBedToTree #-} + +bedToTree :: BEDLike b + => (a -> a -> a) + -> [(b, a)] + -> BEDTree a +bedToTree f xs = + M.fromList + . map ((head *** IM.fromAscListWith f) . unzip) + . groupBy ((==) `on` fst) + . map (\(bed, x) -> (chrom bed, (IM.IntervalCO (chromStart bed) (chromEnd bed), x))) + . V.toList + $ xs' + where + xs' = V.create $ do + v <- V.unsafeThaw . V.fromList $ xs + I.sortBy (compareBed `on` fst) v + return v +{-# INLINE bedToTree #-} + +-- | split a bed region into k consecutive subregions, discarding leftovers +splitBed :: BEDLike b => Int -> b -> [b] +splitBed k bed = map (uncurry (asBed chr)) . bins k $ (s, e) + where + chr = chrom bed + s = chromStart bed + e = chromEnd bed +{-# INLINE splitBed #-} + +-- | split a bed region into consecutive fixed size subregions, discarding leftovers +splitBedBySize :: BEDLike b => Int -> b -> [b] +splitBedBySize k bed = map (uncurry (asBed chr)) . binBySize k $ (s, e) + where + chr = chrom bed + s = chromStart bed + e = chromEnd bed +{-# INLINE splitBedBySize #-} + +-- | split a bed region into consecutive fixed size subregions, including leftovers +splitBedBySizeLeft :: BEDLike b => Int -> b -> [b] +splitBedBySizeLeft k bed = map (uncurry (asBed chr)) . binBySizeLeft k $ (s, e) + where + chr = chrom bed + s = chromStart bed + e = chromEnd bed +{-# INLINE splitBedBySizeLeft #-} + +-- | a type to imply that underlying data structure is sorted +newtype Sorted b = Sorted {fromSorted :: b} + +compareBed :: (BEDLike b1, BEDLike b2) => b1 -> b2 -> Ordering +compareBed x y = compare x' y' + where + x' = (chrom x, chromStart x, chromEnd x) + y' = (chrom y, chromStart y, chromEnd y) +{-# INLINE compareBed #-} + +-- | sort BED, first by chromosome (alphabetical order), then by chromStart, last by chromEnd +sortBed :: BEDLike b => [b] -> Sorted (V.Vector b) +sortBed beds = Sorted $ V.create $ do + v <- V.unsafeThaw . V.fromList $ beds + I.sortBy compareBed v + return v +{-# INLINE sortBed #-} + +-- | return records in A that are overlapped with records in B +intersectBed :: (BEDLike b1, BEDLike b2) => [b1] -> [b2] -> [b1] +intersectBed a b = intersectSortedBed a b' + where + b' = sortBed b +{-# INLINE intersectBed #-} + +-- | return records in A that are overlapped with records in B +intersectSortedBed :: (BEDLike b1, BEDLike b2) + => [b1] -> Sorted (V.Vector b2) -> [b1] +intersectSortedBed a (Sorted b) = filter (not . null . f) a + where + f bed = let chr = chrom bed + interval = IM.IntervalCO (chromStart bed) $ chromEnd bed + in IM.intersecting (M.lookupDefault IM.empty chr tree) interval + tree = sortedBedToTree (\_ _ -> False) . Sorted . zip (V.toList b) . repeat $ False +{-# INLINE intersectSortedBed #-} + +intersectBedWith :: (BEDLike b1, BEDLike b2) + => ([b2] -> a) + -> [b1] + -> [b2] + -> [(b1, a)] +intersectBedWith fn a = intersectSortedBedWith fn a . sortBed +{-# INLINE intersectBedWith #-} + +intersectSortedBedWith :: (BEDLike b1, BEDLike b2) + => ([b2] -> a) + -> [b1] + -> Sorted (V.Vector b2) + -> [(b1, a)] +intersectSortedBedWith fn a (Sorted b) = map f a + where + f bed = let chr = chrom bed + interval = IM.IntervalCO (chromStart bed) $ chromEnd bed + in (bed, fn $ map snd $ IM.intersecting (M.lookupDefault IM.empty chr tree) interval) + tree = sortedBedToTree const . Sorted . V.toList . V.zip b $ b +{-# INLINE intersectSortedBedWith #-} + +isOverlapped :: (BEDLike b1, BEDLike b2) => b1 -> b2 -> Bool +isOverlapped bed1 bed2 = chr == chr' && not (e <= s' || e' <= s) + where + chr = chrom bed1 + s = chromStart bed1 + e = chromEnd bed1 + chr' = chrom bed2 + s' = chromStart bed2 + e' = chromEnd bed2 + +mergeBed :: (BEDLike b, Monad m) => [b] -> Source m b +mergeBed = mergeSortedBed . sortBed +{-# INLINE mergeBed #-} + +mergeBedWith :: (BEDLike b, Monad m) + => ([b] -> b) -> [b] -> Source m b +mergeBedWith f = mergeSortedBedWith f . sortBed +{-# INLINE mergeBedWith #-} + +mergeSortedBed :: (BEDLike b, Monad m) => Sorted (V.Vector b) -> Source m b +mergeSortedBed = mergeSortedBedWith f + where + f xs = asBed (chrom $ head xs) lo hi + where + lo = minimum . map chromStart $ xs + hi = maximum . map chromEnd $ xs +{-# INLINE mergeSortedBed #-} + +mergeSortedBedWith :: (BEDLike b, Monad m) + => ([b] -> b) -> Sorted (V.Vector b) -> Source m b +mergeSortedBedWith mergeFn (Sorted beds) = do + (_, r) <- V.foldM' f acc0 . V.tail $ beds + yield $ mergeFn r + where + x0 = V.head beds + acc0 = ((chrom x0, chromStart x0, chromEnd x0), [x0]) + f ((chr,lo,hi), acc) bed + | chr /= chr' || s' > hi = yield (mergeFn acc) >> + return ((chr',s',e'), [bed]) + | e' > hi = return ((chr',lo,e'), bed:acc) + | otherwise = return ((chr,lo,hi), bed:acc) + where + chr' = chrom bed + s' = chromStart bed + e' = chromEnd bed +{-# INLINE mergeSortedBedWith #-} + +-- * BED6 format + +-- | BED6 format, as described in http://genome.ucsc.edu/FAQ/FAQformat.html#format1.7 +data BED = BED + { _chrom :: !B.ByteString + , _chromStart :: {-# UNPACK #-} !Int + , _chromEnd :: {-# UNPACK #-} !Int + , _name :: !(Maybe B.ByteString) + , _score :: !(Maybe Double) + , _strand :: !(Maybe Bool) -- ^ True: "+", False: "-" + } deriving (Eq, Show) + +instance Default BED where + def = BED + { _chrom = "" + , _chromStart = 0 + , _chromEnd = 0 + , _name = Nothing + , _score = Nothing + , _strand = Nothing + } + +instance BEDLike BED where + asBed chr s e = BED chr s e Nothing Nothing Nothing + + fromLine l = evalState (f (B.split '\t' l)) 1 + where + f :: [B.ByteString] -> State Int BED + f [] = do i <- get + if i <= 3 then error "Read BED fail: Incorrect number of fields" + else return def + f (x:xs) = do + i <- get + put (i+1) + bed <- f xs + case i of + 1 -> return $ bed {_chrom = x} + 2 -> return $ bed {_chromStart = readInt x} + 3 -> return $ bed {_chromEnd = readInt x} + 4 -> return $ bed {_name = guard' x} + 5 -> return $ bed {_score = getScore x} + 6 -> return $ bed {_strand = getStrand x} + _ -> return def + + guard' x | x == "." = Nothing + | otherwise = Just x + getScore x | x == "." = Nothing + | otherwise = Just . readDouble $ x + getStrand str | str == "-" = Just False + | str == "+" = Just True + | otherwise = Nothing + {-# INLINE fromLine #-} + + toLine (BED f1 f2 f3 f4 f5 f6) = B.intercalate "\t" [ f1 + , (B.pack.show) f2 + , (B.pack.show) f3 + , fromMaybe "." f4 + , score' + , strand' + ] + where + strand' | f6 == Just True = "+" + | f6 == Just False = "-" + | otherwise = "." + score' = case f5 of + Just x -> (B.pack.show) x + _ -> "." + {-# INLINE toLine #-} + + chrom = _chrom + chromStart = _chromStart + chromEnd = _chromEnd + +-- | retreive sequences +fetchSeq :: BioSeq DNA a => Genome -> Conduit BED IO (Either String (DNA a)) +fetchSeq g = do gH <- lift $ gHOpen g + table <- lift $ readIndex gH + conduitWith gH table + lift $ gHClose gH + where + conduitWith h index' = do + bed <- await + case bed of + Just (BED chr start end _ _ isForward) -> do + dna <- lift $ getSeq h index' (chr, start, end) + case isForward of + Just False -> yield $ fmap rc dna + _ -> yield dna + conduitWith h index' + _ -> return () +{-# INLINE fetchSeq #-} + +fetchSeq' :: BioSeq DNA a => Genome -> [BED] -> IO [Either String (DNA a)] +fetchSeq' g beds = CL.sourceList beds $= fetchSeq g $$ CL.consume +{-# INLINE fetchSeq' #-} + +-- * BED3 format + +data BED3 = BED3 !B.ByteString !Int !Int deriving (Eq, Show) + +instance Default BED3 where + def = BED3 "" 0 0 + +instance BEDLike BED3 where + asBed = BED3 + + fromLine l = case B.split '\t' l of + (a:b:c:_) -> BED3 a (readInt b) $ readInt c + _ -> error "Read BED fail: Incorrect number of fields" + {-# INLINE fromLine #-} + + toLine (BED3 f1 f2 f3) = B.intercalate "\t" [f1, (B.pack.show) f2, (B.pack.show) f3] + {-# INLINE toLine #-} + + chrom (BED3 f1 _ _) = f1 + chromStart (BED3 _ f2 _) = f2 + chromEnd (BED3 _ _ f3) = f3 diff --git a/src/Bio/Data/BigWig.hs b/src/Bio/Data/BigWig.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Data/BigWig.hs @@ -0,0 +1,32 @@ +module Bio.Data.BigWig + ( module Data.BBI.BigWig + , accumScores + , normalizedScores + ) where + +import Data.BBI.BigWig +import Data.Conduit +import qualified Data.Conduit.List as CL +import Bio.Data.Bed + +-- | for each BED region, return total scores of all wig records overlapped with it +accumScores :: BEDLike b => BWFile -> Conduit b IO Double +accumScores bwF = CL.mapM (helper bwF) +{-# INLINE accumScores #-} + +-- | for each BED region, return normalized scores (divided by the length) of +-- all wig records overlapped with it +normalizedScores :: BEDLike b => BWFile -> Conduit b IO Double +normalizedScores bwF = CL.mapM $ \bed -> do + x <- helper bwF bed + return $ x / (fromIntegral . size) bed +{-# INLINE normalizedScores #-} + +helper :: BEDLike b => BWFile -> b -> IO Double +helper bwF bed = queryBWFile bwF (chr, start, end) $$ CL.fold f 0.0 + where + f acc (_, s, e, v) = acc + fromIntegral (e - s) * v + chr = chrom bed + start = chromStart bed + end = chromEnd bed +{-# INLINE helper #-} diff --git a/src/Bio/Data/Fasta.hs b/src/Bio/Data/Fasta.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Data/Fasta.hs @@ -0,0 +1,67 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +-------------------------------------------------------------------------------- +-- | +-- Module : $Header$ +-- Copyright : (c) 2014 Kai Zhang +-- License : MIT + +-- Maintainer : kai@kzhang.org +-- Stability : experimental +-- Portability : portable + +-- Functions for processing Fasta files +-------------------------------------------------------------------------------- + +module Bio.Data.Fasta + ( FastaLike(..) + , fastaReader + ) where + +import Bio.Motif +import Bio.Seq +import Control.Monad.State.Lazy +import qualified Data.ByteString.Char8 as B +import Data.Conduit +import qualified Data.Conduit.List as CL +import System.IO + +class FastaLike f where + fromFastaRecord :: ( B.ByteString -- ^ record header + , [B.ByteString] -- ^ record body + ) + -> f + + readFasta :: FilePath -> Source IO f + readFasta fl = mapOutput fromFastaRecord . fastaReader $ fl + + -- | non-stream version, read whole file in memory + readFasta' :: FilePath -> IO [f] + readFasta' fl = readFasta fl $$ CL.consume +-- writeFasta :: FilePath -> [f] -> IO () + {-# MINIMAL fromFastaRecord #-} + +instance BioSeq s a => FastaLike (s a) where + fromFastaRecord (_, xs) = fromBS . B.concat $ xs + {-# INLINE fromFastaRecord #-} + +instance FastaLike Motif where + fromFastaRecord (name, mat) = Motif name (toPWM mat) + {-# INLINE fromFastaRecord #-} + +fastaReader :: FilePath -> Source IO (B.ByteString, [B.ByteString]) +fastaReader fl = do handle <- liftIO $ openFile fl ReadMode + loop [] handle + where + loop :: [B.ByteString] -> Handle -> Source IO (B.ByteString, [B.ByteString]) + loop acc h = do + eof <- liftIO $ hIsEOF h + if eof then liftIO (hClose h) >> output acc else do + l <- liftIO $ B.hGetLine h + case () of + _ | B.null l -> loop acc h -- empty line, go to next line + | B.head l == '>' -> output acc >> loop [B.tail l] h + | otherwise -> loop (acc ++ [l]) h + output (x:xs) = yield (x, xs) + output _ = return () +{-# INLINE fastaReader #-} diff --git a/src/Bio/GO.hs b/src/Bio/GO.hs new file mode 100644 --- /dev/null +++ b/src/Bio/GO.hs @@ -0,0 +1,53 @@ +{-# LANGUAGE OverloadedStrings #-} +module Bio.GO + ( GO(..) + , GOId + , GOMap + , getParentById + ) where + +import qualified Data.Text as T +import Data.Text.Encoding (decodeUtf8) +import qualified Data.ByteString.Char8 as B +import Data.Foldable (foldl') +import qualified Data.HashMap.Strict as M +import Data.Maybe + +data GO = GO + { _oboId :: !GOId + , _label :: !T.Text + , _subProcessOf :: !(Maybe GOId) + , _oboNS :: !T.Text + } + +instance Show GO where + show (GO a b c d) = T.unpack $ T.intercalate "\t" + [decodeUtf8 a, b, decodeUtf8 $ fromMaybe "" c, d] + +type GOId = B.ByteString + +type GOMap = M.HashMap GOId GO + +getParentById :: GOId -> GOMap -> Maybe GO +getParentById gid goMap = M.lookup gid goMap >>= _subProcessOf + >>= (`M.lookup` goMap) +{-# INLINE getParentById #-} + +{- +buildGOTree :: GOMap -> [Tree GO] +buildGOTree goMap = map build roots + where + build = unfoldTree $ \x -> (x, M.lookupDefault [] (_oboId x) childrenMap) + + childrenMap = foldl' f M.empty goMap + where + f m go = case _subProcessOf go of + Just x -> M.insertWith (++) x [go] m + _ -> m + roots = mapMaybe f . M.keys $ childrenMap + where + f x = do g <- M.lookup x goMap + if isNothing $ _subProcessOf g + then Nothing + else Just g + -} diff --git a/src/Bio/GO/GREAT.hs b/src/Bio/GO/GREAT.hs new file mode 100644 --- /dev/null +++ b/src/Bio/GO/GREAT.hs @@ -0,0 +1,132 @@ +{-# LANGUAGE BangPatterns #-} +module Bio.GO.GREAT + ( AssocRule(..) + , getRegulatoryDomains + , enrichedTerms + ) where + +import Control.Monad.Primitive +import qualified Data.ByteString.Char8 as B +import Data.Conduit +import Data.Default.Class +import qualified Data.HashMap.Strict as M +import qualified Data.IntervalMap as IM +import qualified Data.Vector as V +import Data.Vector.Algorithms.Intro (sortBy) +import Data.List (sort, group) +import Data.Ord (comparing) +import Data.IntervalMap.Interval (lowerBound, upperBound) +import Data.Maybe + +import Bio.Data.Bed +import Bio.GO +import Bio.Utils.Functions + +-- | how to associate genomic regions with genes +data AssocRule = BasalPlusExtension Int Int Int + | TwoNearest + | OneNearest + +instance Default AssocRule where + def = BasalPlusExtension 5000 1000 50000 + +type Gene = ( B.ByteString -- ^ chromosome + , Int -- ^ tss + , Bool -- ^ is forward stranded + , [GOId] + ) + +-- | given a gene list and the rule, compute the rulatory domain for each gene +getRegulatoryDomains :: AssocRule -> [Gene] -> BEDTree [GOId] +getRegulatoryDomains (BasalPlusExtension up dw ext) gs = + bedToTree undefined $ flip map basal $ \(BED3 chr s e, go) -> + let intervals = fromJust . M.lookup chr $ basalTree + s' = case IM.intersecting intervals (IM.OpenInterval (s - ext) s) of + [] -> s - ext + x -> min s (maximum $ map (upperBound . fst) x) + e' = case IM.intersecting intervals (IM.OpenInterval e (e + ext)) of + [] -> e + ext + x -> max e (minimum $ map (lowerBound . fst) x) + in (BED3 chr s' e', go) + where + basalTree = bedToTree (error "encounter redundant regions") basal + basal = flip map gs $ \(chr,tss,str,go) -> + if str then (BED3 chr (tss - up) (tss + dw), go) + else (BED3 chr (tss - dw) (tss + up), go) +getRegulatoryDomains _ _ = undefined +{-# INLINE getRegulatoryDomains #-} + +-- | how many times a particular GO term is hit by given regions +countTerms :: (BEDLike b, Monad m) + => BEDTree [GOId] + -> Sink b m (Int, M.HashMap GOId Int) +countTerms tree = go 0 M.empty + where + go !n !termCount = do + x <- await + case x of + Just bed ->do + let chr = chrom bed + s = chromStart bed + e = chromEnd bed + terms = nub' . concatMap snd . IM.intersecting + (M.lookupDefault IM.empty chr tree) $ IM.IntervalCO s e + go (n+1) $ foldl (\acc t -> M.insertWith (+) t 1 acc) termCount terms + _ -> return (n, termCount) +{- + expand = nub' . foldl f [] + where + f acc i = traceBack i $ i : acc + traceBack i acc = case getParentById i goMap of + Just g -> traceBack (_oboId g) $ (_oboId g) : acc + _ -> acc +-} +{-# INLINE countTerms #-} + +nub' :: [B.ByteString] -> [B.ByteString] +nub' = map head . group . sort +{-# INLINE nub' #-} + +{- +-- | since GO terms are organized as a tree, a node is considered as being hit +-- if any of its children is hit. During this step, we traverse GO tree to include +-- parents if the regions that hit parents are different from those hiting their +-- children +expandTermCounts :: M.HashMap GOId Int -> GOMap -> M.HashMap GOId Int +expandTermCounts counts goMap = M.fromList . prune . expand $ counts + where + prune m = loop goTree + where + loop (Node p children : xs) = + let children' = filter (\(Node (_,c) _) -> c /= 0) children + l = length children' + in case () of + _ | l == 0 -> p : loop xs + | snd p == (snd . rootLabel . head) children' -> loop $ children' ++ xs + | otherwise -> p : (loop $ children' ++ xs) + loop _ = [] + + goTree :: [Tree (GOId, Int)] + goTree = map (fmap (\x -> (_oboId x, M.lookupDefault 0 (_oboId x) m))) . buildGOTree $ goMap + expand m = M.foldlWithKey' f M.empty m + where + f acc i c = traceBack c i $ M.insertWith (+) i c acc + traceBack c i m = case getParentById i goMap of + Just g -> traceBack c (_oboId g) $ M.insertWith (+) (_oboId g) c m + _ -> m +{-# INLINE expandTermCounts #-} +-} + +enrichedTerms :: (BEDLike b1, BEDLike b2, PrimMonad m) + => Source m b1 + -> Source m b2 + -> BEDTree [GOId] + -> m (V.Vector (GOId, (Double, Double))) +enrichedTerms fg bg tree = do + (n1, table1) <- fg $$ countTerms tree + (n2, table2) <- bg $$ countTerms tree + v <- V.unsafeThaw $ V.fromList $ M.toList $ flip M.mapWithKey table1 $ \t c -> + let k = fromMaybe (error "x") $ M.lookup t table2 + in (1 - hyperquick c k n1 n2, fromIntegral (c*n2) / fromIntegral (n1*k)) + sortBy (comparing snd) v + V.unsafeFreeze v diff --git a/src/Bio/GO/Parser.hs b/src/Bio/GO/Parser.hs new file mode 100644 --- /dev/null +++ b/src/Bio/GO/Parser.hs @@ -0,0 +1,37 @@ +{-# LANGUAGE OverloadedStrings #-} +module Bio.GO.Parser + ( readOWL + , readOWLAsMap + ) where + +import Control.Arrow ((&&&)) +import Control.Applicative ((<$>)) +import qualified Data.ByteString.Lazy.Char8 as L +import qualified Data.HashMap.Strict as M +import Data.Maybe +import qualified Data.Text as T +import Data.Text.Encoding (encodeUtf8) +import Text.XML.Expat.Proc +import Text.XML.Expat.Tree + +import Bio.GO + +readOWL :: FilePath -> IO [GO] +readOWL fl = do + c <- L.readFile fl + let (xml, _) = parse defaultParseOptions c + goTerms = findChildren "owl:Class" (xml :: Node T.Text T.Text) + return . map pickle $ goTerms + where + pickle x = + let id' = encodeUtf8 . f $ findChild "oboInOwl:id" x + label = f $ findChild "rdfs:label" x + parent = ( encodeUtf8 . T.replace "_" ":" . snd + . T.breakOnEnd "/" . snd . head . getAttributes + ) <$> findChild "rdfs:subClassOf" x + namespace = f $ findChild "oboInOwl:hasOBONamespace" x + in GO id' label parent namespace + f = getText . head . getChildren . fromJust + +readOWLAsMap :: FilePath -> IO GOMap +readOWLAsMap fl = M.fromList . map (_oboId &&& id) <$> readOWL fl diff --git a/src/Bio/Motif.hs b/src/Bio/Motif.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Motif.hs @@ -0,0 +1,373 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE BangPatterns #-} +module Bio.Motif + ( PWM(..) + , size + , subPWM + , rcPWM + , gcContentPWM + , Motif(..) + , Bkgd(..) + , toPWM + , ic + , scores + , scores' + , score + , optimalScore + , CDF + , cdf + , cdf' + , scoreCDF + , pValueToScore + , pValueToScoreExact + , toIUPAC + , readMEME + , toMEME + , fromMEME + , writeMEME + , writeFasta + + -- * References + -- $references + ) where + +import Prelude hiding (sum) +import Control.Arrow ((&&&)) +import Control.Monad.State.Lazy +import qualified Data.ByteString.Char8 as B +import Data.Conduit +import Data.Double.Conversion.ByteString (toShortest) +import Data.Default.Class +import Data.List (sortBy, foldl') +import Data.Maybe (fromJust, isNothing) +import Data.Ord (comparing) +import qualified Data.Vector as V +import qualified Data.Vector.Mutable as VM +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Algorithms.Intro as I +import qualified Data.Matrix.Unboxed as M +import Text.Printf (printf) + +import Bio.Seq +import Bio.Utils.Misc (readDouble, readInt) +import Bio.Utils.Functions (binarySearchBy) + +-- | k x 4 position weight matrix for motifs +data PWM = PWM + { _nSites :: !(Maybe Int) -- ^ number of sites used to generate this matrix + , _mat :: !(M.Matrix Double) + } deriving (Show, Read) + +size :: PWM -> Int +size (PWM _ mat) = M.rows mat + +-- | information content of a poistion in pwm +ic :: PWM -> Int -> Double +ic = undefined + +-- | extract sub-PWM given starting position and length, zero indexed +subPWM :: Int -> Int -> PWM -> PWM +subPWM i l (PWM n mat) = PWM n $ M.subMatrix (i,0) (i+l,3) mat +{-# INLINE subPWM #-} + +-- | reverse complementary of PWM +rcPWM :: PWM -> PWM +rcPWM (PWM n mat) = PWM n . M.fromVector d . U.reverse . M.flatten $ mat + where + d = M.dim mat +{-# INLINE rcPWM #-} + +-- | GC content of PWM +gcContentPWM :: PWM -> Double +gcContentPWM (PWM _ mat) = loop 0 0 / fromIntegral m + where + m = M.rows mat + loop !acc !i | i >= m = acc + | otherwise = + let acc' = acc + M.unsafeIndex mat (i,1) + M.unsafeIndex mat (i,2) + in loop acc' (i+1) + +data Motif = Motif + { _name :: !B.ByteString + , _pwm :: !PWM + } deriving (Show, Read) + +-- | background model which consists of single nucletide frequencies, and di-nucletide +-- frequencies +newtype Bkgd = BG (Double, Double, Double, Double) + +instance Default Bkgd where + def = BG (0.25, 0.25, 0.25, 0.25) + +-- | convert pwm to consensus sequence, see D. R. Cavener (1987). +toIUPAC :: PWM -> DNA IUPAC +toIUPAC (PWM _ pwm) = fromBS . B.pack . map f $ M.toRows pwm + where + f v | snd a > 0.5 && snd a > 2 * snd b = fst a + | snd a + snd b > 0.75 = iupac (fst a, fst b) + | otherwise = 'N' + where + [a, b, _, _] = sortBy (flip (comparing snd)) $ zip "ACGT" $ U.toList v + iupac x = case sort' x of + ('A', 'C') -> 'M' + ('G', 'T') -> 'K' + ('A', 'T') -> 'W' + ('C', 'G') -> 'S' + ('C', 'T') -> 'Y' + ('A', 'G') -> 'R' + _ -> undefined + sort' (x, y) | x > y = (y, x) + | otherwise = (x, y) + +-- | get scores of a long sequences at each position +scores :: Bkgd -> PWM -> DNA a -> [Double] +scores bg p@(PWM _ pwm) dna = go $! toBS dna + where + go s | B.length s >= len = scoreHelp bg p (B.take len s) : go (B.tail s) + | otherwise = [] + len = M.rows pwm +{-# INLINE scores #-} + +-- | a streaming version of scores +scores' :: Monad m => Bkgd -> PWM -> DNA a -> Source m Double +scores' bg p@(PWM _ pwm) dna = go 0 + where + go i | i < n - len + 1 = do yield $ scoreHelp bg p $ B.take len $ B.drop i s + go (i+1) + | otherwise = return () + s = toBS dna + n = B.length s + len = M.rows pwm +{-# INLINE scores' #-} + +score :: Bkgd -> PWM -> DNA a -> Double +score bg p dna = scoreHelp bg p $! toBS dna +{-# INLINE score #-} + +-- | the best possible score for a pwm +optimalScore :: Bkgd -> PWM -> Double +optimalScore (BG (a, c, g, t)) (PWM _ pwm) = foldl' (+) 0 . map f . M.toRows $ pwm + where f xs = let (i, s) = U.maximumBy (comparing snd) . + U.zip (U.fromList ([0..3] :: [Int])) $ xs + in case i of + 0 -> log $ s / a + 1 -> log $ s / c + 2 -> log $ s / g + 3 -> log $ s / t + _ -> undefined +{-# INLINE optimalScore #-} + +scoreHelp :: Bkgd -> PWM -> B.ByteString -> Double +scoreHelp (BG (a, c, g, t)) (PWM _ pwm) dna = fst . B.foldl f (0,0) $ dna + where + f (!acc,!i) x = (acc + sc, i+1) + where + sc = case x of + 'A' -> log $! matchA / a + 'C' -> log $! matchC / c + 'G' -> log $! matchG / g + 'T' -> log $! matchT / t + 'N' -> 0 + 'V' -> log $! (matchA + matchC + matchG) / (a + c + g) + 'H' -> log $! (matchA + matchC + matchT) / (a + c + t) + 'D' -> log $! (matchA + matchG + matchT) / (a + g + t) + 'B' -> log $! (matchC + matchG + matchT) / (c + g + t) + 'M' -> log $! (matchA + matchC) / (a + c) + 'K' -> log $! (matchG + matchT) / (g + t) + 'W' -> log $! (matchA + matchT) / (a + t) + 'S' -> log $! (matchC + matchG) / (c + g) + 'Y' -> log $! (matchC + matchT) / (c + t) + 'R' -> log $! (matchA + matchG) / (a + g) + _ -> error "Bio.Motif.score: invalid nucleotide" + matchA = addSome $ M.unsafeIndex pwm (i,0) + matchC = addSome $ M.unsafeIndex pwm (i,1) + matchG = addSome $ M.unsafeIndex pwm (i,2) + matchT = addSome $ M.unsafeIndex pwm (i,3) + addSome !y | y == 0 = pseudoCount + | otherwise = y + pseudoCount = 0.0001 +{-# INLINE scoreHelp #-} + +-- | the probability that a kmer is generated by background (0-orderd model) +pBkgd :: Bkgd -> B.ByteString -> Double +pBkgd (BG (a, c, g, t)) = B.foldl' f 1 + where + f acc x = acc * sc + where + sc = case x of + 'A' -> a + 'C' -> c + 'G' -> g + 'T' -> t + _ -> undefined +{-# INLINE pBkgd #-} + + +-- | unlike pValueToScore, this version compute the exact score but much slower +-- and is inpractical for long motifs +pValueToScoreExact :: Double -- ^ desirable p-Value + -> Bkgd + -> PWM + -> Double +pValueToScoreExact p bg pwm = go 0 0 . sort' . map ((scoreHelp bg pwm &&& pBkgd bg) . B.pack) . replicateM l $ "ACGT" + where + sort' xs = U.create $ do v <- U.unsafeThaw . U.fromList $ xs + I.sortBy (flip (comparing fst)) v + return v + go !acc !i vec | acc > p = fst $ vec U.! (i - 1) + | otherwise = go (acc + snd (vec U.! i)) (i+1) vec + l = size pwm + +-- | calculate the minimum motif mathching score that would produce a kmer with +-- p-Value less than the given number. This score would then be used to search +-- for motif occurrences with significant p-Value +pValueToScore :: Double -> Bkgd -> PWM -> Double +pValueToScore p bg pwm = cdf' (scoreCDF bg pwm) $ 1 - p + +newtype CDF = CDF (V.Vector (Double, Double)) + +-- P(X <= x) +cdf :: CDF -> Double -> Double +cdf (CDF v) x = let i = binarySearchBy cmp v x + in case () of + _ | i >= n -> snd $ V.last v + | i == 0 -> if x == fst (V.head v) then snd (V.head v) else 0.0 + | otherwise -> let (a, a') = v V.! (i-1) + (b, b') = v V.! i + α = (b - x) / (b - a) + β = (x - a) / (b - a) + in α * a' + β * b' + where + cmp (a,_) = compare a + n = V.length v +{-# INLINE cdf #-} + +-- the inverse of cdf +cdf' :: CDF -> Double -> Double +cdf' (CDF v) p + | p > 1 || p < 0 = error "p must be in [0,1]" + | otherwise = let i = binarySearchBy cmp v p + in case () of + _ | i >= n -> 1/0 + | i == 0 -> if p == snd (V.head v) then fst (V.head v) else undefined + | otherwise -> let (a, a') = v V.! (i-1) + (b, b') = v V.! i + α = (b' - p) / (b' - a') + β = (p - a') / (b' - a') + in α * a + β * b + where + cmp (_,a) = compare a + n = V.length v +{-# INLINE cdf' #-} + +-- approximate the cdf of motif matching scores +scoreCDF :: Bkgd -> PWM -> CDF +scoreCDF (BG (a,c,g,t)) pwm = toCDF $ loop (V.singleton 1, const 0) 0 + where + loop (prev,scFn) i + | i < n = + let (lo,hi) = minMax (1/0,-1/0) 0 + nBin' = min 100000 $ ceiling $ (hi - lo) / precision + step = (hi - lo) / fromIntegral nBin' + idx x = let j = truncate $ (x - lo) / step + in if j >= nBin' then nBin' - 1 else j + v = V.create $ do + new <- VM.replicate nBin' 0 + V.sequence_ $ flip V.imap prev $ \x p -> + when (p /= 0) $ do + let idx_a = idx $ sc + log' (M.unsafeIndex (_mat pwm) (i,0)) - log a + idx_c = idx $ sc + log' (M.unsafeIndex (_mat pwm) (i,1)) - log c + idx_g = idx $ sc + log' (M.unsafeIndex (_mat pwm) (i,2)) - log g + idx_t = idx $ sc + log' (M.unsafeIndex (_mat pwm) (i,3)) - log t + sc = scFn x + new `VM.read` idx_a >>= VM.write new idx_a . (a * p + ) + new `VM.read` idx_c >>= VM.write new idx_c . (c * p + ) + new `VM.read` idx_g >>= VM.write new idx_g . (g * p + ) + new `VM.read` idx_t >>= VM.write new idx_t . (t * p + ) + return new + in loop (v, \x -> (fromIntegral x + 0.5) * step + lo) (i+1) + | otherwise = (prev, scFn) + where + minMax (l,h) x + | x >= V.length prev = (l,h) + | prev V.! x /= 0 = + let sc = scFn x + s1 = sc + log' (M.unsafeIndex (_mat pwm) (i,0)) - log a + s2 = sc + log' (M.unsafeIndex (_mat pwm) (i,1)) - log c + s3 = sc + log' (M.unsafeIndex (_mat pwm) (i,2)) - log g + s4 = sc + log' (M.unsafeIndex (_mat pwm) (i,3)) - log t + in minMax (foldr min l [s1,s2,s3,s4],foldr max h [s1,s2,s3,s4]) (x+1) + | otherwise = minMax (l,h) (x+1) + toCDF (v, scFn) = CDF $ V.imap (\i x -> (scFn i, x)) $ V.scanl1 (+) v + precision = 0.002 + n = size pwm + log' x | x == 0 = log 0.001 + | otherwise = log x +{-# INLINE scoreCDF #-} + +-- | get pwm from a matrix +toPWM :: [B.ByteString] -> PWM +toPWM x = PWM Nothing . M.fromLists . map (map readDouble.B.words) $ x + +-- | pwm to bytestring +fromPWM :: PWM -> B.ByteString +fromPWM = B.unlines . map (B.unwords . map toShortest) . M.toLists . _mat + +writeFasta :: FilePath -> [Motif] -> IO () +writeFasta fl motifs = B.writeFile fl contents + where + contents = B.unlines . concatMap f $ motifs + f x = [">" `B.append` _name x, fromPWM $ _pwm x] + +readMEME :: FilePath -> IO [Motif] +readMEME = liftM fromMEME . B.readFile + +writeMEME :: FilePath -> [Motif] -> Bkgd -> IO () +writeMEME fl xs bg = B.writeFile fl $ toMEME xs bg + +toMEME :: [Motif] -> Bkgd -> B.ByteString +toMEME xs (BG (a,c,g,t)) = B.intercalate "" $ header : map f xs + where + header = B.pack $ printf "MEME version 4\n\nALPHABET= ACGT\n\nstrands: + -\n\nBackground letter frequencies\nA %f C %f G %f T %f\n\n" a c g t + f (Motif nm pwm) = + let x = "MOTIF " `B.append` nm + y = B.pack $ printf "letter-probability matrix: alength= 4 w= %d nsites= %d E= 0" (size pwm) sites + z = B.unlines . map (B.unwords . ("":) . map toShortest) . M.toLists . _mat $ pwm + sites | isNothing (_nSites pwm) = 0 + | otherwise = fromJust $ _nSites pwm + in B.unlines [x,y,z] +{-# INLINE toMEME #-} + +fromMEME :: B.ByteString -> [Motif] +fromMEME meme = evalState (go $ B.lines meme) (0, []) + where + go :: [B.ByteString] -> State (Int, [B.ByteString]) [Motif] + go (x:xs) + | "MOTIF" `B.isPrefixOf` x = put (1, [B.drop 6 x]) >> go xs + | otherwise = do + (st, str) <- get + case st of + 1 -> do when (startOfPwm x) $ put (2, str ++ [B.words x !! 7]) + go xs + 2 -> let x' = B.dropWhile (== ' ') x + in if B.null x' + then do put (0, []) + r <- go xs + return (toMotif str : r) + else put (2, str ++ [x']) >> go xs + _ -> go xs + go [] = do (st, str) <- get + return [toMotif str | st == 2] + startOfPwm = B.isPrefixOf "letter-probability matrix:" + toMotif (name:n:xs) = Motif name pwm + where + pwm = PWM (Just $ readInt n) $ M.fromLists . map (map readDouble.B.words) $ xs + toMotif _ = error "error" +{-# INLINE fromMEME #-} + +-- $references +-- +-- * Douglas R. Cavener. (1987) Comparison of the consensus sequence flanking +-- translational start sites in Drosophila and vertebrates. +-- /Nucleic Acids Research/ 15 (4): 1353–1361. +-- diff --git a/src/Bio/Motif/Alignment.hs b/src/Bio/Motif/Alignment.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Motif/Alignment.hs @@ -0,0 +1,116 @@ +{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE BangPatterns #-} + +module Bio.Motif.Alignment + ( alignment + , alignmentBy + , mergePWM + , buildTree + , progressiveMerging + ) where + +import AI.Clustering.Hierarchical +import qualified Data.Vector.Generic as G +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as U +import qualified Data.Matrix.Unboxed as M + +import Bio.Motif +import Bio.Utils.Functions + +-- | penalty function takes the gaps number as input, return penalty value +type PenalFn = Int -> Double + +type DistanceFn = (G.Vector v Double, G.Vector v (Double, Double)) => v Double -> v Double -> Double + +alignment :: PWM -> PWM -> (Double, (PWM, PWM, Int)) +alignment = alignmentBy jsd quadPenal + +-- | linear penalty +linPenal :: PenalFn +linPenal n = fromIntegral n * 0.3 + +-- | quadratic penalty +quadPenal :: PenalFn +quadPenal n = fromIntegral (n ^ (2 :: Int)) * 0.15 + +-- | cubic penalty +cubePenal :: PenalFn +cubePenal n = fromIntegral (n ^ (3 :: Int)) * 0.01 + +-- | exponentail penalty +expPenal :: PenalFn +expPenal n = fromIntegral (2^n :: Int) * 0.01 + +-- internal gaps are not allowed, larger score means larger distance, so the smaller the better +alignmentBy :: DistanceFn -> PenalFn -> PWM -> PWM -> (Double, (PWM, PWM, Int)) +alignmentBy fn pFn m1 m2 + | fst forwardAlign <= fst reverseAlign = (fst forwardAlign, (m1, m2, snd forwardAlign)) + | otherwise = (fst reverseAlign, (m1, m2', snd reverseAlign)) + where + forwardAlign | d1 < d2 = (d1,i1) + | otherwise = (d2,-i2) + where + (d1,i1) = loop opti2 (1/0,-1) s2 s1 0 + (d2,i2) = loop opti1 (1/0,-1) s1 s2 0 + reverseAlign | d1 < d2 = (d1,i1) + | otherwise = (d2,-i2) + where + (d1,i1) = loop opti2 (1/0,-1) s2' s1 0 + (d2,i2) = loop opti1 (1/0,-1) s1 s2' 0 + + loop opti (min',i') a b@(_:xs) !i + | currentBest >= min' = (min',i') + | d < min' = loop opti (d,i) a xs (i+1) + | otherwise = loop opti (min',i') a xs (i+1) + where + d = (G.sum sc + gapP) / fromIntegral (U.length sc + nGaps) + currentBest = opti U.! i + sc = U.fromList $ zipWith fn a b + nGaps = n1 + n2 - 2 * U.length sc + gapP = pFn nGaps + loop _ (min',i') _ _ _ = (min',i') + + opti1 = optimalSc n1 n2 + opti2 = optimalSc n2 n1 + + optimalSc x y = U.fromList $ scanr1 f $ go 0 + where + f v min' = min v min' + go i | a == 0 = [] + | otherwise = pFn b / fromIntegral (a + b) : go (i+1) + where + a = min x (y-i) + b = i + abs (x - (y-i)) + + s1 = M.toRows . _mat $ m1 + s2 = M.toRows . _mat $ m2 + s2' = M.toRows . _mat $ m2' + m2' = rcPWM m2 + n1 = length s1 + n2 = length s2 +{-# INLINE alignmentBy #-} + +mergePWM :: (PWM, PWM, Int) -> PWM +mergePWM (m1, m2, i) | i >= 0 = PWM Nothing (M.fromRows $ take i s1 ++ zipWith f (drop i s1) s2 ++ drop (n1 - i) s2) + | otherwise = PWM Nothing (M.fromRows $ take (-i) s2 ++ zipWith f (drop (-i) s2) s1 ++ drop (n2 + i) s1) + where + f = G.zipWith (\x y -> (x+y)/2) + s1 = M.toRows . _mat $ m1 + s2 = M.toRows . _mat $ m2 + n1 = length s1 + n2 = length s2 + +progressiveMerging :: Dendrogram Motif -> PWM +progressiveMerging t = case t of + Branch _ _ left right -> f (progressiveMerging left) $ progressiveMerging right + Leaf a -> _pwm a + where + f a b = mergePWM $! snd $ alignment a b + +-- | build a guide tree from a set of motifs +buildTree :: [Motif] -> Dendrogram Motif +buildTree motifs = hclust Average (V.fromList motifs) δ + where + δ (Motif _ x) (Motif _ y) = fst $ alignment x y diff --git a/src/Bio/Motif/Search.hs b/src/Bio/Motif/Search.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Motif/Search.hs @@ -0,0 +1,205 @@ +{-# LANGUAGE BangPatterns #-} + +module Bio.Motif.Search + ( findTFBS + , findTFBS' + , findTFBSSlow + , maxMatchSc + , MotifCompo(..) + , spacingConstraint + ) where + +import Bio.Seq (DNA, toBS) +import qualified Bio.Seq as Seq (length) +import Bio.Motif +import Control.Parallel.Strategies (parMap, rdeepseq) +import Data.Conduit +import qualified Data.HashSet as S +import Data.List (foldl', intercalate) +import Data.Ord (comparing) +import qualified Data.ByteString.Char8 as B +import qualified Data.Vector.Unboxed as U +import qualified Data.Matrix.Unboxed as M + +-- | given a user defined threshold, look for TF binding sites on a DNA +-- sequence, using look ahead search. This function doesn't search for binding +-- sites on the reverse strand +findTFBS :: Monad m => Bkgd -> PWM -> DNA a -> Double -> Source m Int +findTFBS bg pwm dna thres = loop 0 + where + loop !i | i >= l - n + 1 = return () + | otherwise = do let (d, _) = lookAheadSearch bg pwm sigma dna i thres + if d == n - 1 + then yield i >> loop (i+1) + else loop (i+1) + sigma = optimalScoresSuffix bg pwm + l = Seq.length dna + n = size pwm +{-# INLINE findTFBS #-} + +-- | a parallel version of findTFBS +findTFBS' :: Bkgd -> PWM -> DNA a -> Double -> [Int] +findTFBS' bg pwm dna th = concat $ parMap rdeepseq f [0,step..l-n+1] + where + f x = loop x + where + loop i | i >= x+step || i >= l-n+1 = [] + | otherwise = let d = fst $ lookAheadSearch bg pwm sigma dna i th + in if d == n-1 + then i : loop (i+1) + else loop (i+1) + sigma = optimalScoresSuffix bg pwm + l = Seq.length dna + n = size pwm + step = 500000 +{-# INLINE findTFBS' #-} + +-- | use naive search +findTFBSSlow :: Monad m => Bkgd -> PWM -> DNA a -> Double -> Source m Int +findTFBSSlow bg pwm dna thres = scores' bg pwm dna $= + loop 0 + where + loop i = do v <- await + case v of + Just v' -> if v' >= thres then yield i >> loop (i+1) + else loop (i+1) + _ -> return () +{-# INLINE findTFBSSlow #-} + +-- | the largest possible match scores starting from every position of a DNA sequence +maxMatchSc :: Bkgd -> PWM -> DNA a -> Double +maxMatchSc bg pwm dna = loop (-1/0) 0 + where + loop !max' !i | i >= l - n + 1 = max' + | otherwise = if d == n - 1 then loop sc (i+1) + else loop max' (i+1) + where + (d, sc) = lookAheadSearch bg pwm sigma dna i max' + sigma = optimalScoresSuffix bg pwm + l = Seq.length dna + n = size pwm +{-# INLINE maxMatchSc #-} + +optimalScoresSuffix :: Bkgd -> PWM -> U.Vector Double +optimalScoresSuffix (BG (a, c, g, t)) (PWM _ pwm) = + U.fromList . tail . map (last sigma -) $ sigma + where + sigma = scanl f 0 $ M.toRows pwm + f !acc xs = let (i, s) = U.maximumBy (comparing snd) . + U.zip (U.fromList ([0..3] :: [Int])) $ xs + in acc + case i of + 0 -> log $ s / a + 1 -> log $ s / c + 2 -> log $ s / g + 3 -> log $ s / t + _ -> undefined +{-# INLINE optimalScoresSuffix #-} + +lookAheadSearch :: Bkgd -- ^ background nucleotide distribution + -> PWM -- ^ pwm + -> U.Vector Double -- ^ best possible match score of suffixes + -> DNA a -- ^ DNA sequence + -> Int -- ^ starting location on the DNA + -> Double -- ^ threshold + -> (Int, Double) -- ^ (d, sc_d), the largest d such that sc_d > th_d +lookAheadSearch (BG (a, c, g, t)) pwm sigma dna start thres = loop (0, -1) 0 + where + loop (!acc, !th_d) !d + | acc < th_d = (d-2, acc) + | otherwise = if d >= n + then (d-1, acc) + else loop (acc + sc, thres - sigma U.! d) (d+1) + where + sc = case toBS dna `B.index` (start + d) of + 'A' -> log $! matchA / a + 'C' -> log $! matchC / c + 'G' -> log $! matchG / g + 'T' -> log $! matchT / t + 'N' -> 0 + 'V' -> log $! (matchA + matchC + matchG) / (a + c + g) + 'H' -> log $! (matchA + matchC + matchT) / (a + c + t) + 'D' -> log $! (matchA + matchG + matchT) / (a + g + t) + 'B' -> log $! (matchC + matchG + matchT) / (c + g + t) + 'M' -> log $! (matchA + matchC) / (a + c) + 'K' -> log $! (matchG + matchT) / (g + t) + 'W' -> log $! (matchA + matchT) / (a + t) + 'S' -> log $! (matchC + matchG) / (c + g) + 'Y' -> log $! (matchC + matchT) / (c + t) + 'R' -> log $! (matchA + matchG) / (a + g) + _ -> error "Bio.Motif.Search.lookAheadSearch: invalid nucleotide" + matchA = addSome $ M.unsafeIndex (_mat pwm) (d,0) + matchC = addSome $ M.unsafeIndex (_mat pwm) (d,1) + matchG = addSome $ M.unsafeIndex (_mat pwm) (d,2) + matchT = addSome $ M.unsafeIndex (_mat pwm) (d,3) + addSome !x | x == 0 = pseudoCount + | otherwise = x + pseudoCount = 0.0001 + n = size pwm +{-# INLINE lookAheadSearch #-} + +data MotifCompo = MotifCompo + { _motif1 :: Motif + , _motif2 :: Motif + , _sameDirection :: Bool + , _spacing :: Int + , _score :: Double + } + +instance Show MotifCompo where + show (MotifCompo m1 m2 isSame sp sc) = intercalate "\t" + [B.unpack $ _name m1, B.unpack $ _name m2, show isSame, show sp, show sc] + +-- | search for spacing constraint between two TFs +spacingConstraint :: Motif -- ^ motif 1 + -> Motif -- ^ motif 2 + -> Bkgd -- ^ backgroud nucleotide distribution + -> Double -- ^ p-Value for motif finding + -> Int -- ^ half window size + -> Int -- ^ max distance to search + -> DNA a -> [MotifCompo] +spacingConstraint m1@(Motif _ pwm1) m2@(Motif _ pwm2) bg th w k dna = same ++ oppose + where + rs = let rs' = [-k, -k+2*w+1 .. 0] + in rs' ++ map (*(-1)) (reverse rs') + -- on the same orientation + same = zipWith f rs $ zipWith (+) nFF nRR + where + nFF = map (nOverlap forward1 forward2 w) rs + nRR = map (nOverlap reverse1 reverse2 w) $ reverse rs + f r c = MotifCompo m1 m2 True r $ fromIntegral c / n + oppose = zipWith f rs $ zipWith (+) nFR nRF + where + nFR = map (nOverlap forward1 reverse2 w) rs + nRF = map (nOverlap reverse1 forward2 w) $ reverse rs + f r c = MotifCompo m1 m2 False r $ fromIntegral c / n + forward1 = findTFBS' bg pwm1 dna th1 + forward2 = S.fromList $ findTFBS' bg pwm2 dna th2 + reverse1 = map (+ (s1 - 1)) $ findTFBS' bg (rcPWM pwm1) dna th1' + reverse2 = S.fromList $ map (+ (s2 - 1)) $ findTFBS' bg (rcPWM pwm2) dna th2' + th1 = pValueToScore th bg pwm1 + th1' = pValueToScore th bg $ rcPWM pwm1 + th2 = pValueToScore th bg pwm2 + th2' = pValueToScore th bg $ rcPWM pwm2 + s1 = size pwm1 + s2 = size pwm2 + n = sqrt $ fromIntegral $ (length forward1 + length reverse1) * (S.size forward2 + S.size reverse2) + + nOverlap :: [Int] -> S.HashSet Int -> Int -> Int -> Int + nOverlap xs ys w' i = foldl' f 0 xs + where + f acc x | any (`S.member` ys) [x + i - w' .. x + i + w'] = acc + 1 + | otherwise = acc +{-# INLINE spacingConstraint #-} + +{- +computePValue :: Double -> [Int] -> [(Int, Double)] +computePValue p xs = zip xs $ map (pValue n p) xs + where + n = foldl' (+) 0 xs +{-# INLINE computePValue #-} + +pValue :: Int -> Double -> Int -> Double +pValue n p x | n > 2000 = complCumulative (poisson (fromIntegral n* p)) $ fromIntegral x + | otherwise = complCumulative (binomial n p) $ fromIntegral x +{-# INLINE pValue #-} +-} diff --git a/src/Bio/RealWorld/BioGRID.hs b/src/Bio/RealWorld/BioGRID.hs new file mode 100644 --- /dev/null +++ b/src/Bio/RealWorld/BioGRID.hs @@ -0,0 +1,91 @@ +{-# LANGUAGE OverloadedStrings #-} + +module Bio.RealWorld.BioGRID + ( TAB2(..) + , fetchByGeneNames + ) where + +import Network.HTTP.Conduit +import Data.List +import qualified Data.ByteString.Lazy.Char8 as BL +import qualified Data.ByteString.Char8 as B +import qualified Data.Text as T + +accessKey :: String +accessKey = "accessKey=6168b8d02b2aa2e9a45af6f3afac4461" + +base :: String +base = "http://webservice.thebiogrid.org/" + +-- | BioGRID tab2 format +data TAB2 = TAB2 + { _biogridId :: B.ByteString + , _entrezIdA :: B.ByteString + , _entrezIdB :: B.ByteString + , _biogridIdA :: B.ByteString + , _biogridIdB :: B.ByteString + , _systematicNameA :: T.Text + , _systematicNameB :: T.Text + , _symbolA :: T.Text + , _symbolB :: T.Text + , _synonymsA :: [T.Text] + , _synonymsB :: [T.Text] + , _experimentalSystemName :: T.Text + , _experimentalSystemType :: T.Text + , _firstAuthor :: T.Text + , _pubmedId :: B.ByteString + , _organismIdA :: B.ByteString + , _organismIdB :: B.ByteString + , _throughput :: T.Text + , _score :: Maybe Double + , _ptm :: T.Text + , _phenotypes :: [T.Text] + , _qualifications :: [T.Text] + , _tags :: [T.Text] + , _source :: T.Text + } deriving (Show) + +parseAsTab2 :: BL.ByteString -> TAB2 +parseAsTab2 l = TAB2 (BL.toStrict $ xs!!0) + (BL.toStrict $ xs!!1) + (BL.toStrict $ xs!!2) + (BL.toStrict $ xs!!3) + (BL.toStrict $ xs!!4) + (T.pack $ BL.unpack $ xs!!5) + (T.pack $ BL.unpack $ xs!!6) + (T.pack $ BL.unpack $ xs!!7) + (T.pack $ BL.unpack $ xs!!8) + (T.splitOn "|" $ T.pack $ BL.unpack $ xs!!9) + (T.splitOn "|" $ T.pack $ BL.unpack $ xs!!10) + (T.pack $ BL.unpack $ xs!!11) + (T.pack $ BL.unpack $ xs!!12) + (T.pack $ BL.unpack $ xs!!13) + (BL.toStrict $ xs!!14) + (BL.toStrict $ xs!!15) + (BL.toStrict $ xs!!16) + (T.pack $ BL.unpack $ xs!!17) + (getScore $ BL.unpack $ xs!!18) + (T.pack $ BL.unpack $ xs!!19) + (T.splitOn "|" $ T.pack $ BL.unpack $ xs!!20) + (T.splitOn "|" $ T.pack $ BL.unpack $ xs!!21) + (T.splitOn "|" $ T.pack $ BL.unpack $ xs!!22) + (T.pack $ BL.unpack $ xs!!23) + where + xs = BL.split '\t' l + getScore "-" = Nothing + getScore x = Just $ read x + +-- | retreive first 10,000 records +fetchByGeneNames :: [String] -> IO [TAB2] +fetchByGeneNames genes = do + initReq <- parseUrl $ intercalate "&" [url, geneList, tax, accessKey] + let request = initReq { method = "GET" + , requestHeaders = [("Content-type", "text/plain")] + } + r <- withManager $ \manager -> httpLbs request manager + return $ map parseAsTab2 $ BL.lines $ responseBody r + where + url = base ++ "/interactions/?searchNames=ture&includeInteractors=false" + geneList = "geneList=" ++ intercalate "|" genes + tax = "taxId=9606" +{-# INLINE fetchByGeneNames #-} diff --git a/src/Bio/RealWorld/ENCODE.hs b/src/Bio/RealWorld/ENCODE.hs new file mode 100644 --- /dev/null +++ b/src/Bio/RealWorld/ENCODE.hs @@ -0,0 +1,153 @@ +{-# LANGUAGE OverloadedStrings #-} +-------------------------------------------------------------------------------- +-- | +-- Module : $Header$ +-- Description : Search and download data from ENCODE project +-- Copyright : (c) Kai Zhang +-- License : MIT + +-- Maintainer : kai@kzhang.org +-- Stability : experimental +-- Portability : portable + +-- Search and download data from ENCODE project +-------------------------------------------------------------------------------- + +module Bio.RealWorld.ENCODE + ( KeyWords(..) + , search + + -- * common search + , terms + , cellIs + , organismIs + , assayIs + + -- * specific search + , getFile + , queryById + , openUrl + , jsonFromUrl + + -- * Inspection + , (|@) + , (|!) + , as + , (&) + + -- * printing + , showResult + ) where + +import Control.Applicative ((<$>)) +import Data.Aeson +import Data.Aeson.Types +import Data.Aeson.Encode.Pretty (encodePretty) +import qualified Data.HashMap.Lazy as M +import qualified Data.ByteString.Lazy.Char8 as B +import qualified Data.ByteString.Char8 as BS +import qualified Data.Sequence as S +import qualified Data.Foldable as F +import qualified Data.Text as T +import qualified Data.Vector as V +import Network.HTTP.Conduit +import Data.Default.Class +import Data.Monoid + +import Bio.RealWorld.ID + +data KeyWords = KeyWords (S.Seq String) -- ^ terms + (S.Seq String) -- ^ constraints + +instance Default KeyWords where + def = KeyWords S.empty $ S.fromList ["frame=object", "limit=all"] + +instance Show KeyWords where + show (KeyWords x y) = f x ++ g y + where + f x' | S.null x' = "" + | otherwise = "searchTerm=" ++ F.foldr1 (\a b -> b ++ ('+':a)) x' ++ "&" + g y' | S.null y' = "" + | otherwise = F.foldr1 (\a b -> b ++ ('&':a)) y' + +instance Monoid KeyWords where + mempty = KeyWords S.empty S.empty + mappend (KeyWords a b) (KeyWords a' b') = KeyWords (a S.>< a') (b S.>< b') + +base :: String +base = "https://www.encodeproject.org/" + +-- | general search using keywords and a set of constraints. Example: +-- search ["chip", "sp1"] ["type=experiment"] +search :: KeyWords -> IO (Either String [Value]) +search kw = do + initReq <- parseUrl url + let request = initReq { method = "GET" + , requestHeaders = [("accept", "application/json")] + } + r <- withManager $ \manager -> httpLbs request manager + return $ (eitherDecode . responseBody) r >>= + parseEither (withObject "ENCODE_JSON" (.: "@graph")) + where + url = base ++ "search/?" ++ show kw + +showResult :: Value -> IO () +showResult = B.putStrLn . encodePretty + +terms :: [String] -> KeyWords +terms xs = KeyWords (S.fromList xs) S.empty + +assayIs :: String -> KeyWords +assayIs x = KeyWords S.empty $ + S.fromList ["type=experiment", "assay_term_name=" ++ x] + +organismIs :: String -> KeyWords +organismIs x = KeyWords S.empty $ + S.fromList ["replicates.library.biosample.donor.organism.scientific_name=" ++ x] + +cellIs :: String -> KeyWords +cellIs x = KeyWords S.empty $ S.fromList ["biosample_term_name=" ++ x] + +-- | accession +queryById :: EncodeAcc -> IO (Either String Value) +queryById acc = jsonFromUrl $ "experiments/" ++ BS.unpack (fromID acc) + +getFile :: FilePath -> String -> IO () +getFile out url = openUrl (base ++ url) "application/octet-stream" >>= + B.writeFile out + +openUrl :: String -> String -> IO B.ByteString +openUrl url datatype = do + initReq <- parseUrl url + let request = initReq { method = "GET" + , requestHeaders = [("accept", BS.pack datatype)] + } + r <- withManager $ \manager -> httpLbs request manager + return $ responseBody r + +jsonFromUrl :: String -> IO (Either String Value) +jsonFromUrl url = eitherDecode <$> openUrl (base ++ url) "application/json" + + +(|@) :: Value -> T.Text -> Value +(|@) (Object obj) key = M.lookupDefault (error errMsg) key obj + where + errMsg = "No such key: " ++ T.unpack key ++ " In: " ++ show obj +(|@) _ _ = error "not an object" +{-# INLINE (|@) #-} + +(|!) :: Value -> Int -> Value +(|!) (Array ar) i = ar V.! i +(|!) _ _ = error "not an array" +{-# INLINE (|!) #-} + +(&) :: a -> (a -> b) -> b +(&) = flip ($) +{-# INLINE (&) #-} + +as :: FromJSON a => Value -> a +as = getResult . fromJSON + where + getResult (Error e) = error e + getResult (Success x) = x +{-# INLINE as #-} diff --git a/src/Bio/RealWorld/Ensembl.hs b/src/Bio/RealWorld/Ensembl.hs new file mode 100644 --- /dev/null +++ b/src/Bio/RealWorld/Ensembl.hs @@ -0,0 +1,41 @@ +{-# LANGUAGE OverloadedStrings #-} +module Bio.RealWorld.Ensembl + ( lookup + ) where + +import Prelude hiding (lookup) +import Data.Aeson +import Data.List.Split (chunksOf) +import qualified Data.ByteString.Char8 as B +import qualified Data.HashMap.Strict as M +import Network.HTTP.Conduit + +import Bio.RealWorld.ID (BioID(..), EnsemblID) + +base :: String +base = "http://rest.ensembl.org/" + +lookup :: [EnsemblID] -> IO (Either String Object) +lookup xs = do + rs <- mapM lookupHelp $ chunksOf 1000 xs + return $ foldl1 f rs + where + f a b = do + a' <- a + b' <- b + return $ M.union a' b' + +lookupHelp :: [EnsemblID] -> IO (Either String Object) +lookupHelp xs = do + initReq <- parseUrl url + let request = initReq { method = "POST" + , requestHeaders = [("Content-type", "application/json")] + , requestBody = body + } + r <- withManager $ \manager -> httpLbs request manager + return . eitherDecode . responseBody $ r + where + url = base ++ "/lookup/id/" + ids = B.pack $ show $ map fromID xs + body = RequestBodyBS $ B.intercalate "" ["{ \"ids\" :", ids, "}"] +{-# INLINE lookupHelp #-} diff --git a/src/Bio/RealWorld/ID.hs b/src/Bio/RealWorld/ID.hs new file mode 100644 --- /dev/null +++ b/src/Bio/RealWorld/ID.hs @@ -0,0 +1,27 @@ +module Bio.RealWorld.ID where + +import qualified Data.ByteString.Char8 as B + +class BioID a where + fromID :: a -> B.ByteString + toID :: B.ByteString -> a + +newtype UniprotID = UniprotID B.ByteString deriving (Show, Eq) + +newtype UCSCID = UCSCID B.ByteString deriving (Show, Eq) + +newtype GOID = GOID B.ByteString deriving (Show, Eq) + +-- | ENCODE Accession +newtype EncodeAcc = EncodeAcc B.ByteString deriving (Show, Eq) + +-- | Ensembl ID +newtype EnsemblID = EnsemblID B.ByteString deriving (Show, Eq) + +instance BioID EncodeAcc where + fromID (EncodeAcc x) = x + toID = EncodeAcc + +instance BioID EnsemblID where + fromID (EnsemblID x) = x + toID = EnsemblID diff --git a/src/Bio/RealWorld/UCSC.hs b/src/Bio/RealWorld/UCSC.hs new file mode 100644 --- /dev/null +++ b/src/Bio/RealWorld/UCSC.hs @@ -0,0 +1,87 @@ +{-# LANGUAGE OverloadedStrings #-} +-------------------------------------------------------------------------------- +-- | +-- Module : $Header$ +-- Description : Search and download data from ENCODE project +-- Copyright : (c) Kai Zhang +-- License : MIT + +-- Maintainer : kai@kzhang.org +-- Stability : experimental +-- Portability : portable + +-- resources from UCSC +-------------------------------------------------------------------------------- + +module Bio.RealWorld.UCSC + ( UCSCGene(..) + , getTSS + , getJunction + , readUCSCGenes + , readUCSCGenes' + ) where + +import Control.Monad.IO.Class (liftIO) +import qualified Data.ByteString.Char8 as B +import Data.Conduit +import qualified Data.Conduit.List as CL +import qualified Data.Vector.Unboxed as U +import System.IO + +import Bio.RealWorld.ID +import Bio.Utils.Misc (readInt) + +data UCSCGene = UCSCGene + { _geneName :: !B.ByteString + , _chrom :: !B.ByteString + , _strand :: !Bool + , _transcript :: !(Int, Int) + , _cds :: !(Int, Int) + , _exons :: !(U.Vector (Int, Int)) + , _introns :: !(U.Vector (Int, Int)) + , _proteinId :: !UniprotID + , _alignId :: !UCSCID + } deriving (Show) + +-- | get Transcription Start Site +getTSS :: UCSCGene -> (B.ByteString, Int) +getTSS g = (_chrom g, fst $ _transcript g) + +-- | get exon-intron junctions +getJunction :: UCSCGene -> (B.ByteString, U.Vector Int) +getJunction g = (_chrom g, U.map fst $ _introns g) + +-- | read genes from UCSC "knownGenes.tsv" +readUCSCGenes :: FilePath -> Source IO UCSCGene +readUCSCGenes fl = do + handle <- liftIO $ openFile fl ReadMode + _ <- liftIO $ B.hGetLine handle -- header + loop handle + where + loop h = do + eof <- liftIO $ hIsEOF h + if eof + then liftIO $ hClose h + else do + l <- liftIO $ B.hGetLine h + yield $ readGeneFromLine l + loop h +{-# INLINE readUCSCGenes #-} + +readUCSCGenes' :: FilePath -> IO [UCSCGene] +readUCSCGenes' fl = readUCSCGenes fl $$ CL.consume +{-# INLINE readUCSCGenes' #-} + +readGeneFromLine :: B.ByteString -> UCSCGene +readGeneFromLine xs = + let [f1,f2,f3,f4,f5,f6,f7,_,f9,f10,f11,f12] = B.split '\t' xs + str | f3 == "+" = True + | otherwise = False + trans = (readInt f4, readInt f5) + cds = (readInt f6, readInt f7) + exonStarts = map readInt . init . B.split ',' $ f9 + exonEnds = map readInt . init . B.split ',' $ f10 + exons = U.fromList $ zip exonStarts exonEnds + introns = U.fromList $ zip exonEnds $ tail exonStarts + in UCSCGene f1 f2 str trans cds exons introns (UniprotID f11) (UCSCID f12) +{-# INLINE readGeneFromLine #-} diff --git a/src/Bio/Seq.hs b/src/Bio/Seq.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Seq.hs @@ -0,0 +1,143 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE FlexibleContexts #-} +module Bio.Seq + ( + -- * Alphabet + Basic + , IUPAC + , Ext + -- * Sequence types + , DNA + , RNA + , Peptide + , BioSeq'(..) + , BioSeq(..) + -- * DNA related functions + , rc + , gcContent + , nucleotideFreq + ) where + +import Prelude hiding (length) +import qualified Data.ByteString.Char8 as B +import qualified Data.HashMap.Strict as M +import qualified Data.HashSet as S +import Data.Char8 (toUpper) +import Data.Monoid (Monoid(..)) + +-- | Alphabet defined by http://www.chem.qmul.ac.uk/iupac/ +-- | Standard unambiguous alphabet +data Basic + +-- | full IUPAC alphabet, including ambiguous letters +data IUPAC + +-- | extended alphabet +data Ext + +-- | DNA sequence +newtype DNA alphabet = DNA B.ByteString + +-- | RNA sequence +newtype RNA alphabet = RNA B.ByteString + +-- | Peptide sequence +newtype Peptide alphabet = Peptide B.ByteString + +instance Show (DNA a) where + show (DNA s) = B.unpack s + +instance Monoid (DNA a) where + mempty = DNA B.empty + mappend (DNA x) (DNA y) = DNA (x `mappend` y) + mconcat dnas = DNA . B.concat . map toBS $ dnas + +class BioSeq' s where + toBS :: s a -> B.ByteString + + slice :: Int -> Int -> s a -> s a + + length :: s a -> Int + length = B.length . toBS + {-# MINIMAL toBS, slice #-} + + +instance BioSeq' DNA where + toBS (DNA s) = s + slice i l (DNA s) = DNA . B.take l . B.drop i $ s + +instance BioSeq' RNA where + toBS (RNA s) = s + slice i l (RNA s) = RNA . B.take l . B.drop i $ s + +instance BioSeq' Peptide where + toBS (Peptide s) = s + slice i l (Peptide s) = Peptide . B.take l . B.drop i $ s + +class BioSeq' s => BioSeq s a where + alphabet :: s a -> S.HashSet Char + fromBS :: B.ByteString -> s a + +instance BioSeq DNA Basic where + alphabet _ = S.fromList "ACGT" + fromBS = DNA . B.map (f . toUpper) + where + f x | x `S.member` alphabet (undefined :: DNA Basic) = x + | otherwise = error $ "Bio.Seq.fromBS: unknown character: " ++ [x] + +instance BioSeq DNA IUPAC where + alphabet _ = S.fromList "ACGTNVHDBMKWSYR" + fromBS = DNA . B.map (f . toUpper) + where + f x | x `S.member` alphabet (undefined :: DNA IUPAC) = x + | otherwise = error $ "Bio.Seq.fromBS: unknown character: " ++ [x] + +instance BioSeq DNA Ext where + alphabet = undefined + fromBS = undefined + +instance BioSeq RNA Basic where + alphabet _ = S.fromList "ACGU" + fromBS = RNA . B.map (f . toUpper) + where + f x | x `S.member` alphabet (undefined :: RNA Basic) = x + | otherwise = error $ "Bio.Seq.fromBS: unknown character: " ++ [x] + +-- | O(n) Reverse complementary of DNA sequence +rc :: DNA alphabet -> DNA alphabet +rc (DNA s) = DNA . B.map f . B.reverse $ s + where + f x = case x of + 'A' -> 'T' + 'C' -> 'G' + 'G' -> 'C' + 'T' -> 'A' + _ -> x + +-- | O(n) Compute GC content +gcContent :: DNA alphabet -> Double +gcContent = (\(a,b) -> a / fromIntegral b) . B.foldl' f (0.0,0::Int) . toBS + where + f (!x,!n) c = + let x' = case c of + 'A' -> x + 'C' -> x + 1 + 'G' -> x + 1 + 'T' -> x + 'H' -> x + 0.25 + 'D' -> x + 0.25 + 'V' -> x + 0.75 + 'B' -> x + 0.75 + 'S' -> x + 1 + 'W' -> x + _ -> x + 0.5 -- "NMKYR" + in (x', n+1) + +-- | O(n) Compute single nucleotide frequency +nucleotideFreq :: BioSeq DNA a => DNA a -> M.HashMap Char Int +nucleotideFreq dna = B.foldl' f m0 . toBS $ dna + where + m0 = M.fromList . zip (S.toList $ alphabet dna) . repeat $ 0 + f m x = M.adjust (+1) x m +{-# INLINE nucleotideFreq #-} diff --git a/src/Bio/Seq/IO.hs b/src/Bio/Seq/IO.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Seq/IO.hs @@ -0,0 +1,128 @@ +{-# LANGUAGE BangPatterns #-} + +module Bio.Seq.IO + ( Genome + , GenomeH + , gHOpen + , gHClose + , pack + , getSeqs + , getSeq + , readIndex + , getChrom + , getChrSizes + , mkIndex + ) where + +import Bio.Seq +import Bio.Utils.Misc (readInt) +import Control.Applicative ((<$>)) +import qualified Data.ByteString.Char8 as B +import qualified Data.HashMap.Lazy as M +import Data.List.Split +import System.IO + +-- | The first 2048 bytes are header. Header consists of a magic string, +-- followed by chromosome information. Example: +-- \nCHR1 START SIZE +newtype Genome = G FilePath + +newtype GenomeH = GH Handle + +headerSize :: Int +headerSize = 2048 + +magic :: String +magic = "" + +pack :: FilePath -> IO Genome +pack fl = withFile fl ReadMode f + where f h = do l <- hGetLine h + if l == magic + then return $ G fl + else error "Bio.Seq.Query.pack: Incorrect format" + +gHOpen :: Genome -> IO GenomeH +gHOpen (G fl) = do h <- openFile fl ReadMode + return $ GH h + +gHClose :: GenomeH -> IO () +gHClose (GH h) = hClose h + +type IndexTable = M.HashMap B.ByteString (Int, Int) + +type Query = (B.ByteString, Int, Int) -- (chr, start, end), zero-based index, half-close-half-open + +getSeqs :: BioSeq s a => Genome -> [Query] -> IO [Either String (s a)] +getSeqs g querys = do gH <- gHOpen g + index <- readIndex gH + r <- mapM (getSeq gH index) querys + gHClose gH + return r +{-# INLINE getSeqs #-} + +getSeq :: BioSeq s a => GenomeH -> IndexTable -> Query -> IO (Either String (s a)) +getSeq (GH h) index (chr, start, end) = case M.lookup chr index of + Just (chrStart, chrSize) -> if end > chrSize + then return $ Left $ "Bio.Seq.getSeq: out of index: " ++ + show end ++ ">" ++ show chrSize + else do + hSeek h AbsoluteSeek $ fromIntegral $ headerSize + chrStart + start + (Right . fromBS) <$> B.hGet h (end - start) + _ -> return $ Left $ "Bio.Seq.getSeq: Cannot find " ++ show chr +{-# INLINE getSeq #-} + +getChrom :: Genome -> B.ByteString -> IO (Maybe (DNA IUPAC)) +getChrom g chr = do + chrSize <- getChrSizes g + case lookup chr chrSize of + Just s -> do [Right dna] <- getSeqs g [(chr, 0, s)] + return $ Just dna + _ -> return Nothing +{-# INLINE getChrom #-} + +getChrSizes :: Genome -> IO [(B.ByteString, Int)] +getChrSizes g = do gh <- gHOpen g + table <- readIndex gh + gHClose gh + return . map (\(k, (_, l)) -> (k, l)) . M.toList $ table +{-# INLINE getChrSizes #-} + +readIndex :: GenomeH -> IO IndexTable +readIndex (GH h) = do header <- B.hGetLine h >> B.hGetLine h + return $ M.fromList . map f . chunksOf 3 . B.words $ header + where + f [k, v, l] = (k, (readInt v, readInt l)) + f _ = error "error" +{-# INLINE readIndex #-} + +-- | indexing a genome. +mkIndex :: [FilePath] -- ^ fasta files representing individual chromosomes + -> FilePath -- ^ output file + -> IO () +mkIndex fls outFl = do + outH <- openFile outFl WriteMode + hPutStr outH $ magic ++ "\n" ++ replicate 2024 '#' -- header size: 1024 + chrs <- mapM (write outH) fls + hSeek outH AbsoluteSeek 24 + B.hPutStrLn outH $ mkHeader chrs + hClose outH + where + write handle fl = do h <- openFile fl ReadMode + fastaHeader <- B.hGetLine h + n <- loop 0 h + hClose h + return (B.tail fastaHeader, n) + where + loop !n h' = do eof <- hIsEOF h' + if eof then return n + else do l <- B.hGetLine h' + B.hPutStr handle l + loop (n + B.length l) h' +{-# INLINE mkIndex #-} + +mkHeader :: [(B.ByteString, Int)] -> B.ByteString +mkHeader xs = B.unwords.fst $ foldl f ([], 0) xs + where + f (s, i) (s', i') = (s ++ [s', B.pack $ show i, B.pack $ show i'], i + i') +{-# INLINE mkHeader #-} diff --git a/src/Bio/Utils/Functions.hs b/src/Bio/Utils/Functions.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Utils/Functions.hs @@ -0,0 +1,167 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE FlexibleContexts #-} +-------------------------------------------------------------------------------- +-- | +-- Module : $Header$ +-- Copyright : (c) 2014 Kai Zhang +-- License : MIT + +-- Maintainer : kai@kzhang.org +-- Stability : experimental +-- Portability : portable + +-- some useful functions +-------------------------------------------------------------------------------- + +module Bio.Utils.Functions ( + ihs + , ihs' + , scale + , filterFDR + , hyperquick + , kld + , jsd + , binarySearch + , binarySearchBy + , binarySearchByBounds +) where + +import Data.Bits (shiftR) +import Data.List (foldl') +import Data.Ord (comparing) +import qualified Data.Vector as V +import qualified Data.Vector.Generic as G +import qualified Data.Vector.Unboxed as U +import Statistics.Sample (meanVarianceUnb) +import Statistics.Function (sortBy) + +-- | inverse hyperbolic sine transformation +ihs :: Double -- ^ θ, determine the shape of the function + -> Double -- ^ input + -> Double +ihs !θ !x | θ == 0 = x + | otherwise = log (θ * x + sqrt (θ * θ * x * x + 1)) / θ +{-# INLINE ihs #-} + +-- | inverse hyperbolic sine transformation with θ = 1 +ihs' :: Double -> Double +ihs' = ihs 1 + +-- | scale data to zero mean and 1 variance +scale :: G.Vector v Double => v Double -> v Double +scale xs = G.map (\x -> (x - m) / sqrt s) xs + where + (m,s) = meanVarianceUnb xs +{-# INLINE scale #-} + +-- | given the p-values, filter data by controling FDR +filterFDR :: G.Vector v (a, Double) => Double -> v (a, Double) -> v (a, Double) +filterFDR α xs = go n . sortBy (comparing snd) $ xs + where + go rank v | rank <= 0 = G.empty + | snd (v `G.unsafeIndex` (rank-1)) <= fromIntegral rank * α / fromIntegral n = G.slice 0 rank v + | otherwise = go (rank-1) v + n = G.length xs +{-# INLINE filterFDR #-} + +hyperquick :: Int -> Int -> Int -> Int -> Double +hyperquick x m _n _N = loop (m-2) s s (2*e) + where + loop !k !ak !bk !epsk + | k < _N - (_n-x) - 1 && epsk > e = + let ck = ak / bk + k' = k + 1 + jjm = invJm _n x _N k' + bk' = bk * jjm + 1 + ak' = ak * jjm + espk' = fromIntegral (_N - (_n - x) - 1 - k') * (ck - ak' / bk') + in loop k' ak' bk' espk' + | otherwise = 1 - (ak / bk - epsk / 2) + s = foldl' (\s' k -> 1 + s' * invJm _n x _N k) 1.0 [x..m-2] + invJm _n x _N m = ( 1 - fromIntegral x / fromIntegral (m+1) ) / + ( 1 - fromIntegral (_n-1-x) / fromIntegral (_N-1-m) ) + e = 1e-20 + +{- +hyperquick' :: Int -> Int -> Int -> Int -> Double +hyperquick' x m _n _N = loop (m-2) s s (2*e) + where + loop !k !ak !bk !epsk + | k < _N - (_n-x) - 1 && epsk > e = + let ck = ak / bk + k' = k + 1 + jjm = invJm _n x _N k' + bk' = bk * jjm + 1 + ak' = ak * jjm + espk' = fromIntegral (_N - (_n - x) - 1 - k') * (ck - ak' / bk') + in loop k' ak' bk' espk' + | otherwise = 1 - (ak / bk - epsk / 2) + s = foldl' (\s' k -> 1 + s' * invJm _n x _N k) 1.0 [x..m-2] + invJm _n x _N m = ( 1 - fromIntegral x / fromIntegral (m+1) ) / + ( 1 - fromIntegral (_n-1-x) / fromIntegral (_N-1-m) ) + e = 1e-200 + -} + +-- | compute the Kullback-Leibler divergence between two valid (not check) probability distributions. +-- kl(X,Y) = \sum_i P(x_i) log_2(P(x_i)\/P(y_i)). +kld :: (G.Vector v Double, G.Vector v (Double, Double)) => v Double -> v Double -> Double +kld xs ys | G.length xs /= G.length ys = error "Incompitable dimensions" + | otherwise = G.foldl' f 0 . G.zip xs $ ys + where + f acc (x, y) | x == 0 = acc + | otherwise = acc + x * (logBase 2 x - logBase 2 y) +{-# SPECIALIZE kld :: U.Vector Double -> U.Vector Double -> Double #-} +{-# SPECIALIZE kld :: V.Vector Double -> V.Vector Double -> Double #-} + +-- | Jensen-Shannon divergence: JS(X,Y) = 1\/2 KL(X,(X+Y)\/2) + 1\/2 KL(Y,(X+Y)\/2). +jsd :: (G.Vector v Double, G.Vector v (Double, Double)) => v Double -> v Double -> Double +jsd xs ys = 0.5 * kld xs zs + 0.5 * kld ys zs + where zs = G.zipWith (\x y -> (x + y) / 2) xs ys +{-# SPECIALIZE jsd :: U.Vector Double -> U.Vector Double -> Double #-} +{-# SPECIALIZE jsd :: V.Vector Double -> V.Vector Double -> Double #-} + +-- | O(log n). return the position of the first element that is >= query +binarySearch :: (G.Vector v e, Ord e) + => v e -> e -> Int +binarySearch vec e = binarySearchByBounds compare vec e 0 $ G.length vec - 1 +{-# INLINE binarySearch #-} + +binarySearchBy :: G.Vector v e + => (e -> a -> Ordering) -> v e -> a -> Int +binarySearchBy cmp vec e = binarySearchByBounds cmp vec e 0 $ G.length vec - 1 +{-# INLINE binarySearchBy #-} + +binarySearchByBounds :: G.Vector v e + => (e -> a -> Ordering) -> v e -> a -> Int -> Int -> Int +binarySearchByBounds cmp vec e = loop + where + loop !l !u + | u < l = l + | otherwise = case cmp (vec G.! k) e of + LT -> loop (k+1) u + EQ -> loop l (k-1) + GT -> loop l (k-1) + where k = u + l `shiftR` 1 +{-# INLINE binarySearchByBounds #-} + +{- +empiricalCDF :: G.Vector v Double => v Double -> v Double +empiricalCDF xs = runST $ do + let n = G.length xs + indices = groupBy ( (==) `on` ((xs G.!).snd) ) $ zip [1.0..] $ sortBy (compare `on` (xs G.!)) [0..n-1] + updates mv (v,ys) = mapM_ (flip (GM.unsafeWrite mv) v.snd) ys + xs' <- G.thaw xs + mapM_ (updates xs'. ((flip (/) (fromIntegral n).fst.last) &&& id)) indices + G.unsafeFreeze xs' +{-# INLINE empiricalCDF #-} + +cdf :: G.Vector v Double => v Double -> v Double +cdf xs = let f = empiricalCDF xs + n = fromIntegral $ G.length xs + δ = 1 / (4 * (n**0.25) * sqrt (pi * log n)) + in G.map (\ x -> case () of + _ | x < δ -> δ + | x > 1 - δ -> 1 - δ + | otherwise -> x) f +{-# INLINE cdf #-} +-} diff --git a/src/Bio/Utils/Misc.hs b/src/Bio/Utils/Misc.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Utils/Misc.hs @@ -0,0 +1,43 @@ +module Bio.Utils.Misc + ( readInt + , readDouble + , bins + , binBySize + , binBySizeLeft + ) where + +import Data.ByteString.Char8 (ByteString) +import Data.ByteString.Lex.Fractional (readSigned, readExponential) +import Data.ByteString.Lex.Integral (readDecimal) +import Data.Maybe (fromMaybe) + +readInt :: ByteString -> Int +readInt x = fst . fromMaybe errMsg . readSigned readDecimal $ x + where + errMsg = error $ "readInt: Fail to cast ByteString to Int:" ++ show x +{-# INLINE readInt #-} + +readDouble :: ByteString -> Double +readDouble x = fst . fromMaybe errMsg . readSigned readExponential $ x + where + errMsg = error $ "readDouble: Fail to cast ByteString to Double:" ++ show x +{-# INLINE readDouble #-} + +-- | divide a given half-close-half-open region into fixed size +-- half-close-half-open intervals, discarding leftovers +binBySize :: Int -> (Int, Int) -> [(Int, Int)] +binBySize step (start, end) = let xs = [start, start + step .. end] + in zip xs . tail $ xs +{-# INLINE binBySize #-} + +-- | Including leftovers, the last bin may not have desired size. +binBySizeLeft :: Int -> (Int, Int) -> [(Int, Int)] +binBySizeLeft step (start, end) = let xs = [start, start + step .. end-1] ++ [end] + in zip xs . tail $ xs +{-# INLINE binBySizeLeft #-} + +-- | divide a given region into k equal size sub-regions, discarding leftovers +bins :: Int -> (Int, Int) -> [(Int, Int)] +bins binNum (start, end) = let k = (end - start) `div` binNum + in take binNum . binBySize k $ (start, end) +{-# INLINE bins #-} diff --git a/src/Bio/Utils/Overlap.hs b/src/Bio/Utils/Overlap.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Utils/Overlap.hs @@ -0,0 +1,136 @@ +module Bio.Utils.Overlap + ( overlapFragment + , overlapNucl + , coverage + ) where + +import qualified Data.ByteString.Char8 as B +import qualified Data.IntervalMap.Strict as IM +import qualified Data.HashMap.Strict as M +import qualified Data.Vector.Unboxed as V +import qualified Data.Vector.Unboxed.Mutable as VM +import Data.List +import Data.Function +import Bio.Data.Bed +import Control.Monad +import Data.Conduit +import qualified Data.Conduit.List as CL +import Control.Monad.Trans.Class (lift) + +-- | convert lines of a BED file into a data structure - A hashmap of which the +-- | chromosomes, and values are interval maps. +toMap :: [(B.ByteString, (Int, Int))] -> M.HashMap B.ByteString (IM.IntervalMap Int Int) +toMap input = M.fromList.map create.groupBy ((==) `on` (fst.fst)) $ zip input [0..] + where + f ((_, x), i) = (toInterval x, i) + create xs = (fst.fst.head $ xs, IM.fromDistinctAscList.map f $ xs) +{-# INLINE toMap #-} + +{- +-- | coverages of bins +-- FIXME: Too ugly +coverage :: [BED] -- ^ genomic locus in BED format + -> [BED] -- ^ reads in BED format + -> (V.Vector Double, Int) +coverage bin tags = getResult (V.create (VM.replicate (n+1) 0 >>= go tags)) + where + getResult v = (V.zipWith normalize (V.slice 0 n v) featWidth, v V.! n) + go ts v = do + forM_ ts (\t -> do + let set = M.lookup (t^.chrom) featMap + s = t^.chromStart + e = t^.chromEnd + b = (s, e) + l = s - e + 1 + intervals = case set of + Just iMap -> IM.intersecting iMap.toInterval $ b + _ -> [] + forM_ intervals (\interval -> do + let i = snd interval + nucl = overlap b . fst $ interval + VM.write v i . (+nucl) =<< VM.read v i + ) + VM.write v n . (+l) =<< VM.read v n + ) + return v + featMap = toMap.map (\x -> (x^.chrom, (x^.chromStart, x^.chromEnd))) $ bin + featWidth = V.fromList.map (\x -> x^.chromEnd - x^.chromStart) $ bin + n = length bin + overlap (l, u) (IM.ClosedInterval l' u') + | l' >= l = if u' <= u then u'-l'+1 else u-l'+1 + | otherwise = if u' <= u then u'-l+1 else u-l+1 + overlap _ _ = 0 + normalize a b = fromIntegral a / fromIntegral b + -} + +coverage :: [BED] -- ^ genomic locus in BED format + -> Source IO BED -- ^ reads in BED format + -> IO (V.Vector Double, Int) +coverage bin tags = liftM getResult $ tags $$ sink + where + sink :: Sink BED IO (V.Vector Int) + sink = do + v <- lift $ VM.replicate (n+1) 0 + CL.mapM_ $ \t -> do + let set = M.lookup (_chrom t) featMap + s = _chromStart t + e = _chromEnd t + b = (s, e) + l = e - s + 1 + intervals = case set of + Just iMap -> IM.intersecting iMap.toInterval $ b + _ -> [] + forM_ intervals (\interval -> do + let i = snd interval + nucl = overlap b . fst $ interval + VM.write v i . (+nucl) =<< VM.read v i + ) + VM.write v n . (+l) =<< VM.read v n + lift $ V.freeze v + getResult v = (V.zipWith normalize (V.slice 0 n v) featWidth, v V.! n) + featMap = toMap.map (\x -> (_chrom x, (_chromStart x, _chromEnd x))) $ bin + featWidth = V.fromList.map (\x -> _chromEnd x - _chromStart x) $ bin + n = length bin + overlap (l, u) (IM.ClosedInterval l' u') + | l' >= l = if u' <= u then u'-l'+1 else u-l'+1 + | otherwise = if u' <= u then u'-l+1 else u-l+1 + overlap _ _ = 0 + normalize a b = fromIntegral a / fromIntegral b + +overlapFragment, overlapNucl :: + [(Int, Int)] -- ^ Ascending order list + -> [(Int, Int)] -- ^ tags in any order + -> V.Vector Int +overlapFragment xs ts = V.create (VM.replicate n 0 >>= go ts) + where + n = length xs + iMap = IM.fromAscList $ zip (map toInterval xs) [0..] + go ts' v = do + forM_ ts' (\x -> do + let indices = snd.unzip.IM.intersecting iMap.toInterval $ x + forM_ indices (\i -> VM.write v i . (+1) =<< VM.read v i) + ) + return v + +overlapNucl xs ts = V.create (VM.replicate n 0 >>= go ts) + where + n = length xs + iMap = IM.fromAscList $ zip (map toInterval xs) [0..] + go ts' v = do + forM_ ts' (\x -> do + let intervals = IM.intersecting iMap.toInterval $ x + forM_ intervals (\interval -> do + let i = snd interval + nucl = overlap x . fst $ interval + VM.write v i . (+nucl) =<< VM.read v i + ) + ) + return v + overlap (l, u) (IM.ClosedInterval l' u') + | l' >= l = if u' <= u then u'-l'+1 else u-l'+1 + | otherwise = if u' <= u then u'-l+1 else u-l+1 + overlap _ _ = 0 + +toInterval :: (a, a) -> IM.Interval a +toInterval (l, u) = IM.ClosedInterval l u +{-# INLINE toInterval #-} diff --git a/src/Bio/Utils/Types.hs b/src/Bio/Utils/Types.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Utils/Types.hs @@ -0,0 +1,37 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE FlexibleContexts #-} +module Bio.Utils.Types + ( Sorted + , ordering + , fromSorted + , toSorted + , unsafeToSorted + ) where + +import qualified Data.Foldable as F +import Data.Ord () + +data Sorted f a where + Sorted :: (F.Foldable f, Ord a) + => { ordering :: !Ordering + , fromSorted :: !(f a) + } + -> Sorted f a + +deriving instance Show (f a) => Show (Sorted f a) + +-- | if the data has been sorted, wrap it into Sorted type +toSorted :: (F.Foldable f, Ord a) => f a -> Sorted f a +toSorted xs = Sorted o xs + where + o = snd . F.foldl' g (const EQ, EQ) $ xs + g (func, ord) x + | ord == EQ = (compare x, ord') + | ord' == ord || ord' == EQ = (compare x, ord) + | otherwise = error "data is not sorted" + where + ord' = func x + +unsafeToSorted :: (F.Foldable f, Ord a) => Ordering -> f a -> Sorted f a +unsafeToSorted = Sorted diff --git a/tests/Tests/Bed.hs b/tests/Tests/Bed.hs new file mode 100644 --- /dev/null +++ b/tests/Tests/Bed.hs @@ -0,0 +1,19 @@ +{-# LANGUAGE OverloadedStrings #-} + +module Tests.Bed (tests) where + +import Bio.Data.Bed +import Test.Tasty +import Test.Tasty.HUnit +import qualified Data.Vector as V + +tests :: TestTree +tests = testGroup "Test: Bio.Data.Bed" + [ testCase "sortBed" sortBedTest + ] + +sortBedTest :: Assertion +sortBedTest = do beds <- readBed' "tests/data/peaks.bed" :: IO [BED] + let (Sorted actual) = sortBed beds + expect <- fmap V.fromList $ readBed' "tests/data/peaks.sorted.bed" + expect @=? actual diff --git a/tests/Tests/ChIPSeq.hs b/tests/Tests/ChIPSeq.hs new file mode 100644 --- /dev/null +++ b/tests/Tests/ChIPSeq.hs @@ -0,0 +1,29 @@ +module Tests.ChIPSeq (tests) where + +import Bio.Data.Bed +import Bio.ChIPSeq +import Data.Conduit +import qualified Data.Conduit.List as CL +import Test.Tasty +import Test.Tasty.HUnit +import qualified Data.Vector as V +import Text.Printf + +peaks :: IO [BED3] +peaks = readBed' "tests/data/peaks.bed" + +tags :: Source IO BED +tags = readBed "tests/data/example.bed" + +tests :: TestTree +tests = testGroup "Test: Bio.ChIPSeq" + [ testCase "rpkm" testRPKM + ] + +testRPKM :: Assertion +testRPKM = do regions <- peaks + r1 <- tags $$ rpkmBed regions + r2 <- CL.sourceList regions $= rpkmBam "tests/data/example.bam" $$ CL.consume + let r1' = map (printf "%0.6f") . V.toList $ r1 :: [String] + r2' = map (printf "%0.6f") r2 + r1' @=? r2' diff --git a/tests/Tests/Motif.hs b/tests/Tests/Motif.hs new file mode 100644 --- /dev/null +++ b/tests/Tests/Motif.hs @@ -0,0 +1,61 @@ +{-# LANGUAGE OverloadedStrings #-} + +module Tests.Motif (tests) where + +import Test.Tasty +import Test.Tasty.HUnit +import System.Random +import Data.Default.Class +import qualified Data.Conduit.List as CL +import qualified Data.ByteString.Char8 as B +import Data.Conduit + +import Bio.Seq +import Bio.Motif +import Bio.Motif.Search +import Bio.Data.Fasta + +dna :: DNA Basic +dna = fromBS $ B.pack $ map f $ take 5000 $ randomRs (0, 3) (mkStdGen 2) + where + f :: Int -> Char + f x = case x of + 0 -> 'A' + 1 -> 'C' + 2 -> 'G' + 3 -> 'T' + _ -> undefined + +motifs :: IO [Motif] +motifs = readFasta' "tests/data/motifs.fasta" + +tests :: TestTree +tests = testGroup "Test: Bio.Motif" + [ --testCase "IUPAC converting" toIUPACTest + testCase "TFBS scanning" findTFBSTest + , testCase "pValue calculation" pValueTest + ] + + +{- +toIUPACTest :: Assertion +toIUPACTest = assertEqual "toIUPAC check" expect actual + where + expect = "SAA" + actual = show . toIUPAC $ pwm + -} + +findTFBSTest :: Assertion +findTFBSTest = do + ms <- motifs + let (Motif _ pwm) = head ms + expect <- findTFBS def pwm dna (0.6 * optimalScore def pwm) $$ CL.consume + actual <- findTFBSSlow def pwm dna (0.6 * optimalScore def pwm) $$ CL.consume + assertEqual "findTFBS" expect actual + +pValueTest :: Assertion +pValueTest = do + ms <- motifs + let expect = map (pValueToScoreExact 1e-4 def . _pwm) ms + actual = map (pValueToScore 1e-4 def . _pwm) ms + assertEqual "pValueToScore" expect actual diff --git a/tests/Tests/Seq.hs b/tests/Tests/Seq.hs new file mode 100644 --- /dev/null +++ b/tests/Tests/Seq.hs @@ -0,0 +1,19 @@ +{-# LANGUAGE OverloadedStrings #-} +module Tests.Seq (tests) where + +import Bio.Seq +import qualified Data.HashMap.Strict as M +import Test.Tasty +import Test.Tasty.HUnit + +tests :: TestTree +tests = testGroup "Test: Bio.Seq" + [ testCase "nucleotideFreq" testNuclFreq + ] + +testNuclFreq :: Assertion +testNuclFreq = do + let dna = fromBS "ACTTCCCGGGD" :: DNA IUPAC + test = M.lookupDefault undefined 'C' $ nucleotideFreq dna + expected = 4 + test @=? expected diff --git a/tests/test.hs b/tests/test.hs new file mode 100644 --- /dev/null +++ b/tests/test.hs @@ -0,0 +1,13 @@ +import qualified Tests.Bed as Bed +import qualified Tests.Motif as Motif +import qualified Tests.ChIPSeq as ChIPSeq +import qualified Tests.Seq as Seq +import Test.Tasty + +main :: IO () +main = defaultMain $ testGroup "Main" + [ Bed.tests + , Seq.tests + , ChIPSeq.tests + , Motif.tests + ]