cobot-tools (empty) → 0.1.0.1
raw patch · 18 files changed
+905/−0 lines, 18 filesdep +QuickCheckdep +arraydep +basesetup-changed
Dependencies added: QuickCheck, array, base, cobot, cobot-tools, containers, data-msgpack, deepseq, hspec, lens, mtl, neat-interpolation, text
Files
- ChangeLog.md +11/−0
- LICENSE +30/−0
- README.md +38/−0
- Setup.hs +2/−0
- cobot-tools.cabal +86/−0
- src/Bio/Tools/Sequence/Primers/Constants.hs +64/−0
- src/Bio/Tools/Sequence/Primers/Optimization.hs +391/−0
- src/Bio/Tools/Sequence/Primers/Properties.hs +18/−0
- src/Bio/Tools/Sequence/Primers/Types.hs +26/−0
- src/Bio/Tools/Sequence/ViennaRNA/Cofold.hs +5/−0
- src/Bio/Tools/Sequence/ViennaRNA/Fold.hs +5/−0
- src/Bio/Tools/Sequence/ViennaRNA/Internal/Cofold.hs +37/−0
- src/Bio/Tools/Sequence/ViennaRNA/Internal/Fold.hs +31/−0
- src/Bio/Tools/Sequence/ViennaRNA/Internal/RNALike.hs +17/−0
- src/Bio/Tools/Sequence/ViennaRNA/Native/vienna_rna_wrapper.c +47/−0
- test/Spec.hs +15/−0
- test/SpecPrimers.hs +48/−0
- test/SpecViennaRNA.hs +34/−0
+ ChangeLog.md view
@@ -0,0 +1,11 @@+# Changelog for cobot-tools++## [Unreleased]++## [0.1.0.1] - 2019-07-10+### Added+- `extra-lib-dirs` parameter to `stack.yaml`.++## [0.1.0.0] - 2019-07-01+### Added+- Tool for primer design.
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Author name here (c) 2019++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Author name here nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,38 @@+# cobot-tools+Tools for computational biology++# Software required to use cobot-tools++## ViennaRNA++ViennaRNA is a library that already contains lots of useful algorithms that work+with RNA sequences. We use it in order to make tools from cobot-tools better.++Installation guide for Linux can be found [here](https://www.tbi.univie.ac.at/RNA/).++Some issues occur when you try to install ViennaRNA for MacOS. If you want to avoid+them, use the following guide.++### Installation of ViennaRNA for MacOS++Installation of ViennaRNA for MacOS can be done in following steps:++ 1. Download library. Here is the [link](https://www.tbi.univie.ac.at/RNA/download/sourcecode/2_4_x/ViennaRNA-2.4.13.tar.gz) + for version 2.4.13.+ 2. Unpack it.+ 3. Go to the directory, where archive has been unpacked.+ 4. Run command `./configure --without-perl --without-python`.+ 5. Remove field `libRNA_la_LDFLAGS` and flag `-static` from the file `src/ViennaRNA/Makefile.am`.+ 6. Due to us having changed .am-file, it is needed to install following packages: autoconf and automake. In order to do it, run command: `brew install autoconf automake`.+ 7. Run in the root of library command `autoreconf`.+ 8. Compile library using command `make -j8` (`j8` - flag that allows us to compile code using 8 jobs).+ 9. Run command `make install`.+ +# Tools of cobot-tools++## Sequence.Primer.Optimisation.designPrimer++Given 'DNA' sequence and position in that sequence designs forward primer for that sequence. +Primer will start at the given position.++
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ cobot-tools.cabal view
@@ -0,0 +1,86 @@+-- This file has been generated from package.yaml by hpack version 0.28.2.+--+-- see: https://github.com/sol/hpack+--+-- hash: 62f7da4d659cbc72389c4ea8d12693721e66abd1462ef2f2493c83c1601bf809++name: cobot-tools+version: 0.1.0.1+synopsis: Biological data file formats and IO+description: Please see the README on GitHub at <https://github.com/less-wrong/cobot-tools#readme>+category: Bio+homepage: https://github.com/less-wrong/cobot-tools#readme+bug-reports: https://github.com/less-wrong/cobot-tools/issues+author: Pavel Yakovlev, Bogdan Neterebskii, Alexander Sadovnikov+maintainer: pavel@yakovlev.me+copyright: 2018-2019, Less Wrong Bio+license: BSD3+license-file: LICENSE+build-type: Simple+cabal-version: >= 1.10+extra-source-files:+ ChangeLog.md+ README.md++source-repository head+ type: git+ location: https://github.com/less-wrong/cobot-tools++library+ exposed-modules:+ Bio.Tools.Sequence.Primers.Constants+ Bio.Tools.Sequence.Primers.Optimization+ Bio.Tools.Sequence.Primers.Properties+ Bio.Tools.Sequence.Primers.Types+ Bio.Tools.Sequence.ViennaRNA.Cofold+ Bio.Tools.Sequence.ViennaRNA.Fold+ Bio.Tools.Sequence.ViennaRNA.Internal.Cofold+ Bio.Tools.Sequence.ViennaRNA.Internal.Fold+ Bio.Tools.Sequence.ViennaRNA.Internal.RNALike+ other-modules:+ Paths_cobot_tools+ hs-source-dirs:+ src+ default-extensions: DeriveGeneric DeriveFunctor DeriveFoldable DeriveAnyClass FlexibleInstances InstanceSigs MultiParamTypeClasses RecordWildCards ScopedTypeVariables OverloadedStrings TypeFamilies DataKinds ConstraintKinds TypeOperators TemplateHaskell FlexibleContexts+ c-sources:+ src/Bio/Tools/Sequence/ViennaRNA/Native/vienna_rna_wrapper.c+ extra-libraries:+ RNA+ build-depends:+ array >=0.5 && <0.6+ , base >=4.7 && <5+ , cobot+ , containers >=0.5.7.1 && <0.7+ , data-msgpack >=0.0.9 && <0.1+ , deepseq >=1.4 && <1.5+ , lens >=4.16 && <5.0+ , mtl >=2.2.1 && <2.3.0+ , text+ default-language: Haskell2010++test-suite cobot-tools-test+ type: exitcode-stdio-1.0+ main-is: Spec.hs+ other-modules:+ SpecPrimers+ SpecViennaRNA+ Paths_cobot_tools+ hs-source-dirs:+ test+ default-extensions: OverloadedStrings TypeFamilies+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ build-depends:+ QuickCheck >=2.9.2 && <2.13+ , array >=0.5 && <0.6+ , base >=4.7 && <5+ , cobot+ , cobot-tools+ , containers >=0.5.7.1 && <0.7+ , data-msgpack >=0.0.9 && <0.1+ , deepseq >=1.4 && <1.5+ , hspec >=2.4.1 && <2.7+ , lens >=4.16 && <5.0+ , mtl >=2.2.1 && <2.3.0+ , neat-interpolation >=0.3+ , text+ default-language: Haskell2010
+ src/Bio/Tools/Sequence/Primers/Constants.hs view
@@ -0,0 +1,64 @@+module Bio.Tools.Sequence.Primers.Constants+ ( annealingTemp+ , bindingRate+ , eTgtEFoldRel+ , maxPrimerGCContent+ , maxPrimerLength+ , maxPrimerMeltingTemp+ , minPrimerGCContent+ , minPrimerLength+ , minPrimerMeltingTemp+ , topN+ ) where++-- | Minimum recommended primer's length.+--+minPrimerLength :: Int+minPrimerLength = 17++-- | Maximum recommended primer's length.+--+maxPrimerLength :: Int+maxPrimerLength = 35++-- | Temperature under which annealing of primers happens.+--+annealingTemp :: Double+annealingTemp = 62++-- | Minimum recommended primer's melting temperature.+--+minPrimerMeltingTemp :: Int+minPrimerMeltingTemp = 55++-- | Maximum recommended primer's melting temperature.+--+maxPrimerMeltingTemp :: Int+maxPrimerMeltingTemp = 75++-- | Minimum recommended primer's GC-content.+--+minPrimerGCContent :: Float+minPrimerGCContent = 40++-- | Maximum recommended primer's GC-content.+--+maxPrimerGCContent :: Float+maxPrimerGCContent = 60++-- | Rate of nucleotides from primer that should bind to target area of source sequence.+--+bindingRate :: Float+bindingRate = 0.75++-- | Minimum allowed value of primer's tagret interaction energy to+-- primer's folding energy relation.+--+eTgtEFoldRel :: Float+eTgtEFoldRel = 2.7++-- | Number of candidates that are taken if none of candidates satisfies+-- GC-content condition.+--+topN :: Int+topN = 10
+ src/Bio/Tools/Sequence/Primers/Optimization.hs view
@@ -0,0 +1,391 @@+{-# LANGUAGE ViewPatterns #-}++module Bio.Tools.Sequence.Primers.Optimization+ ( designPrimer+ ) where++import Bio.Chain (fromList)+import Bio.Chain.Alignment (AffineGap (..),+ AffineGap2,+ AlignmentResult (..),+ LocalAlignment (..),+ Operation (..), align)+import Bio.NucleicAcid.Chain (NucleicAcidChain (..))+import Bio.NucleicAcid.Nucleotide (Complementary (..),+ DNA (..), symbol)+import Bio.Tools.Sequence.Primers.Constants as Constants (annealingTemp,+ bindingRate,+ eTgtEFoldRel,+ maxPrimerGCContent,+ maxPrimerLength,+ maxPrimerMeltingTemp,+ minPrimerGCContent,+ minPrimerLength,+ minPrimerMeltingTemp,+ topN)+import Bio.Tools.Sequence.Primers.Types (Primer,+ ScoredPrimer (..), eTgt,+ eTgtEFold, gcContent,+ meltingTemp, seq')+import Bio.Tools.Sequence.ViennaRNA.Cofold (cofold)+import Bio.Tools.Sequence.ViennaRNA.Fold (fold)+import Control.Lens ((&), (.~), (^.))+import Control.Monad (when)+import Control.Monad.Except (MonadError, throwError)+import Data.List (sortOn)+import Data.Maybe (catMaybes)+import Data.Text (Text)++-- | Given 'DNA' sequence and position in that sequence designs forward primer+-- for that sequence. Primer will start at the given position. @isCyclic@ marks+-- whether the sequence that we design primer for is cyclic or not.+--+-- Flag @isCyclic@ also defines type of algorithm that will be used to design primer.+--+-- If @isCyclic@ is set to False, ViennaRNA will be used to check that primer+-- has no off-target interactions with given sequence. This is pretty accurate method,+-- but when it is used to process long sequences (more than 1000 bps), it's quite slow.+--+-- If @isCyclic@ is set to True, algorithm that uses several heuristics will be+-- used to check that primer has no off-target interactions with given sequence.+-- This method is less accurate than using ViennaRNA to calculate energy of primer's interaction+-- with sequence, but much more faster.+--+-- We did such segregation of algorithms, because cyclic DNA sequences (plasmids)+-- are very long and non-cyclic sequences that we work with are no longer than 1000 bps.+--+designPrimer :: MonadError Text m => Bool -> [DNA] -> Int -> m ScoredPrimer+designPrimer isCyclic dna' pos = do+ when (pos >= length dna') $ throwError "Bio.Tools.Sequence.Primers.Optimization: given position is out of range."++ withE <- energyFilter . gcContentFilter . lengthFilter $ candidates++ case gcOnEndFilter withE of+ [] -> throwError badPrimersError+ l -> pure $ head $ sortOn (fmap negate . (^. gcContent)) l+ where+ dna = if isCyclic then dna' <> take Constants.maxPrimerLength dna' else dna'+ candidates = genTemperatureCandidates (drop pos dna) []++ badPrimersError :: Text+ badPrimersError = "Bio.Tools.Sequence.Primers.Optimization: all primers designed from given position are seriously flawed."++ -- | Generates candidate primers with melting temperature in needed range.+ --+ genTemperatureCandidates :: [DNA] -> [DNA] -> [ScoredPrimer]+ genTemperatureCandidates [] _ = []+ genTemperatureCandidates (x : xs) base = res+ where+ newBase = base <> [x]+ newPrimer = toScored newBase+ primerWithTemp = calcMeltingTemp newPrimer++ Just mTemp = primerWithTemp ^. meltingTemp++ res | mTemp < Constants.minPrimerMeltingTemp = genTemperatureCandidates xs newBase+ | mTemp > Constants.maxPrimerMeltingTemp = []+ | otherwise = primerWithTemp : genTemperatureCandidates xs newBase++ toScored :: Primer -> ScoredPrimer+ toScored p = ScoredPrimer p Nothing Nothing Nothing Nothing++ -- | Filters 'ScoredPrimer's based on their length.+ --+ lengthFilter :: [ScoredPrimer] -> [ScoredPrimer]+ lengthFilter = filter lengthPred+ where+ lengthPred :: ScoredPrimer -> Bool+ lengthPred sp = Constants.minPrimerLength <= ls && ls <= Constants.maxPrimerLength+ where+ ls = length $ sp ^. seq'++ -- | Leaves only 'ScoredPrimer's that have GC on their 3' end.+ -- If no such primers are found, filtering doesn't happen.+ --+ gcOnEndFilter :: [ScoredPrimer] -> [ScoredPrimer]+ gcOnEndFilter sps | null filtered = sps+ | otherwise = filtered+ where+ filtered = filter gcOnEndPred sps++ gcOnEndPred :: ScoredPrimer -> Bool+ gcOnEndPred = (`elem` [DC, DG]) . last . (^. seq')++ -- | Leaves primers whose GC-content is in needed range. If no such primers+ -- are found, leaves @Constants.topN@ primers with GC-content nearest to needed range.+ --+ gcContentFilter :: [ScoredPrimer] -> [ScoredPrimer]+ gcContentFilter sps | null filtered = take Constants.topN sorted+ | otherwise = filtered+ where+ primersWithGC = fmap calcGCContent sps+ filtered = filter gcContentPred primersWithGC++ sorted = sortOn scoringFunc primersWithGC++ gcContentPred :: ScoredPrimer -> Bool+ gcContentPred sp = Constants.minPrimerGCContent <= gcc && gcc <= Constants.maxPrimerGCContent+ where+ Just gcc = sp ^. gcContent++ scoringFunc :: ScoredPrimer -> Float+ scoringFunc sp = min (abs $ Constants.minPrimerGCContent - gcc) (abs $ Constants.maxPrimerGCContent - gcc)+ where+ Just gcc = sp ^. gcContent++ -- | Filters primers using energy characteristics.+ --+ energyFilter :: MonadError Text m => [ScoredPrimer] -> m [ScoredPrimer]+ energyFilter sps = fmap eTgtEFoldFilter . eTgtFilter $ sps++ -- | Leaves only primers whose relation of energy of interaction with target+ -- to energy of forming a secondary structure is higher then @Constants.eTgtEFoldRel@.+ --+ eTgtEFoldFilter :: [ScoredPrimer] -> [ScoredPrimer]+ eTgtEFoldFilter s = filter eTgtEFoldPred . fmap calcTargetFold $ s+ where+ eTgtEFoldPred :: ScoredPrimer -> Bool+ eTgtEFoldPred sp = etef >= Constants.eTgtEFoldRel+ where+ Just etef = sp ^. eTgtEFold++ -- | Leaves only primers that bind to target on source sequence.+ --+ eTgtFilter :: MonadError Text m => [ScoredPrimer] -> m [ScoredPrimer]+ eTgtFilter s | null res = throwError badPositionError+ | otherwise = pure res+ where+ res = catMaybes $ fmap (calcTarget isCyclic dna pos) s++ badPositionError :: Text+ badPositionError = "Bio.Tools.Sequence.Primers.Optimization: all primers designed from given position don't bind to target."++-- | Calculates melting temperature for 'ScoredPrimer'.+-- Temperature is calculated using the following formula: 4 * (C + G) + 2 * (A + T).+--+calcMeltingTemp :: ScoredPrimer -> ScoredPrimer+calcMeltingTemp sp = sp & meltingTemp .~ Just temperature+ where+ temperature = sum $ fmap tempForNuc (sp ^. seq')++ tempForNuc :: DNA -> Int+ tempForNuc DC = 4+ tempForNuc DG = 4+ tempForNuc DA = 2+ tempForNuc DT = 2++-- | Calculates GC-content for 'ScoredPrimer'.+-- GC-content is calulcated using the following formula: (G + C) / (A + T + G + C).+--+calcGCContent :: ScoredPrimer -> ScoredPrimer+calcGCContent sp = sp & gcContent .~ Just content+ where+ nGC = length $ filter (\x -> x == DC || x == DG) $ sp ^. seq'+ content = fromIntegral nGC / fromIntegral (length $ sp ^. seq')++-- | Calculates relation of 'ScoredPrimers's target interaction energy to+-- its folding energy.+--+calcTargetFold :: ScoredPrimer -> ScoredPrimer+calcTargetFold sp = sp & eTgtEFold .~ Just tgtToFold+ where+ foldingEnergy = fst $ fold Constants.annealingTemp $ sp ^. seq'++ -- @calcTargetFold@ is used only after 'ScoredPrimer's target interaction energy+ -- has been calculated+ Just tgt = sp ^. eTgt+ tgtToFold = abs $ tgt / toEps foldingEnergy++-- | Calculates energy of interaction of given 'DNA' and given 'DNA' sequence's+-- complementary strand.+--+cofoldEnergy :: [DNA] -> [DNA] -> Float+cofoldEnergy sp s = abs $ fst $ cofold Constants.annealingTemp (reverse sp, fmap cNA s)++-- | Local alignment algorithm that aligns arguments based on their complementarity+-- and doesn't allow gaps on the second argument of the alignment (that is primer in our terms).+--+-- We don't allow gaps on the primer during alignment, because this situation+-- is not physical: it means that hairpins can appear on the plasmid.+--+localAlignment :: LocalAlignment AffineGap2 DNA DNA+localAlignment = LocalAlignment complementScoring (plasmidGap, primerGap)+ where+ plasmidGap :: AffineGap+ plasmidGap = AffineGap (-5) (-1)++ primerGap :: AffineGap+ primerGap = AffineGap (-1000) (-1000)++ complementScoring :: DNA -> DNA -> Int+ complementScoring DA (cNA -> DT) = 3+ complementScoring DT (cNA -> DA) = 3+ complementScoring DC (cNA -> DG) = 5+ complementScoring DG (cNA -> DC) = 5+ complementScoring _ _ = -3++-- | Calculates energy of interaction of primer @sp@ with sequence @dna@.+-- If @sp@ doesn't bind to target area on @dna@ (target area starts at @sourceInd@),+-- then Nothing is returned.+--+-- There are different algorithms to check off-target interaction depending of whether+-- @dna@ is cyclic or not. It is defined by the first parameter.+--+calcTarget :: Bool -> [DNA] -> Int -> ScoredPrimer -> Maybe ScoredPrimer+-- algorithm for cyclic sequences+calcTarget True dna sourceInd sp | null offTargets || tgtE > maxOffTarget = res+ | otherwise = Nothing+ where+ primerSeq = sp ^. seq'+ primerLen = length primerSeq++ -- position in the @dna@, where primer reaches half of its length+ primerLenHalf = sourceInd + primerLen `div` 2++ -- cut all @matchingSeqs@ into k-mers of all sizes in range [@primerLenHalf@; @primerLen@].+ offTargets = concatMap seqToKMers matchingSeqs++ tgtE = cofoldEnergy primerSeq primerSeq+ maxOffTarget = maximum $ fmap (cofoldEnergy primerSeq) offTargets++ res = pure $ sp & eTgt .~ Just tgtE++ seqToKMers :: [DNA] -> [[DNA]]+ seqToKMers s = concatMap (toKMers s) [primerLenHalf .. primerLen]+ where+ toKMers :: [a] -> Int -> [[a]]+ toKMers l k | length l < k = []+ | otherwise = take k l : toKMers (tail l) k++ -- | All parts of @dna@ that align the best on given @sp@ using algorithm @localAlignment@.+ -- We consider complementary strands of these parts to be most energetically preferrable+ -- off-targets for our primer.+ --+ matchingSeqs :: [[DNA]]+ matchingSeqs = fmap matchingSeq [invertedDna, dnaRevComp]+ where+ -- we invert the cyclic @dna@ in such way that target area is being cut in half,+ -- so that it won't be considered as off-target interaction+ invertedDna | [x] <- invertPoints = drop x dna <> drop Constants.maxPrimerLength (take x dna)+ | [x, y] <- invertPoints = take (y - x) . drop x $ dna+ | otherwise = error "This branch is never visited."++ -- here we check that @primerLenHalf@ is not in the part of sequence,+ -- that was appended at the end to create cyclic sequence.+ -- If it's not in this part, invert at @primerLenHalf@.+ -- Otherwise consider everything in between @modPos@ and (@modPos@ + (length @dna@ - @Constants.maxPrimerLength@))+ -- as inverted sequence.+ modPos = primerLenHalf `mod` (length dna - Constants.maxPrimerLength)+ invertPoints | modPos >= Constants.maxPrimerLength = [primerLenHalf]+ | otherwise = [modPos, modPos + (length dna - Constants.maxPrimerLength)]++ -- off-target interaction could happen in reversed direction+ dnaRevComp = reverse $ fmap cNA dna++ matchingSeq :: [DNA] -> [DNA]+ matchingSeq dna' = res'+ where+ ar = alignment $ alignmentFunc dna' primerSeq++ traceStart = toCoord $ last ar+ traceEnd = toCoord $ head ar++ -- (traceStart, traceEnd) is inclusive range that describes off-target area+ -- in @dna'@. We want to extend that range by @primerLen@ in both directions+ -- to consider different possibilities+ (l, r) = (max 0 (traceEnd - primerLen), min (length dna' - 1) (traceStart + primerLen))+ res' = take (r - l + 1) $ drop l dna'++ toCoord :: Operation Int Int -> Int+ toCoord = getI++ alignmentFunc :: [DNA] -> [DNA] -> AlignmentResult (NucleicAcidChain Int DNA) (NucleicAcidChain Int DNA)+ alignmentFunc plasmid primer = alRes+ where+ plasmidC = toChain plasmid+ primerC = toChain primer++ alRes = align localAlignment plasmidC primerC++ toChain :: [DNA] -> NucleicAcidChain Int DNA+ toChain = NucleicAcidChain . fromList+-- algorithm for linear sequences+calcTarget False dna sourceInd sp | Just e <- primerTargetEnergy = pure $ sp & eTgt .~ Just e+ | otherwise = Nothing+ where+ primerTargetEnergy :: Maybe Float+ primerTargetEnergy | abs e < revE = Nothing+ | Just cnt <- actualBindingRateM, checkBindingRate cnt = Just e+ | otherwise = Nothing+ where+ primerSeq = sp ^. seq'+ spStr = symbol <$> primerSeq++ compPrimerSeq = fmap cNA dna++ (e, bindStr) = cofold Constants.annealingTemp (reverse primerSeq, compPrimerSeq)+ primerLength = length spStr++ -- energy of primer's interaction revesed dna strand+ revE = cofoldEnergy primerSeq (reverse $ fmap cNA dna)++ -- inclusive range in which target binding site is being contained in @bindStr@+ (lInd, rInd) = (primerLength + sourceInd, lInd + primerLength - 1)++ bindWithInds = zip bindStr [0..]+ actualBindingRateM = calcMatchesInBindingSite [] 0 bindWithInds++ -- | Calculates number of nucleotides in target binding site that bind+ -- to primer. Also checks that last nucleotide of primer binds to last nucleotide+ -- of target binding site. If this condition is not satisfied, Nothing is returned.+ --+ -- @bindStr@ is used for these calculations. @bindStr@ represents interaction between primer+ -- and sequence in form of a dot plot. Dot plot also contains balanced bracket sequence+ -- that shows how nucleotides bind to each other.+ --+ calcMatchesInBindingSite :: [Int] -> Int -> [(Char, Int)] -> Maybe Int+ calcMatchesInBindingSite _ cnt [] = Just cnt+ calcMatchesInBindingSite stack cnt (x : xs) -- if we encounter close bracket, then we check that its position is in target binding site+ -- and position of open bracket that corresponds to it is in primer+ | c == closeBracket, lInd <= i && i < rInd && i' < primerLength = calcMatchesInBindingSite l (cnt + 1) xs+ -- close bracket in target corresponds to non-primer interaction, we skip this bracket+ | c == closeBracket, lInd <= i && i < rInd = calcMatchesInBindingSite l cnt xs+ -- next two conditions check that 3' end of primer binds to end of the target binding site+ | c == closeBracket, i == rInd, i' == 0 = Just $ cnt + 1+ | c == closeBracket, i == rInd, i' /= 0 = Nothing+ -- close bracket is out of target range, we skip it+ | c == closeBracket = calcMatchesInBindingSite l cnt xs+ -- open brackets are put on stack+ | c == openBracket = calcMatchesInBindingSite (i : stack) cnt xs+ -- we ignore dots+ | otherwise = calcMatchesInBindingSite stack cnt xs+ where+ (c, i) = x++ -- since the bracket sequence is balanced, we will get to this pattern-matching only if+ -- there is something on top of the stack+ (i' : l) = stack++ openBracket :: Char+ openBracket = '('++ closeBracket :: Char+ closeBracket = ')'++ -- | Checks that not less then @Constants.bindingRate@ * 100 precents of primer's nucleotides+ -- interact with target.+ --+ checkBindingRate :: Int -> Bool+ checkBindingRate cnt = fromIntegral cnt / fromIntegral primerLength >= Constants.bindingRate++--------------------------------------------------------------------------------+-- Utility functions.+--------------------------------------------------------------------------------++-- | This function is used to avoid dividing by zero.+--+toEps :: Float -> Float+toEps x | abs x < eps = eps+ | otherwise = x+ where+ eps = 0.01
+ src/Bio/Tools/Sequence/Primers/Properties.hs view
@@ -0,0 +1,18 @@+module Bio.Tools.Sequence.Primers.Properties+ ( calculateTargetToFold+ ) where++import Bio.NucleicAcid.Nucleotide (Complementary (..))+import Bio.Tools.Sequence.Primers.Constants (annealingTemp)+import Bio.Tools.Sequence.Primers.Types (Primer)+import Bio.Tools.Sequence.ViennaRNA.Cofold (cofold)+import Bio.Tools.Sequence.ViennaRNA.Fold (fold)++-- | For given primer calculates relation of this primer's energy of interaction+-- with target to its energy of forming secondary structure.+--+calculateTargetToFold :: Primer -> Float+calculateTargetToFold primer = tgtEnergy / foldingEnergy+ where+ tgtEnergy = abs $ fst $ cofold annealingTemp (reverse primer, fmap cNA primer)+ foldingEnergy = abs $ fst $ fold annealingTemp primer
+ src/Bio/Tools/Sequence/Primers/Types.hs view
@@ -0,0 +1,26 @@+{-# LANGUAGE TemplateHaskell #-}++module Bio.Tools.Sequence.Primers.Types where++import Bio.NucleicAcid.Nucleotide (DNA)+import Control.Lens (makeLenses)+++-- | Primer is just an alias for sequence of nucleotides.+--+type Primer = [DNA]++-- | Primer with its characteristics.+--+data ScoredPrimer = ScoredPrimer { _seq' :: Primer -- ^ primer's sequence+ , _meltingTemp :: Maybe Int -- ^ melting temperature+ , _gcContent :: Maybe Float -- ^ GC-content+ , _eTgtEFold :: Maybe Float -- ^ relation of energy of interaction with+ -- target to energy of primer's forming+ -- secondary structure+ , _eTgt :: Maybe Float -- ^ energy of interaction with target if+ -- primer binds to target. Otherwise set to 0+ }+ deriving (Eq, Show)++makeLenses ''ScoredPrimer
+ src/Bio/Tools/Sequence/ViennaRNA/Cofold.hs view
@@ -0,0 +1,5 @@+module Bio.Tools.Sequence.ViennaRNA.Cofold+ ( cofold+ ) where++import Bio.Tools.Sequence.ViennaRNA.Internal.Cofold (cofold)
+ src/Bio/Tools/Sequence/ViennaRNA/Fold.hs view
@@ -0,0 +1,5 @@+module Bio.Tools.Sequence.ViennaRNA.Fold+ ( fold+ ) where++import Bio.Tools.Sequence.ViennaRNA.Internal.Fold (fold)
+ src/Bio/Tools/Sequence/ViennaRNA/Internal/Cofold.hs view
@@ -0,0 +1,37 @@+module Bio.Tools.Sequence.ViennaRNA.Internal.Cofold+ ( cofold+ ) where++import Bio.NucleicAcid.Nucleotide (symbol)+import Bio.Tools.Sequence.ViennaRNA.Internal.RNALike (RNALike (..))+import Foreign.C.String (CString,+ newCString,+ peekCString)+import Foreign.C.Types (CDouble (..),+ CFloat (..))+import Foreign.Marshal.Alloc (free)+import System.IO.Unsafe (unsafePerformIO)++foreign import ccall "vrna_cofold_temperature" vrna_cofold_temperature :: CString -> CString -> CDouble -> CFloat++-- TODO make it for list of pairs to allocate CString of only one size+-- TODO (need to handle different sizes of oligs, try to take just max size (because CString ends with \0)+vRnaCofoldString :: Double -> (String, String) -> (Float, String)+vRnaCofoldString temperature (rnaSequence1, rnaSequence2) = unsafePerformIO $ do+ let inputSeq = rnaSequence1 ++ "&" ++ rnaSequence2+ cRnaSequence <- newCString' inputSeq+ let energy = realToFrac $ vrna_cofold_temperature cRnaSequence cRnaSequence (realToFrac temperature)+ structureRes <- evalResult energy cRnaSequence+ free' cRnaSequence+ return (energy, structureRes)+ where+ newCString' = newCString+ free' = free+ evalResult energy cRnaSequence = energy `seq` peekCString cRnaSequence++-- | Calculates cofolding energy and interaction dot-plot between two 'RNALike' strands.+--+cofold :: RNALike a => Double -> ([a], [a]) -> (Float, String)+cofold temperature = vRnaCofoldString temperature . toStringsPair+ where+ toStringsPair (nucs1, nucs2) = (symbol . toRNA <$> nucs1, symbol . toRNA <$> nucs2)
+ src/Bio/Tools/Sequence/ViennaRNA/Internal/Fold.hs view
@@ -0,0 +1,31 @@+module Bio.Tools.Sequence.ViennaRNA.Internal.Fold+ ( fold+ ) where++import Bio.NucleicAcid.Nucleotide (symbol)+import Bio.Tools.Sequence.ViennaRNA.Internal.RNALike (RNALike (..))+import Foreign.C.String (CString,+ newCString,+ peekCString)+import Foreign.C.Types (CDouble (..),+ CFloat (..))+import Foreign.Marshal.Alloc (free)+import System.IO.Unsafe (unsafePerformIO)++foreign import ccall "vrna_fold_temperature" vrna_fold_temperature :: CString -> CString -> CDouble -> CFloat++vRnaFoldString :: Double -> String -> (Float, String)+vRnaFoldString temperature rnaSequence = unsafePerformIO $ do+ cRnaSequence <- newCString rnaSequence+ cStructure <- newCString rnaSequence -- use the same rnaSequence just to get string of the same size+ let energy = realToFrac $ vrna_fold_temperature cRnaSequence cStructure (realToFrac temperature)+ structureRes <- energy `seq` peekCString cStructure+ free cStructure+ free cRnaSequence+ return (energy, structureRes)++-- | Calculates folding energy of given 'RNALike' strand. Also returns secondary+-- structure of that strand in dot-plot form.+--+fold :: RNALike a => Double -> [a] -> (Float, String)+fold temperature = vRnaFoldString temperature . (symbol . toRNA <$>)
+ src/Bio/Tools/Sequence/ViennaRNA/Internal/RNALike.hs view
@@ -0,0 +1,17 @@+module Bio.Tools.Sequence.ViennaRNA.Internal.RNALike+ ( RNALike (..)+ ) where++import Bio.NucleicAcid.Nucleotide (DNA, RNA)+import qualified Bio.NucleicAcid.Nucleotide as N (toRNA)++-- | Class that describes objects that can be converted to RNA.+--+class RNALike a where+ toRNA :: a -> RNA++instance RNALike DNA where+ toRNA = N.toRNA++instance RNALike RNA where+ toRNA = id
+ src/Bio/Tools/Sequence/ViennaRNA/Native/vienna_rna_wrapper.c view
@@ -0,0 +1,47 @@+#include <ViennaRNA/model.h>+#include <ViennaRNA/fold_compound.h>+#include <ViennaRNA/mfe.h>++float+vrna_fold_temperature(+ const char *string,+ char *structure,+ double temperature) {++ float mfe;+ vrna_fold_compound_t *vc;+ vrna_md_t md;++ vrna_md_set_default(&md);+ md.temperature = temperature;+ vc = vrna_fold_compound(string, &md, 0);+ mfe = vrna_mfe(vc, structure);++ vrna_fold_compound_free(vc);++ return mfe;+}++float+vrna_cofold_temperature(+ const char *seq,+ char *structure,+ double temperature) {++ float mfe;+ vrna_fold_compound_t *vc;+ vrna_md_t md;++ vrna_md_set_default(&md);+ md.temperature = temperature;+ md.min_loop_size = 0; /* set min loop length to 0 */++ /* get compound structure */+ vc = vrna_fold_compound(seq, &md, 0);++ mfe = vrna_mfe_dimer(vc, structure);++ vrna_fold_compound_free(vc);++ return mfe;+}
+ test/Spec.hs view
@@ -0,0 +1,15 @@+import SpecPrimers+import SpecViennaRNA+import System.IO+import Test.Hspec++main :: IO ()+main = do+ hSetBuffering stdout NoBuffering+ hspec $ do+ -- Primers+ testPrimers++ -- ViennaRNA+ foldTest+ cofoldTest
+ test/SpecPrimers.hs view
@@ -0,0 +1,48 @@+{-# LANGUAGE OverloadedStrings #-}++module SpecPrimers where++import Bio.NucleicAcid.Nucleotide (DNA, symbol)+import Bio.Tools.Sequence.Primers.Optimization+import Bio.Tools.Sequence.Primers.Types+import Test.Hspec++testPrimers :: SpecWith ()+testPrimers = describe "Primer optimization test" testPrimersOptimization++testPrimersOptimization :: SpecWith ()+testPrimersOptimization = do+ it "Sequence: GATCCAACTTCAAAGAGTCCTGGCCAGACATGATAAGATACATTGATGAGTTTGGACAAACCACAACTAGAATGCAGTGAAAAAAATGCTTTATTTGTGAAATTTGTGATGC; pos: 0" $+ let s = "GATCCAACTTCAAAGAGTCCTGGCCAGACATGATAAGATACATTGATGAGTTTGGACAAACCACAACTAGAATGCAGTGAAAAAAATGCTTTATTTGTGAAATTTGTGATGC" :: [DNA]+ in fmap symbol . _seq' <$> designPrimer False s 0 `shouldBe` Right "GATCCAACTTCAAAGAGTCCTGGC"+ it "Sequence: GATCCAACTTCAAAGAGTCCTGGCCAGACATGATAAGATACATTGATGAGTTTGGACAAACCACAACTAGAATGCAGTGAAAAAAATGCTTTATTTGTGAAATTTGTGATGC; pos: 61" $+ let s = "GATCCAACTTCAAAGAGTCCTGGCCAGACATGATAAGATACATTGATGAGTTTGGACAAACCACAACTAGAATGCAGTGAAAAAAATGCTTTATTTGTGAAATTTGTGATGC" :: [DNA]+ in fmap symbol . _seq' <$> designPrimer False s 61 `shouldBe` Right "CACAACTAGAATGCAGTGAAAAAAATG"+ it "Sequence: GATCCAACTTCAAAGAGTCCTGGCCAGACATGATAAGATACATTGATGAGTTTGGACAAACCACAACTAGAATGCAGTGAAAAAAATGCTTTATTTGTGAAATTTGTGATGC; pos: 90" $+ let s = "GATCCAACTTCAAAGAGTCCTGGCCAGACATGATAAGATACATTGATGAGTTTGGACAAACCACAACTAGAATGCAGTGAAAAAAATGCTTTATTTGTGAAATTTGTGATGC" :: [DNA]+ in fmap symbol . _seq' <$> designPrimer False s 90 `shouldBe` Right "TTATTTGTGAAATTTGTGATGC"+ it "Sequence: GATCCAACTTCAAAGAGTCCTGGCCAGACATGATAAGATACATTGATGAGTTTGGACAAACCACAACTAGAATGCAGTGAAAAAAATGCTTTATTTGTGAAATTTGTGATGCTATTGCTTTATTTGTAACCATTATAAGCTGCAATAAACAAGTTTCAGGCACCGGGCTTGCGGGTCATGCAC; pos: 90" $+ let s = "GATCCAACTTCAAAGAGTCCTGGCCAGACATGATAAGATACATTGATGAGTTTGGACAAACCACAACTAGAATGCAGTGAAAAAAATGCTTTATTTGTGAAATTTGTGATGCTATTGCTTTATTTGTAACCATTATAAGCTGCAATAAACAAGTTTCAGGCACCGGGCTTGCGGGTCATGCAC" :: [DNA]+ in fmap symbol . _seq' <$> designPrimer False s 90 `shouldBe` Right "TTATTTGTGAAATTTGTGATGCTATTGC"+ it "Shouldn't find any primers. Sequence: GCAATAGCATCACAAATTTCACAAATAAGCAATAGCATCACAAATTTCACAAATAAGCAATAGCATCACAAATTTCACAAATAATTATTTGTGAAATTTGTGATGCTATTGCTTATTTGTGAAATTTGTGATGCTATTGCTTATTTGTGAAATTTGTGATGCTATTGC; pos: 0" $+ let s = "GCAATAGCATCACAAATTTCACAAATAAGCAATAGCATCACAAATTTCACAAATAAGCAATAGCATCACAAATTTCACAAATAATTATTTGTGAAATTTGTGATGCTATTGCTTATTTGTGAAATTTGTGATGCTATTGCTTATTTGTGAAATTTGTGATGCTATTGC" :: [DNA]+ in fmap symbol . _seq' <$> designPrimer False s 0 `shouldBe` Left "Bio.Tools.Sequence.Primers.Optimization: all primers designed from given position don't bind to target."+ it "Shouldn't find any primers. Sequence: ACAAATAATTATTTGTGAAATTTGACAAATAATTATTTGTGAAATTTGACAAATAATTATTTGTGAAATTTGCAAATTTCACAAATAATTATTTGTCAAATTTCACAAATAATTATTTGTCAAATTTCACAAATAATTATTTGT; pos: 61" $+ let s = "GCAATAGCATCACAAATTTCACAAATAAGCAATAGCATCACAAATTTCACAAATAAGCAATAGCATCACAAATTTCACAAATAATTATTTGTGAAATTTGTGATGCTATTGCTTATTTGTGAAATTTGTGATGCTATTGCTTATTTGTGAAATTTGTGATGCTATTGC" :: [DNA]+ in fmap symbol . _seq' <$> designPrimer False s 61 `shouldBe` Left "Bio.Tools.Sequence.Primers.Optimization: all primers designed from given position don't bind to target."+ let wholePlasmid = "CCTGCAGGCAGCTGCGCGCTCGCTCGCTCACTGAGGCCGCCCGGGCGTCGGGCGACCTTTGGTCGCCCGGCCTCAGTGAGCGAGCGAGCGCGCAGAGAGGGAGTGGCCAACTCCATCACTAGGGGTTCCTGCGGCCCGACATTCCGGAGGTACCTCTAGATTGGCAAACAGCTATTATGGGTATTATGGGTGATCTCCAGATGGCTAAACTTTTAAATCATGAATGAAGTAGATATTACCAAATTGCTTTTTCAGCATCCATTTAGATAATCATGTTTTTTGCCTTTAATCTGTTAATGTAGTGAATTACAGAAATACATTTCCTAAATCATTACATCCCCCAAATCGTTAATCTGCTAAAGTACATCTCTGGCTCAAACAAGACTGGTTGTGCATCTCAATTAGTCAGCAACCATAGTCCCGCCCCTAACTCCGCCCATCCCGCCCCTAACTCCGCCCAGTTCCGCCCATTCTCCGCCCCATCGCTGACTAATTTTTTTTATTTATGCAGAGGCCGAGGCCGCCTCGGCCTCTGAGCTATTCCAGAAGTAGTGAGGAGGCTTTTTTGGAGGCCTAGGCTTTTGCAAACGTAACTATAACGGTCCTAAGGTAGCGAAATGGTGAGCAAGGGCGAGGAGCTGTTCACCGGGGTGGTGCCCATCCTGGTCGAGCTGGACGGCGACGTAAACGGCCACAAGTTCAGCGTGTCCGGCGAGGGCGAGGGCGATGCCACCTACGGCAAGCTGACCCTGAAGTTCATCTGCACCACCGGCAAGCTGACCGTGCCCTGGCCCACCCTCGTGACCACCCTGACCTACGGCGTGCAGTGCTTCAGCCGCTACCCCGACCACATGAAGCAGCACGACTTCTTCAAGTCCGCCATGCCCGAAGGCTACGTCCAGGAGCGCACCATCTTCTTCAAGGACGACGGCAACTACAAGACCCGCGCCGAGGTGAAGTTCGAGGGCGACACCCTGGTGAACCGCATCGAGCTGAAGGGCATCGACTTCAAGGAGGACGGCAACATCCTGGGGCACAAGCTGGAGTACAACTACAACAGCCACAACGTCTATATCATGGCCGACAAGCAGAAGAACGGCATCAAGGTGAACTTCAAGATCCGCCACAACATCGAGGACGGCAGCGTGCAGCTCGCCGACCACTACCAGCAGAACACCCCCATCGGCGACGGCCCCGTGCTGCTGCCCGACAACCACTACCTGAGCACCCAGTCCGCCCTGAGCAAAGACCCCAACGAGAAGCGCGATCACATGGTCCTGCTGGAGTTCGTGACCGCCGCCGGGATCACTCTCGGCATGGACGAGCTGTACAAGTGATAGGGATAACAGGGTAATGTCGACAAGCTTGCTAGCACGGGTGGCATCCCTGTGACCCCTCCCCAGTGCCTCTCCTGGCCCTGGAAGTTGCCACTCCAGTGCCCACCAGCCTTGTCCTAATAAAATTAAGTTGCATCATTTTGTCTGACTAGGTGTCCTTCTATAATATTATGGGGTGGAGGGGGGTGGTATGGAGCAAGGGGCAAGTTGGGAAGACAACCTGTAGGGCCTGCGGGGTCTATTGGGAACCAAGCTGGAGTGCAGTGGCACAATCTTGGCTCACTGCAATCTCCGCCTCCTGGGTTCAAGCGATTCTCCTGCCTCAGCCTCCCGAGTTGTTGGGATTCCAGGCATGCATGACCAGGCTCAGCTAATTTTTGTTTTTTTGGTAGAGACGGGGTTTCACCATATTGGCCAGGCTGGTCTCCAACTCCTAATCTCAGGTGATCTACCCACCTTGGCCTCCCAAATTGCTGGGATTACAGGCGTGAACCACTGCTCCCTTCCCTGTCCTTATGCATGGGCCCGTACGATCACTAGTGTACAGCGGCCGCAGGAACCCCTAGTGATGGAGTTGGCCACTCCCTCTCTGCGCGCTCGCTCGCTCACTGAGGCCGGGCGACCAAAGGTCGCCCGACGCCCGGGCTTTGCCCGGGCGGCCTCAGTGAGCGAGCGAGCGCGCAGCTGCCTGCAGGGGCGCCTGATGCGGTATTTTCTCCTTACGCATCTGTGCGGTATTTCACACCGCATACGTCAAAGCAACCATAGTACGCGCCCTGTAGCGGCGCATTAAGCGCGGCGGGTGTGGTGGTTACGCGCAGCGTGACCGCTACACTTGCCAGCGCCTTAGCGCCCGCTCCTTTCGCTTTCTTCCCTTCCTTTCTCGCCACGTTCGCCGGCTTTCCCCGTCAAGCTCTAAATCGGGGGCTCCCTTTAGGGTTCCGATTTAGTGCTTTACGGCACCTCGACCCCAAAAAACTTGATTTGGGTGATGGTTCACGTAGTGGGCCATCGCCCTGATAGACGGTTTTTCGCCCTTTGACGTTGGAGTCCACGTTCTTTAATAGTGGACTCTTGTTCCAAACTGGAACAACACTCAACTCTATCTCGGGCTATTCTTTTGATTTATAAGGGATTTTGCCGATTTCGGTCTATTGGTTAAAAAATGAGCTGATTTAACAAAAATTTAACGCGAATTTTAACAAAATATTAACGTTTACAATTTTATGGTGCACTCTCAGTACAATCTGCTCTGATGCCGCATAGTTAAGCCAGCCCCGACACCCGCCAACACCCGCTGACGCGCCCTGACGGGCTTGTCTGCTCCCGGCATCCGCTTACAGACAAGCTGTGACCGTCTCCGGGAGCTGCATGTGTCAGAGGTTTTCACCGTCATCACCGAAACGCGCGAGACGAAAGGGCCTCGTGATACGCCTATTTTTATAGGTTAATGTCATGATAATAATGGTTTCTTAGACGTCAGGTGGCACTTTTCGGGGAAATGTGCGCGGAACCCCTATTTGTTTATTTTTCTAAATACATTCAAATATGTATCCGCTCATGAGACAATAACCCTGATAAATGCTTCAATAATATTGAAAAAGGAAGAGTATGAGTATTCAACATTTCCGTGTCGCCCTTATTCCCTTTTTTGCGGCATTTTGCCTTCCTGTTTTTGCTCACCCAGAAACGCTGGTGAAAGTAAAAGATGCTGAAGATCAGTTGGGTGCACGAGTGGGTTACATCGAACTGGATCTCAACAGCGGTAAGATCCTTGAGAGTTTTCGCCCCGAAGAACGTTTTCCAATGATGAGCACTTTTAAAGTTCTGCTATGTGGCGCGGTATTATCCCGTATTGACGCCGGGCAAGAGCAACTCGGTCGCCGCATACACTATTCTCAGAATGACTTGGTTGAGTACTCACCAGTCACAGAAAAGCATCTTACGGATGGCATGACAGTAAGAGAATTATGCAGTGCTGCCATAACCATGAGTGATAACACTGCGGCCAACTTACTTCTGACAACGATCGGAGGACCGAAGGAGCTAACCGCTTTTTTGCACAACATGGGGGATCATGTAACTCGCCTTGATCGTTGGGAACCGGAGCTGAATGAAGCCATACCAAACGACGAGCGTGACACCACGATGCCTGTAGCAATGGCAACAACGTTGCGCAAACTATTAACTGGCGAACTACTTACTCTAGCTTCCCGGCAACAATTAATAGACTGGATGGAGGCGGATAAAGTTGCAGGACCACTTCTGCGCTCGGCCCTTCCGGCTGGCTGGTTTATTGCTGATAAATCTGGAGCCGGTGAGCGTGGGTCTCGCGGTATCATTGCAGCACTGGGGCCAGATGGTAAGCCCTCCCGTATCGTAGTTATCTACACGACGGGGAGTCAGGCAACTATGGATGAACGAAATAGACAGATCGCTGAGATAGGTGCCTCACTGATTAAGCATTGGTAACTGTCAGACCAAGTTTACTCATATATACTTTAGATTGATTTAAAACTTCATTTTTAATTTAAAAGGATCTAGGTGAAGATCCTTTTTGATAATCTCATGACCAAAATCCCTTAACGTGAGTTTTCGTTCCACTGAGCGTCAGACCCCGTAGAAAAGATCAAAGGATCTTCTTGAGATCCTTTTTTTCTGCGCGTAATCTGCTGCTTGCAAACAAAAAAACCACCGCTACCAGCGGTGGTTTGTTTGCCGGATCAAGAGCTACCAACTCTTTTTCCGAAGGTAACTGGCTTCAGCAGAGCGCAGATACCAAATACTGTTCTTCTAGTGTAGCCGTAGTTAGGCCACCACTTCAAGAACTCTGTAGCACCGCCTACATACCTCGCTCTGCTAATCCTGTTACCAGTGGCTGCTGCCAGTGGCGATAAGTCGTGTCTTACCGGGTTGGACTCAAGACGATAGTTACCGGATAAGGCGCAGCGGTCGGGCTGAACGGGGGGTTCGTGCACACAGCCCAGCTTGGAGCGAACGACCTACACCGAACTGAGATACCTACAGCGTGAGCTATGAGAAAGCGCCACGCTTCCCGAAGGGAGAAAGGCGGACAGGTATCCGGTAAGCGGCAGGGTCGGAACAGGAGAGCGCACGAGGGAGCTTCCAGGGGGAAACGCCTGGTATCTTTATAGTCCTGTCGGGTTTCGCCACCTCTGACTTGAGCGTCGATTTTTGTGATGCTCGTCAGGGGGGCGGAGCCTATGGAAAAACGCCAGCAACGCGGCCTTTTTACGGTTCCTGGCCTTTTGCTGGCCTTTTGCTCACATGT"+ it "Shouldn't find any primers; Sequence: whole plasmid; pos: 5" $+ fmap symbol . _seq' <$> designPrimer True wholePlasmid 5 `shouldBe` Left "Bio.Tools.Sequence.Primers.Optimization: all primers designed from given position don't bind to target."+ it "Sequence: whole plasmid; pos: 1781" $+ fmap symbol . _seq' <$> designPrimer True wholePlasmid 1781 `shouldBe` Right "TGATCTACCCACCTTGGCCTCCC"+ it "Sequence: whole plasmid; pos: 4628 (last index in linear sequence)" $+ fmap symbol . _seq' <$> designPrimer True wholePlasmid 4628 `shouldBe` Right "TCCTGCAGGCAGCTGCGCGC"+ it "Should fail, because index is out of range; Sequence: whole plasmid; pos: 5000" $+ fmap symbol . _seq' <$> designPrimer True wholePlasmid 5000 `shouldBe` Left "Bio.Tools.Sequence.Primers.Optimization: given position is out of range."+ let wholePlasmid2 = "GGACGTCCGTCGACGCGCGAGCGAGCGAGTGACTCCGGCGGGCCCGCAGCCCGCTGGAAACCAGCGGGCCGGAGTCACTCGCTCGCTCGCGCGTCTCTCCCTCACCGGTTGAGGTAGTGATCCCCAAGGATGCGCAGATCAATAATTATCATTAGTTAATGCCCCAGTAATCAAGTATCGGGTATATACCTCAAGGCGCAATGTATTGAATGCCATTTACCGGGCGGACCGACTGGCGGGTTGCTGGGGGCGGGTAACTGCAGTTATTACTGCATACAAGGGTATCATTGCGGTTATCCCTGAAAGGTAACTGCAGTTACCCACCTCATAAATGCCATTTGACGGGTGAACCGTCATGTAGTTCACATAGTATACGGTTCATGCGGGGGATAACTGCAGTTACTGCCATTTACCGGGCGGACCGTAATACGGGTCATGTACTGGAATACCCTGAAAGGATGAACCGTCATGTAGATGCATAATCAGTAGCGATAATGGTACCACTACGCCAAAACCGTCATGTAGTTACCCGCACCTATCGCCAAACTGAGTGCCCCTAAAGGTTCAGAGGTGGGGTAACTGCAGTTACCCTCAAACAAAACCGTGGTTTTAGTTGCCCTGAAAGGTTTTACAGCATTGTTGAGGCGGGGTAACTGCGTTTACCCGCCATCCGCACATGCCACCCTCCAGATATATTCGTCTCGAGCAAATCACTTGGCAGTCTAGCGGACCTCTGCGGTAGGTGCGACAAAACTGGAGGTATCTTCTGTGGCCCTGGCTAGGTCGGAGGTACGTCGCGCACTTGTACTAGTACCGGCTCTCGGGACCGGACTAGTGGTAGACGGACGACCCGATGGACGACTCGCGGCTCACGTGGCACAAGGACCTGGTGCTCTTGCGGTTGTTCTAGGACTTGGCCGGGTTCTCTATGTTGTCGCCGTTCGACCTCCTCAAGCACGTCCCGTTGGACCTCTCCCTCACGTACCTCCTCTTCACGTCGAAGCTCCTCCGGTCCCTTCACAAGCTCTTGTGGCTCGCCTGGTGGCTCAAGACCTTCGTCATGCACCTGCCGCTGGTCACGCTCTCGTTGGGAACGGACTTGCCGCCGTCGACGTTCCTGCTGTAGTTGTCGATGCTCACGACCACGGGAAAGCCGAAGCTCCCGTTCTTGACGCTCGACCTGCACTGGACGTTGTAGTTCTTGCCGGCGACGCTCGTCAAGACGTTCTTGTCGCGGCTGTTGTTTCACCACACATCGACGTGGCTCCCGATGTCTGACCGGCTCTTGGTCTTCTCGACGCTCGGGCGGCACGGGAAGGGGACGCCGTCTCACTCGCACAGGGTCTGGTCGTTCGACTGGTCTCGGCTCTGGCACAAGGGGCTGCACCTGATGCACTTATCGTGGCTCCGGCTCTGGTAGGACCTGTTGTAGTGGGTCTCGTGGGTCAGGAAGTTGCTGAAGTGGTCTCAACACCCGCCGCTCCTGCGGTTCGGGCCGGTCAAGGGGACCGTCCACCACGACTTGCCGTTTCACCTACGGAAGACGCCGCCGTCGTAGCACTTGCTCTTCACCTAGCACTGTCGGCGGGTGACGCACCTCTGGCCGCACTTCTAGTGGCACCACCGGCCGCTTGTGTTATAGCTCCTCTGGCTCGTGTGGCTCGTCTTCGCCTTGCAGTAGGCCTAATAGGGGGTGGTGTTGATGTTGCGGCGGTAGTTGTTCATGTTGGTGCTGTAGCGGGACGACCTCGACCTGCTCGGAGACCACGACTTATCGATGCACTGGGGGTAGACGTAGCGGCTGTTCCTCATGTGGTTGTAGAAGGACTTCAAGCCGTCGCCGATGCACAGGCCGACCCCGTCTCACAAGGTGTTCCCGTCTTCGCGGGACCACGACGTCATGGACTCTCACGGGGACCACCTGTCTCGGTGGACGGACGAATCGTGGTTCAAGTGGTAGATGTTGTTGTACAAGACGCGGCCGAAGGTGCTCCCGCCGTCTCTGTCGACGGTCCCGCTGTCGCCGCCTGGGGTGCACTGGCTTCACCTCCCGTGGTCGAAGGACTGGCCGTAGTAGTCGACCCCGCTCCTCACGCGGTACTTCCCGTTCATGCCGTAGATGTGGTTTCACTCGGCCATGCACTTGACCTAGTTCCTCTTTTGGTTCGACTGGACTGACTTAAGCAGCTGTTAGTTGGAGACCTAATGTTTTAAACACTTTCTAACTGACCATAAGAATTGATACAACGAGGAAAATGCGATACACCTATGCGACGAAATTACGGAAACATAGTACGATAACGAAGGGCATACCGAAAGTAAAAGAGGAGGAACATATTTAGGACCAACGACAGAGAAATACTCCTCAACACCGGGCAACAGTCCGTTGCACCGCACCACACGTGACACAAACGACTGCGTTGGGGGTGACCAACCCCGTAACGGTGGTGGACAGTCGAGGAAAGGCCCTGAAAGCGAAAGGGGGAGGGATAACGGTGCCGCCTTGAGTAGCGGCGGACGGAACGGGCGACGACCTGTCCCCGAGCCGACAACCCGTGACTGTTAAGGCACCACAACAGCCCCTTTAGTAGCAGGAAAGGAACCGACGAGCGGACACAACGGTGGACCTAAGACGCGCCCTGCAGGAAGACGATGCAGGGAAGCCGGGAGTTAGGTCGCCTGGAAGGAAGGGCGCCGGACGACGGCCGAGACGCCGGAGAAGGCGCAGAAGCGGAAGCGGGAGTCTGCTCAGCCTAGAGGGAAACCCGGCGGAGGGGCGGACCGACGAGCTCTCTAGCCCACCGTAGGGACACTGGGGAGGGGTCACGGAGAGGACCGGGACCTTCAACGGTGAGGTCACGGGTGGTCGGAACAGGATTATTTTAATTCAACGTAGTAAAACAGACTGATCCACAGGAAGATATTATAATACCCCACCTCCCCCCACCATACCTCGTTCCCCGTTCAACCCTTCTGTTGGACATCCCGGACGCCCCAGATAACCCTTGGTTCGACCTCACGTCACCGTGTTAGAACCGAGTGACGTTAGAGGCGGAGGACCCAAGTTCGCTAAGAGGACGGAGTCGGAGGGCTCAACAACCCTAAGGTCCGTACGTACTGGTCCGAGTCGATTAAAAACAAAAAAACCATCTCTGCCCCAAAGTGGTATAACCGGTCCGACCAGAGGTTGAGGATTAGAGTCCACTAGATGGGTGGAACCGGAGGGTTTAACGACCCTAATGTCCGCACTTGGTGACGAGGGAAGGGACAGGAATCCTTGGGGATCACTACCTCAACCGGTGAGGGAGAGACGCGCGAGCGAGCGAGTGACTCCGGCCCGCTGGTTTCCAGCGGGCTGCGGGCCCGAAACGGGCCCGCCGGAGTCACTCGCTCGCTCGCGCGTCGACGGACGTCCCCGCGGACTACGCCATAAAAGAGGAATGCGTAGACACGCCATAAAGTGTGGCGTATGCAGTTTCGTTGGTATCATGCGCGGGACATCGCCGCGTAATTCGCGCCGCCCACACCACCAATGCGCGTCGCACTGGCGATGTGAACGGTCGCGGAATCGCGGGCGAGGAAAGCGAAAGAAGGGAAGGAAAGAGCGGTGCAAGCGGCCGAAAGGGGCAGTTCGAGATTTAGCCCCCGAGGGAAATCCCAAGGCTAAATCACGAAATGCCGTGGAGCTGGGGTTTTTTGAACTAAACCCACTACCAAGTGCATCACCCGGTAGCGGGACTATCTGCCAAAAAGCGGGAAACTGCAACCTCAGGTGCAAGAAATTATCACCTGAGAACAAGGTTTGACCTTGTTGTGAGTTGAGATAGAGCCCGATAAGAAAACTAAATATTCCCTAAAACGGCTAAAGCCAGATAACCAATTTTTTACTCGACTAAATTGTTTTTAAATTGCGCTTAAAATTGTTTTATAATTGCAAATGTTAAAATACCACGTGAGAGTCATGTTAGACGAGACTACGGCGTATCAATTCGGTCGGGGCTGTGGGCGGTTGTGGGCGACTGCGCGGGACTGCCCGAACAGACGAGGGCCGTAGGCGAATGTCTGTTCGACACTGGCAGAGGCCCTCGACGTACACAGTCTCCAAAAGTGGCAGTAGTGGCTTTGCGCGCTCTGCTTTCCCGGAGCACTATGCGGATAAAAATATCCAATTACAGTACTATTATTACCAAAGAATCTGCAGTCCACCGTGAAAAGCCCCTTTACACGCGCCTTGGGGATAAACAAATAAAAAGATTTATGTAAGTTTATACATAGGCGAGTACTCTGTTATTGGGACTATTTACGAAGTTATTATAACTTTTTCCTTCTCATACTCATAAGTTGTAAAGGCACAGCGGGAATAAGGGAAAAAACGCCGTAAAACGGAAGGACAAAAACGAGTGGGTCTTTGCGACCACTTTCATTTTCTACGACTTCTAGTCAACCCACGTGCTCACCCAATGTAGCTTGACCTAGAGTTGTCGCCATTCTAGGAACTCTCAAAAGCGGGGCTTCTTGCAAAAGGTTACTACTCGTGAAAATTTCAAGACGATACACCGCGCCATAATAGGGCATAACTGCGGCCCGTTCTCGTTGAGCCAGCGGCGTATGTGATAAGAGTCTTACTGAACCAACTCATGAGTGGTCAGTGTCTTTTCGTAGAATGCCTACCGTACTGTCATTCTCTTAATACGTCACGACGGTATTGGTACTCACTATTGTGACGCCGGTTGAATGAAGACTGTTGCTAGCCTCCTGGCTTCCTCGATTGGCGAAAAAACGTGTTGTACCCCCTAGTACATTGAGCGGAACTAGCAACCCTTGGCCTCGACTTACTTCGGTATGGTTTGCTGCTCGCACTGTGGTGCTACGGACATCGTTACCGTTGTTGCAACGCGTTTGATAATTGACCGCTTGATGAATGAGATCGAAGGGCCGTTGTTAATTATCTGACCTACCTCCGCCTATTTCAACGTCCTGGTGAAGACGCGAGCCGGGAAGGCCGACCGACCAAATAACGACTATTTAGACCTCGGCCACTCGCACCCAGAGCGCCATAGTAACGTCGTGACCCCGGTCTACCATTCGGGAGGGCATAGCATCAATAGATGTGCTGCCCCTCAGTCCGTTGATACCTACTTGCTTTATCTGTCTAGCGACTCTATCCACGGAGTGACTAATTCGTAACCATTGACAGTCTGGTTCAAATGAGTATATATGAAATCTAACTAAATTTTGAAGTAAAAATTAAATTTTCCTAGATCCACTTCTAGGAAAAACTATTAGAGTACTGGTTTTAGGGAATTGCACTCAAAAGCAAGGTGACTCGCAGTCTGGGGCATCTTTTCTAGTTTCCTAGAAGAACTCTAGGAAAAAAAGACGCGCATTAGACGACGAACGTTTGTTTTTTTGGTGGCGATGGTCGCCACCAAACAAACGGCCTAGTTCTCGATGGTTGAGAAAAAGGCTTCCATTGACCGAAGTCGTCTCGCGTCTATGGTTTATGACAAGAAGATCACATCGGCATCAATCCGGTGGTGAAGTTCTTGAGACATCGTGGCGGATGTATGGAGCGAGACGATTAGGACAATGGTCACCGACGACGGTCACCGCTATTCAGCACAGAATGGCCCAACCTGAGTTCTGCTATCAATGGCCTATTCCGCGTCGCCAGCCCGACTTGCCCCCCAAGCACGTGTGTCGGGTCGAACCTCGCTTGCTGGATGTGGCTTGACTCTATGGATGTCGCACTCGATACTCTTTCGCGGTGCGAAGGGCTTCCCTCTTTCCGCCTGTCCATAGGCCATTCGCCGTCCCAGCCTTGTCCTCTCGCGTGCTCCCTCGAAGGTCCCCCTTTGCGGACCATAGAAATATCAGGACAGCCCAAAGCGGTGGAGACTGAACTCGCAGCTAAAAACACTACGAGCAGTCCCCCCGCCTCGGATACCTTTTTGCGGTCGTTGCGCCGGAAAAATGCCAAGGACCGGAAAACGACCGGAAAACGAGTGTACA"+ it "Shouldn't find any primers; Sequence: whole plasmid; pos: 0" $+ fmap symbol . _seq' <$> designPrimer True wholePlasmid2 0 `shouldBe` Left "Bio.Tools.Sequence.Primers.Optimization: all primers designed from given position don't bind to target."+ it "Shouldn't find any primers; Sequence: whole plasmid; pos: 7" $+ fmap symbol . _seq' <$> designPrimer True wholePlasmid2 7 `shouldBe` Left "Bio.Tools.Sequence.Primers.Optimization: all primers designed from given position don't bind to target."+ it "Sequence: whole plasmid; pos: 10" $+ fmap symbol . _seq' <$> designPrimer True wholePlasmid2 10 `shouldBe` Right "CGACGCGCGAGCGAGCG"
+ test/SpecViennaRNA.hs view
@@ -0,0 +1,34 @@+module SpecViennaRNA where++import Bio.NucleicAcid.Nucleotide (DNA)+import Bio.Tools.Sequence.ViennaRNA.Cofold (cofold)+import Bio.Tools.Sequence.ViennaRNA.Fold (fold)+import Test.Hspec++foldSpec :: Spec+foldSpec = it "should work" $ do+ let (energy1, structure1) = fold 37 ("CTGGATCGCAATGACGCTCTTAGGTCTCGT" :: [DNA])+ energy1 `shouldBe` -2.9+ structure1 `shouldBe` "..(((((((......)).....)))))..."++ let (energy2, structure2) = fold 37 ("CTGGATCGCAATGGGTCTCGT" :: [DNA])+ energy2 `shouldBe` -4.4+ structure2 `shouldBe` "..(((((......)))))..."++cofoldSpec :: Spec+cofoldSpec = it "should work" $ do+ let (energy1, structure1) = cofold 37 ("CAAGTACAGTTACAAGAAAGTGGAGGAGGATTAGTACAACCGGGAGGAAGTCTCAG" :: [DNA], "TCTGAATCCACTCGCAGCACAGGAGAGTCTGAGACTTCCTCCCGGTTGTACTAATC" :: [DNA])+ energy1 `shouldBe` -56.40+ structure1 `shouldBe` "......((.((......)).))......((((((((((((((((((((((((((((.........((((...........))))))))))))))))))))))))))))))))"++ let (energy2, structure2) = cofold 37 ("TCTGAATCCACTCGCAGCACAGGAGAGTCTGAGACTTCCTCCCGGTTGTACTAATC" :: [DNA], "ACTCTCCTGTGCTGCGAGTGGATTCAGATTCAGTAACTACGCGATGAGTTGGGTCC" :: [DNA])+ energy2 `shouldBe` -62.30+ structure2 `shouldBe` "((((((((((((((((((((((((((((...((((..((....))..)).))....))))))))))))))))))))))))))))..........((.(((....))).)).."++foldTest :: Spec+foldTest = describe "Predict RNA secondary structure and calculate energy for one sequence" $ do+ foldSpec++cofoldTest :: Spec+cofoldTest = describe "Predict RNA secondary structure and calculate energy for two sequences" $ do+ cofoldSpec