diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -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. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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. + + diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/cobot-tools.cabal b/cobot-tools.cabal new file mode 100644 --- /dev/null +++ b/cobot-tools.cabal @@ -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 +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 diff --git a/src/Bio/Tools/Sequence/Primers/Constants.hs b/src/Bio/Tools/Sequence/Primers/Constants.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Tools/Sequence/Primers/Constants.hs @@ -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 diff --git a/src/Bio/Tools/Sequence/Primers/Optimization.hs b/src/Bio/Tools/Sequence/Primers/Optimization.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Tools/Sequence/Primers/Optimization.hs @@ -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 diff --git a/src/Bio/Tools/Sequence/Primers/Properties.hs b/src/Bio/Tools/Sequence/Primers/Properties.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Tools/Sequence/Primers/Properties.hs @@ -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 diff --git a/src/Bio/Tools/Sequence/Primers/Types.hs b/src/Bio/Tools/Sequence/Primers/Types.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Tools/Sequence/Primers/Types.hs @@ -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 diff --git a/src/Bio/Tools/Sequence/ViennaRNA/Cofold.hs b/src/Bio/Tools/Sequence/ViennaRNA/Cofold.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Tools/Sequence/ViennaRNA/Cofold.hs @@ -0,0 +1,5 @@ +module Bio.Tools.Sequence.ViennaRNA.Cofold + ( cofold + ) where + +import Bio.Tools.Sequence.ViennaRNA.Internal.Cofold (cofold) diff --git a/src/Bio/Tools/Sequence/ViennaRNA/Fold.hs b/src/Bio/Tools/Sequence/ViennaRNA/Fold.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Tools/Sequence/ViennaRNA/Fold.hs @@ -0,0 +1,5 @@ +module Bio.Tools.Sequence.ViennaRNA.Fold + ( fold + ) where + +import Bio.Tools.Sequence.ViennaRNA.Internal.Fold (fold) diff --git a/src/Bio/Tools/Sequence/ViennaRNA/Internal/Cofold.hs b/src/Bio/Tools/Sequence/ViennaRNA/Internal/Cofold.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Tools/Sequence/ViennaRNA/Internal/Cofold.hs @@ -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) diff --git a/src/Bio/Tools/Sequence/ViennaRNA/Internal/Fold.hs b/src/Bio/Tools/Sequence/ViennaRNA/Internal/Fold.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Tools/Sequence/ViennaRNA/Internal/Fold.hs @@ -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 <$>) diff --git a/src/Bio/Tools/Sequence/ViennaRNA/Internal/RNALike.hs b/src/Bio/Tools/Sequence/ViennaRNA/Internal/RNALike.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Tools/Sequence/ViennaRNA/Internal/RNALike.hs @@ -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 diff --git a/src/Bio/Tools/Sequence/ViennaRNA/Native/vienna_rna_wrapper.c b/src/Bio/Tools/Sequence/ViennaRNA/Native/vienna_rna_wrapper.c new file mode 100644 --- /dev/null +++ b/src/Bio/Tools/Sequence/ViennaRNA/Native/vienna_rna_wrapper.c @@ -0,0 +1,47 @@ +#include +#include +#include + +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; +} diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -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 diff --git a/test/SpecPrimers.hs b/test/SpecPrimers.hs new file mode 100644 --- /dev/null +++ b/test/SpecPrimers.hs @@ -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" diff --git a/test/SpecViennaRNA.hs b/test/SpecViennaRNA.hs new file mode 100644 --- /dev/null +++ b/test/SpecViennaRNA.hs @@ -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