bioinformatics-toolkit (empty) → 0.1.0
raw patch · 36 files changed
+4260/−0 lines, 36 filesdep +IntervalMapdep +aesondep +aeson-prettysetup-changed
Dependencies added: IntervalMap, aeson, aeson-pretty, base, bbi, binary, bioinformatics-toolkit, bytestring, bytestring-lexing, clustering, colour, conduit, containers, criterion, data-default-class, double-conversion, hexpat, http-conduit, matrices, mtl, palette, parallel, primitive, random, samtools, shelly, split, statistics, tasty, tasty-golden, tasty-hunit, text, transformers, unordered-containers, vector, vector-algorithms, word8
Files
- LICENSE +21/−0
- README.md +2/−0
- Setup.hs +2/−0
- benchmarks/bench.hs +53/−0
- bioinformatics-toolkit.cabal +163/−0
- data/motifs.fasta +1000/−0
- exe/mkindex.hs +11/−0
- exe/viewSeq.hs +15/−0
- src/Bio/ChIPSeq.hs +251/−0
- src/Bio/ChIPSeq/FragLen.hs +52/−0
- src/Bio/Data/Bam.hs +63/−0
- src/Bio/Data/Bed.hs +418/−0
- src/Bio/Data/BigWig.hs +32/−0
- src/Bio/Data/Fasta.hs +67/−0
- src/Bio/GO.hs +53/−0
- src/Bio/GO/GREAT.hs +132/−0
- src/Bio/GO/Parser.hs +37/−0
- src/Bio/Motif.hs +373/−0
- src/Bio/Motif/Alignment.hs +116/−0
- src/Bio/Motif/Search.hs +205/−0
- src/Bio/RealWorld/BioGRID.hs +91/−0
- src/Bio/RealWorld/ENCODE.hs +153/−0
- src/Bio/RealWorld/Ensembl.hs +41/−0
- src/Bio/RealWorld/ID.hs +27/−0
- src/Bio/RealWorld/UCSC.hs +87/−0
- src/Bio/Seq.hs +143/−0
- src/Bio/Seq/IO.hs +128/−0
- src/Bio/Utils/Functions.hs +167/−0
- src/Bio/Utils/Misc.hs +43/−0
- src/Bio/Utils/Overlap.hs +136/−0
- src/Bio/Utils/Types.hs +37/−0
- tests/Tests/Bed.hs +19/−0
- tests/Tests/ChIPSeq.hs +29/−0
- tests/Tests/Motif.hs +61/−0
- tests/Tests/Seq.hs +19/−0
- tests/test.hs +13/−0
+ LICENSE view
@@ -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.
+ README.md view
@@ -0,0 +1,2 @@+the-bioinformatician-s-toolkit+==============================
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ benchmarks/bench.hs view
@@ -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+ ]
+ bioinformatics-toolkit.cabal view
@@ -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
+ data/motifs.fasta view
@@ -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
+ exe/mkindex.hs view
@@ -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
+ exe/viewSeq.hs view
@@ -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
+ src/Bio/ChIPSeq.hs view
@@ -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 #-}
+ src/Bio/ChIPSeq/FragLen.hs view
@@ -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 #-}
+ src/Bio/Data/Bam.hs view
@@ -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 #-}
+ src/Bio/Data/Bed.hs view
@@ -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
+ src/Bio/Data/BigWig.hs view
@@ -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 #-}
+ src/Bio/Data/Fasta.hs view
@@ -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 #-}
+ src/Bio/GO.hs view
@@ -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+ -}
+ src/Bio/GO/GREAT.hs view
@@ -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
+ src/Bio/GO/Parser.hs view
@@ -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
+ src/Bio/Motif.hs view
@@ -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.+-- <doi:10.1093/nar/15.4.1353 http://nar.oxfordjournals.org/content/15/4/1353>
+ src/Bio/Motif/Alignment.hs view
@@ -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
+ src/Bio/Motif/Search.hs view
@@ -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 #-}+-}
+ src/Bio/RealWorld/BioGRID.hs view
@@ -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 #-}
+ src/Bio/RealWorld/ENCODE.hs view
@@ -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 #-}
+ src/Bio/RealWorld/Ensembl.hs view
@@ -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 #-}
+ src/Bio/RealWorld/ID.hs view
@@ -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
+ src/Bio/RealWorld/UCSC.hs view
@@ -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 #-}
+ src/Bio/Seq.hs view
@@ -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 #-}
+ src/Bio/Seq/IO.hs view
@@ -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:+-- <HASKELLBIOINFORMATICS>\nCHR1 START SIZE+newtype Genome = G FilePath++newtype GenomeH = GH Handle++headerSize :: Int+headerSize = 2048++magic :: String+magic = "<HASKELLBIOINFORMATICS>"++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 #-}
+ src/Bio/Utils/Functions.hs view
@@ -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 #-}+-}
+ src/Bio/Utils/Misc.hs view
@@ -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 #-}
+ src/Bio/Utils/Overlap.hs view
@@ -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 #-}
+ src/Bio/Utils/Types.hs view
@@ -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
+ tests/Tests/Bed.hs view
@@ -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
+ tests/Tests/ChIPSeq.hs view
@@ -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'
+ tests/Tests/Motif.hs view
@@ -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
+ tests/Tests/Seq.hs view
@@ -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
+ tests/test.hs view
@@ -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+ ]