diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -2,6 +2,10 @@ ## [Unreleased] +## [0.1.1.0] - 2019-11-13 +### Added +- Codon optimization tool + ## [0.1.0.1] - 2019-07-10 ### Added - `extra-lib-dirs` parameter to `stack.yaml`. diff --git a/cobot-tools.cabal b/cobot-tools.cabal --- a/cobot-tools.cabal +++ b/cobot-tools.cabal @@ -1,11 +1,13 @@ --- This file has been generated from package.yaml by hpack version 0.28.2. +cabal-version: 1.12 + +-- This file has been generated from package.yaml by hpack version 0.31.2. -- -- see: https://github.com/sol/hpack -- --- hash: 62f7da4d659cbc72389c4ea8d12693721e66abd1462ef2f2493c83c1601bf809 +-- hash: 1d89fad89173ca2b6e911afc77c4c881a041f3b7fee7a1da594289e32d6d1b22 name: cobot-tools -version: 0.1.0.1 +version: 0.1.1.0 synopsis: Biological data file formats and IO description: Please see the README on GitHub at category: Bio @@ -17,10 +19,9 @@ license: BSD3 license-file: LICENSE build-type: Simple -cabal-version: >= 1.10 extra-source-files: - ChangeLog.md README.md + ChangeLog.md source-repository head type: git @@ -28,6 +29,8 @@ library exposed-modules: + Bio.Tools.Sequence.CodonOptimization.Algo + Bio.Tools.Sequence.CodonOptimization.Type Bio.Tools.Sequence.Primers.Constants Bio.Tools.Sequence.Primers.Optimization Bio.Tools.Sequence.Primers.Properties @@ -49,12 +52,15 @@ build-depends: array >=0.5 && <0.6 , base >=4.7 && <5 + , bytestring , cobot , containers >=0.5.7.1 && <0.7 + , data-default , data-msgpack >=0.0.9 && <0.1 , deepseq >=1.4 && <1.5 , lens >=4.16 && <5.0 , mtl >=2.2.1 && <2.3.0 + , regex-tdfa , text default-language: Haskell2010 @@ -62,6 +68,7 @@ type: exitcode-stdio-1.0 main-is: Spec.hs other-modules: + SpecCodonOptimization SpecPrimers SpecViennaRNA Paths_cobot_tools @@ -73,14 +80,17 @@ QuickCheck >=2.9.2 && <2.13 , array >=0.5 && <0.6 , base >=4.7 && <5 + , bytestring , cobot , cobot-tools , containers >=0.5.7.1 && <0.7 + , data-default , 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 + , regex-tdfa , text default-language: Haskell2010 diff --git a/src/Bio/Tools/Sequence/CodonOptimization/Algo.hs b/src/Bio/Tools/Sequence/CodonOptimization/Algo.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Tools/Sequence/CodonOptimization/Algo.hs @@ -0,0 +1,153 @@ +module Bio.Tools.Sequence.CodonOptimization.Algo + ( optimizeAA + , optimizeDNA + , score + , scoreCmp + ) where + +import Bio.NucleicAcid.Nucleotide (symbol) +import Bio.NucleicAcid.Nucleotide.Type (DNA (..)) +import Bio.Protein.AminoAcid.Type (AA (..)) +import Bio.Tools.Sequence.CodonOptimization.Type (CodonConfig (..), CodonScoreConfig (..), + ak2Codon, + ak2MaxFrequCodon, + codon2ak, + codonFrequencies, + forbiddenRegexp, + standardTemperature) +import Bio.Tools.Sequence.ViennaRNA.Fold (fold) +import Data.List (maximumBy, take, foldl') +import Data.Map as Map (lookup) +import Data.Maybe (fromMaybe) +import Text.Regex.TDFA + +-- | 'optimizeDNA' function does translation from [DNA] to [AA] and then calls 'optimizeAA' +optimizeDNA :: CodonConfig -- ^ Config data object. Contains main parameters of codon-optimization and all parameters for scoring function + -> [DNA] -- ^ Initial, not optimized nucleotide sequence + -> [DNA] -- ^ Result, optimized nucleotide sequence +optimizeDNA cfg dna = optimizeAA cfg (translate dna) + where + translate :: [DNA] -> [AA] + translate [] = [] + translate dnaSeq = + case Map.lookup (take 3 dnaSeq) codon2ak of + Just ak -> ak : translate (drop 3 dnaSeq) + _ -> error $ "Unknown codon: " ++ show (take 3 dnaSeq) + +-- | 'optimizeAA' function does codon-optimisation for incoming amino-acid sequence. +-- Incoming amino-acid sequence transformed to nucleotide sequence and optimized used the codon-optimization algorithm. +-- Algorithm described here doi: 10.1007/s11693-010-9062-3 +optimizeAA :: CodonConfig -- ^ Config data object. Contains main parameters of codon-optimization and all parameters for scoring function + -> [AA] -- ^ Initial, not optimized amino-acid sequence + -> [DNA] -- ^ Result, optimized nucleotide sequence +optimizeAA cfg@(CodonConfig initLen winLen _) aa = foldl' concatByScore initial variants + where + lenAA = length aa + variants = generateVariants (drop initLen aa) winLen + initial = concatMap (\ak -> maybe "" fst (Map.lookup ak ak2MaxFrequCodon)) (take initLen aa) + + -- | 'concatByScore' function gets maximum by score variable string and then concat it to result string + concatByScore :: [DNA] -- ^ initial string + -> [[DNA]] -- ^ list of variable string + -> [DNA] -- ^ result string + concatByScore result vars + | length result == 3 * (lenAA - winLen - 1) = result ++ maximumBy (scoreCmp cfg result) vars + | otherwise = result ++ take 3 (maximumBy (scoreCmp cfg result) vars) + + +-- | 'generateVariants' function generates list of all possible variants of nucleotide sequence for amino-acid sequence. +-- It is just recursive execution of 'windowVariants' for all amino-acid sequence. +-- Example: generateVariants [PHE,TRP,GLU,MET] 2 => [[[DT,DT,DT,DT,DG,DG,DG,DA,DA], [DT,DT,DC,DT,DG,DG,DG,DA,DA], [DT,DT,DT,DT,DG,DG,DG,DA,DG], [DT,DT,DC,DT,DG,DG,DG,DA,DG]], +-- [[DT,DG,DG,DG,DA,DA,DA,DT,DG], [DT,DG,DG,DG,DA,DG,DA,DT,DG]]] +-- Returns empty list in case of empty incoming string +generateVariants :: [AA] -- ^ amino-acid sequence + -> Int -- ^ length of window. means how much amino-acid from the right side side will be taken during scoring variant for one codon + -> [[[DNA]]] -- ^ result list. for each amino-acid now there is list of all variants for nucleotide sequence +generateVariants [] _ = [] +generateVariants aa winLen + | length aa == winLen = [] + | otherwise = windowVariants aa winLen : generateVariants (drop 1 aa) winLen + +-- | 'windowVariants' function generates list of all possible variants of nucleotide sequence for amino-acid window. +-- Example: windowVariants [PHE,TRP,GLU,MET] 2 => [[DT,DT,DT,DT,DG,DG,DG,DA,DA], [DT,DT,DC,DT,DG,DG,DG,DA,DA], [DT,DT,DT,DT,DG,DG,DG,DA,DG], [DT,DT,DC,DT,DG,DG,DG,DA,DG]] +-- Returns empty list in case of empty incoming string +windowVariants :: [AA] -> Int -> [[DNA]] +windowVariants sequ winLen = map concat . mapM getCodons . take (winLen + 1) $ sequ + +-- | 'getCodons' function gets list of codons for amino-acid +-- Example: 'getCodons' PRO => [[DC,DC,DT], [DC,DC,DC], [DC,DC,DA], [DC,DC,DG]] +-- Returns empty list in case of unknown amino-acid +getCodons :: AA -> [[DNA]] +getCodons ak = fromMaybe [] (Map.lookup ak ak2Codon) + +-- | 'score' function gets scoring for incoming string. +-- Scoring function is a composite function of several scoring. More about scoring algorithm see here doi: 10.1007/s11693-010-9062-3 +score :: CodonConfig -- ^ Config data object. Contains main parameters of codon-optimization and all parameters for scoring function + -> [DNA] -- ^ nucleotide sequence to score + -> Double -- ^ result score value +score (CodonConfig _ winLen + (CodonScoreConfig codonUsageWeight gcWeight gcFactor gcWindow rnaFoldingWeight + rnaFoldingFactor rnaFoldingWindow forbiddenDNAWeight gcContentDesired)) nkSequ = + scoreCU + scoreGC - scoreMT - realToFrac scoreRNAFold + where + sequLen = length nkSequ + optimizedLen = sequLen - (winLen + 1) * 3 + scoreGC = + if optimizedLen < gcWindow - optimizedLen -- check if we have enough sequence for gcWindow + then 0 + else realToFrac gcWeight * gcScore (drop (sequLen - gcWindow) nkSequ) + scoreCU = realToFrac codonUsageWeight * codonUsage (drop optimizedLen nkSequ) + scoreMT = realToFrac forbiddenDNAWeight * motiveScore nkSequ + scoreRNAFold = scoreRnaf nkSequ + + -- | 'gcScore' function for the GC content. + -- It is negatively counted absolute difference between the desired GC content and the GC content of a test sequence. + gcScore :: [DNA] -> Double + gcScore sequ = -abs (gc / (at + gc) * 100 - realToFrac gcContentDesired) ** gcFactor + where + gc = realToFrac $ length $ filter (\s -> s == DC || s == DG) sequ + at = realToFrac $ length $ filter (\s -> s == DA || s == DT) sequ + + -- | 'codonUsage' function gets higest score for most frequently used codons + codonUsage :: [DNA] -> Double + codonUsage sequ = (cai ** (1 / codonCount)) * 100 + where + codonCount = realToFrac (length sequ) / 3 + cai = countWeight sequ 1.0 + + -- | 'countWeight' function is recursive counting weight for incoming string according to codon usage frequencies + countWeight :: [DNA] -> Double -> Double + countWeight [] acc = acc + countWeight windowSeq acc = countWeight (drop 3 windowSeq) (acc * fromMaybe 0 (countWeightMb windowSeq)) + + countWeightMb :: [DNA] -> Maybe Double + countWeightMb str = do + let codon = take 3 str + codonFreq <- Map.lookup codon codonFrequencies + ak <- Map.lookup codon codon2ak + (_, codonMaxFreq) <- Map.lookup ak ak2MaxFrequCodon + return $ codonFreq / codonMaxFreq + + -- | 'motiveScore' counts score for the occurrence of desired and unwanted DNA motifs. + motiveScore :: [DNA] -> Double + motiveScore sequ = + if any (drop (length sequ - 100) (map symbol sequ) =~) forbiddenRegexp + then 50 + else 0 + + + -- | 'scoreRnaf' counts energy of RNA folding + scoreRnaf :: [DNA] -> Int + scoreRnaf sequ = truncate $ rnaFoldingWeight * (abs result ** rnaFoldingFactor) + where + result = fst $ fold standardTemperature (drop (length sequ - rnaFoldingWindow) sequ) + +-- | 'scoreCmp' is compare function for two strings using 'score' function +scoreCmp :: CodonConfig -> [DNA] -> [DNA] -> [DNA] -> Ordering +scoreCmp cfg optimized str1 str2 + | score1 > score2 = GT + | score1 == score2 = EQ + | otherwise = LT + where + score1 = score cfg (optimized ++ str1) + score2 = score cfg (optimized ++ str2) diff --git a/src/Bio/Tools/Sequence/CodonOptimization/Type.hs b/src/Bio/Tools/Sequence/CodonOptimization/Type.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Tools/Sequence/CodonOptimization/Type.hs @@ -0,0 +1,251 @@ +module Bio.Tools.Sequence.CodonOptimization.Type + ( CodonConfig(..) + , CodonScoreConfig(..) + , ak2Codon + , codon2ak + , forbiddenRegexp + , codonFrequencies + , ak2MaxFrequCodon + , standardTemperature + ) where + +import Bio.NucleicAcid.Nucleotide.Type (DNA (..)) +import Bio.Protein.AminoAcid.Type (AA (..)) +import qualified Data.ByteString.Lazy as BSL (ByteString) +import Data.Default (Default (..), def) +import Data.Map as Map (Map, fromList) + +standardTemperature :: Double +standardTemperature = 37 + +-- | all parameters for codon optimization +data CodonConfig = + CodonConfig + { initLen :: Int -- ^ number of first ak from initial sequence, which will optimised without scoring function + , windowLen :: Int -- ^ length of variation window + , scoreConf :: CodonScoreConfig -- ^ parameters for scoring function + } deriving (Show) + +-- | all parameters for scoring function of codon optimization +data CodonScoreConfig = + CodonScoreConfig + { codonUsageWeight :: Double -- ^ Codon usage weight + , gcWeight :: Double -- ^ GC-content weight + , gcFactor :: Double -- ^ GC_score in the power of F_gc is used + , gcWindow :: Int -- ^ length of the window for GC-score calculation (bp) + , rnaFoldingWeight :: Float -- ^ Weight of the RNA folding score + , rnaFoldingFactor :: Float -- ^ RNA folding score in the power of F_rnaf is used + , rnaFoldingWindow :: Int -- ^ length of the window for RNA folding score calculation (bp) + , forbiddenDNAWeight :: Double -- ^ forbidden DNA motifs score weight + , gcContentDesired :: Int -- ^ desired gc content in percents + } deriving (Show) + +instance Default CodonScoreConfig where + def = CodonScoreConfig 1 0.5 1.4 40 0.001 2.6 100 1 43 + +instance Default CodonConfig where + def = CodonConfig 3 1 def + +forbiddenRegexp :: [BSL.ByteString] +forbiddenRegexp = + [ "ATTTA" + , "ATACTCCCCC" + , "CGATCG" + , "GGGGACTTTGCACTGGAACTTACAACACCCCAGCAAGGACGCG" + , "CCGGCGGGT" + , "TTTATAATTTCTTCTTCCAGAA" + , "CCGTGCTGGCGTCTG" + , "AATAAA.{10,30}CA{30,}(TCTG|TG.CT)" + , "(GGG|CCA)CGCCTATAAA(((C|T)(C|T)A.(T|A)(C|T)(C|T))|(TCA(G|T)T(T|C)))(A|G)G(A|T)(C|T)(G|A|C)" + , "CAGG" + , "(A|C)AGGT(A|G)AGT" + , "AATAAA" + , "GCC(A|G)CCATGG" + ] + +ak2Codon :: Map AA [[DNA]] +ak2Codon = + fromList + [ (PHE, [[DT, DT, DT], [DT, DT, DC]]) + , (TYR, [[DT, DA, DT], [DT, DA, DC]]) + , (CYS, [[DT, DG, DT], [DT, DG, DC]]) + , (TRP, [[DT, DG, DG]]) + , (LEU, [[DT, DT, DA], [DT, DT, DG], [DC, DT, DT], [DC, DT, DC], [DC, DT, DA], [DC, DT, DG]]) + , (PRO, [[DC, DC, DT], [DC, DC, DC], [DC, DC, DA], [DC, DC, DG]]) + , (HIS, [[DC, DA, DT], [DC, DA, DC]]) + , (GLN, [[DC, DA, DA], [DC, DA, DG]]) + , (ILE, [[DA, DT, DT], [DA, DT, DC], [DA, DT, DA]]) + , (MET, [[DA, DT, DG]]) + , (THR, [[DA, DC, DT], [DA, DC, DC], [DA, DC, DA], [DA, DC, DG]]) + , (ASN, [[DA, DA, DT], [DA, DA, DC]]) + , (LYS, [[DA, DA, DA], [DA, DA, DG]]) + , (SER, [[DT, DC, DT], [DT, DC, DC], [DT, DC, DA], [DT, DC, DG], [DA, DG, DT], [DA, DG, DC]]) + , (ARG, [[DC, DG, DT], [DC, DG, DC], [DC, DG, DA], [DC, DG, DG], [DA, DG, DA], [DA, DG, DG]]) + , (VAL, [[DG, DT, DT], [DG, DT, DC], [DG, DT, DA], [DG, DT, DG]]) + , (ALA, [[DG, DC, DT], [DG, DC, DC], [DG, DC, DA], [DG, DC, DG]]) + , (ASP, [[DG, DA, DT], [DG, DA, DC]]) + , (GLU, [[DG, DA, DA], [DG, DA, DG]]) + , (GLY, [[DG, DG, DT], [DG, DG, DC], [DG, DG, DA], [DG, DG, DG]]) + ] + + +ak2MaxFrequCodon :: Map AA ([DNA], Double) +ak2MaxFrequCodon = + fromList + [ (ALA, ([DG, DC, DC], 0.4)) + , (CYS, ([DT, DG, DC], 0.54)) + , (ASP, ([DG, DA, DC], 0.54)) + , (GLU, ([DG, DA, DG], 0.58)) + , (PHE, ([DT, DT, DC], 0.54)) + , (GLY, ([DG, DG, DC], 0.34)) + , (HIS, ([DC, DA, DC], 0.58)) + , (ILE, ([DA, DT, DC], 0.47)) + , (LYS, ([DA, DA, DG], 0.57)) + , (LEU, ([DC, DT, DG], 0.4)) + , (MET, ([DA, DT, DG], 1.0)) + , (ASN, ([DA, DA, DC], 0.53)) + , (PRO, ([DC, DC, DC], 0.32)) + , (GLN, ([DC, DA, DG], 0.73)) + , (ARG, ([DA, DG, DG], 0.21)) + , (SER, ([DA, DG, DC], 0.24)) + , (THR, ([DA, DC, DC], 0.36)) + , (VAL, ([DG, DT, DG], 0.46)) + , (TRP, ([DT, DG, DG], 1.0)) + , (TYR, ([DT, DA, DC], 0.56)) + ] + +codon2ak :: Map [DNA] AA +codon2ak = + fromList + [ ([DT, DT, DT], PHE) + , ([DT, DT, DC], PHE) + , ([DT, DT, DA], LEU) + , ([DT, DT, DG], LEU) + , ([DT, DC, DT], SER) + , ([DT, DC, DC], SER) + , ([DT, DC, DA], SER) + , ([DT, DC, DG], SER) + , ([DT, DA, DT], TYR) + , ([DT, DA, DC], TYR) + , ([DT, DG, DT], CYS) + , ([DT, DG, DC], CYS) + , ([DT, DG, DG], TRP) + , ([DC, DT, DT], LEU) + , ([DC, DT, DC], LEU) + , ([DC, DT, DA], LEU) + , ([DC, DT, DG], LEU) + , ([DC, DC, DT], PRO) + , ([DC, DC, DC], PRO) + , ([DC, DC, DA], PRO) + , ([DC, DC, DG], PRO) + , ([DC, DA, DT], HIS) + , ([DC, DA, DC], HIS) + , ([DC, DA, DA], GLN) + , ([DC, DA, DG], GLN) + , ([DC, DG, DT], ARG) + , ([DC, DG, DC], ARG) + , ([DC, DG, DA], ARG) + , ([DC, DG, DG], ARG) + , ([DA, DT, DT], ILE) + , ([DA, DT, DC], ILE) + , ([DA, DT, DA], ILE) + , ([DA, DT, DG], MET) + , ([DA, DC, DT], THR) + , ([DA, DC, DC], THR) + , ([DA, DC, DA], THR) + , ([DA, DC, DG], THR) + , ([DA, DA, DT], ASN) + , ([DA, DA, DC], ASN) + , ([DA, DA, DA], LYS) + , ([DA, DA, DG], LYS) + , ([DA, DG, DT], SER) + , ([DA, DG, DC], SER) + , ([DA, DG, DA], ARG) + , ([DA, DG, DG], ARG) + , ([DG, DT, DT], VAL) + , ([DG, DT, DC], VAL) + , ([DG, DT, DA], VAL) + , ([DG, DT, DG], VAL) + , ([DG, DC, DT], ALA) + , ([DG, DC, DC], ALA) + , ([DG, DC, DA], ALA) + , ([DG, DC, DG], ALA) + , ([DG, DA, DT], ASP) + , ([DG, DA, DC], ASP) + , ([DG, DA, DA], GLU) + , ([DG, DA, DG], GLU) + , ([DG, DG, DT], GLY) + , ([DG, DG, DC], GLY) + , ([DG, DG, DA], GLY) + , ([DG, DG, DG], GLY) + ] + +-- | taken from https://github.com/Edinburgh-Genome-Foundry/codon-usage-tables/blob/master/codon_usage_data/tables/h_sapiens_9606.csv +codonFrequencies :: Map [DNA] Double +codonFrequencies = + fromList + [ ([DT, DA, DA], 0.30) + , ([DT, DA, DG], 0.24) + , ([DT, DG, DA], 0.47) + , ([DG, DC, DA], 0.23) + , ([DG, DC, DC], 0.40) + , ([DG, DC, DG], 0.11) + , ([DG, DC, DT], 0.27) + , ([DT, DG, DC], 0.54) + , ([DT, DG, DT], 0.46) + , ([DG, DA, DC], 0.54) + , ([DG, DA, DT], 0.46) + , ([DG, DA, DA], 0.42) + , ([DG, DA, DG], 0.58) + , ([DT, DT, DC], 0.54) + , ([DT, DT, DT], 0.46) + , ([DG, DG, DA], 0.25) + , ([DG, DG, DC], 0.34) + , ([DG, DG, DG], 0.25) + , ([DG, DG, DT], 0.16) + , ([DC, DA, DC], 0.58) + , ([DC, DA, DT], 0.42) + , ([DA, DT, DA], 0.17) + , ([DA, DT, DC], 0.47) + , ([DA, DT, DT], 0.36) + , ([DA, DA, DA], 0.43) + , ([DA, DA, DG], 0.57) + , ([DC, DT, DA], 0.07) + , ([DC, DT, DC], 0.20) + , ([DC, DT, DG], 0.40) + , ([DC, DT, DT], 0.13) + , ([DT, DT, DA], 0.08) + , ([DT, DT, DG], 0.13) + , ([DA, DT, DG], 1.00) + , ([DA, DA, DC], 0.53) + , ([DA, DA, DT], 0.47) + , ([DC, DC, DA], 0.28) + , ([DC, DC, DC], 0.32) + , ([DC, DC, DG], 0.11) + , ([DC, DC, DT], 0.29) + , ([DC, DA, DA], 0.27) + , ([DC, DA, DG], 0.73) + , ([DA, DG, DA], 0.21) + , ([DA, DG, DG], 0.21) + , ([DC, DG, DA], 0.11) + , ([DC, DG, DC], 0.18) + , ([DC, DG, DG], 0.20) + , ([DC, DG, DT], 0.08) + , ([DA, DG, DC], 0.24) + , ([DA, DG, DT], 0.15) + , ([DT, DC, DA], 0.15) + , ([DT, DC, DC], 0.22) + , ([DT, DC, DG], 0.05) + , ([DT, DC, DT], 0.19) + , ([DA, DC, DA], 0.28) + , ([DA, DC, DC], 0.36) + , ([DA, DC, DG], 0.11) + , ([DA, DC, DT], 0.25) + , ([DG, DT, DA], 0.12) + , ([DG, DT, DC], 0.24) + , ([DG, DT, DG], 0.46) + , ([DG, DT, DT], 0.18) + , ([DT, DG, DG], 1.00) + , ([DT, DA, DC], 0.56) + , ([DT, DA, DT], 0.44) + ] diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -1,4 +1,5 @@ import SpecPrimers +import SpecCodonOptimization import SpecViennaRNA import System.IO import Test.Hspec @@ -13,3 +14,6 @@ -- ViennaRNA foldTest cofoldTest + + -- CodonOptimization + codonOptimizationSpec diff --git a/test/SpecCodonOptimization.hs b/test/SpecCodonOptimization.hs new file mode 100644 --- /dev/null +++ b/test/SpecCodonOptimization.hs @@ -0,0 +1,205 @@ +{-# LANGUAGE OverloadedStrings #-} + +module SpecCodonOptimization where + +import Bio.Protein.AminoAcid () +import Bio.Tools.Sequence.CodonOptimization.Algo (optimizeAA, optimizeDNA, score, + scoreCmp) +import Bio.Tools.Sequence.CodonOptimization.Type (CodonConfig (..), CodonScoreConfig (..)) +import Data.Default (def) +import Data.List (maximumBy, + minimumBy) +import Test.Hspec (Spec, describe, it, + shouldBe) + +scoreCfg :: CodonScoreConfig +scoreCfg = def + +conf :: CodonConfig +conf = def + +codonOptimizationSpec :: Spec +codonOptimizationSpec = + describe "Codon optimization spec" $ do + optimizeSequence + optimizeDNASequence + optimizeShortSequence + optimizeExtremelyShortSequence + optimizeSequenceWindow3 + optimizeSequenceInit5 + scoreComparing + scoreFun + scoreFunDifferentCodonUsageWeight + scoreFunDifferentGCWeight + scoreFunDifferentGCFactor + scoreFunDifferentGCWindow + scoreFunDifferentFoldingWeight + scoreFunDifferentFoldingFactor + scoreFunDifferentFoldingWindow + scoreFunWithForbiddenSeq + scoreFunDifferentForbiddenSeqWeight + scoreFunDifferentGCDesired + +scoreComparing :: Spec +scoreComparing = + describe "scoreComparing" $ + it "should correct compare by score" $ do + let optimized = "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGC" + let vars = ["AGTACTGGT","AGCACTGGT","TCGACAGGT","AGTACAGGT","AGCACAGGT","TCTACTGGC","TCCACTGGC","TCAACTGGC","TCGACTGGC","AGCACCGGC"] + let cfg = CodonConfig 3 3 (CodonScoreConfig 1 0.5 1.4 40 0.001 2.6 100 1 43) + let resMin = maximumBy (scoreCmp cfg optimized) vars + let resMax = minimumBy (scoreCmp cfg optimized) vars + resMax `shouldBe` "TCGACAGGT" + resMin `shouldBe` "AGCACCGGC" + +optimizeSequence :: Spec +optimizeSequence = + describe "optimizeSequence" $ + it "should correct optimize sequence" $ do + let res = optimizeAA conf "ASTKGPSVFPLAPSSKSTSGGTAALGCLVKDYFPEPVTVSWNSGALTSGVHTFPAVLQSSGLYSLSSVVTVPSSSLGTQTYICNVNHKPSNTKVDKKV" + res `shouldBe` "GCCAGCACCAAGGGCCCCAGCGTGTTTCCTCTGGCCCCTAGCAGCAAGAGCACCAGCGGCGGCACCGCCGCCCTGGGCTGCCTGGTGAAGGACTACTTCCCTGAGCCTGTGACCGTGAGCTGGAACAGCGGCGCCCTGACCAGCGGCGTGCACACCTTCCCTGCCGTGCTGCAGAGCAGCGGCCTGTACAGCCTGAGCAGCGTGGTGACCGTGCCTAGCAGCAGCCTGGGCACCCAGACCTACATCTGCAACGTGAACCACAAGCCTAGCAACACCAAGGTGGACAAGAAGGTG" + + +optimizeDNASequence :: Spec +optimizeDNASequence = + describe "optimizeDNASequence" $ + it "should correct optimize amino-acid sequence" $ do + let res = optimizeDNA conf "GCTAGTACCAAGGGCCCCAGCGTGTTTCCTCTGGCCCCTAGCAGCAAGAGCACCAGCGGCGGCACCGCCGCCCTGGGCTGCCTGGTGAAGGACTACTTCCCTGAGCCTGTGACCGTGAGCTGGAACAGCGGCGCCCTGACCAGCGGCGTGCACACCTTCCCTGCCGTGCTGCAGAGCAGCGGCCTGTACAGCCTGAGCAGCGTGGTGACCGTGCCTAGCAGCAGCCTGGGCACCCAGACCTACATCTGCAACGTGAACCACAAGCCTAGCAACACCAAGGTGGACAAGAAGGTG" + res `shouldBe` "GCCAGCACCAAGGGCCCCAGCGTGTTTCCTCTGGCCCCTAGCAGCAAGAGCACCAGCGGCGGCACCGCCGCCCTGGGCTGCCTGGTGAAGGACTACTTCCCTGAGCCTGTGACCGTGAGCTGGAACAGCGGCGCCCTGACCAGCGGCGTGCACACCTTCCCTGCCGTGCTGCAGAGCAGCGGCCTGTACAGCCTGAGCAGCGTGGTGACCGTGCCTAGCAGCAGCCTGGGCACCCAGACCTACATCTGCAACGTGAACCACAAGCCTAGCAACACCAAGGTGGACAAGAAGGTG" + +optimizeShortSequence :: Spec +optimizeShortSequence = + describe "optimizeShortSequence" $ + it "should correct optimize short sequence" $ do + let res = optimizeAA conf "METDTLLLWVLLLWVPGSTG" + res `shouldBe` "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGCAGCACCGGC" + +optimizeExtremelyShortSequence :: Spec +optimizeExtremelyShortSequence = + describe "optimizeExtremelyShortSequence" $ + it "should correct optimize extremely short sequence" $ do + let res = optimizeAA conf "METDTLL" + res `shouldBe` "ATGGAGACCGACACCCTGCTG" + +optimizeSequenceWindow3 :: Spec +optimizeSequenceWindow3 = + describe "optimizeSequenceWindow3" $ + it "should correct optimize sequence with window 3" $ do + let conf' = CodonConfig 3 3 scoreCfg + let res = optimizeAA conf' "METDTLLLWVLLLWVPGSTG" + res `shouldBe` "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGCTCTACCGGC" + +optimizeSequenceInit5 :: Spec +optimizeSequenceInit5 = + describe "optimizeSequenceInit5" $ + it "should correct optimize sequence with init param = 5" $ do + let conf' = CodonConfig 5 1 scoreCfg + let res = optimizeAA conf' "ASTKGPSVFPLAPSSKSTSGGTAALGCLVKDYFPEPVTVSWNSGALTSGVHTFPAVLQSSGLYSLSSVVTVPSSSLGTQTYICNVNHKPSNTKVDKKV" + res `shouldBe` "GCCAGCACCAAGGGCCCCAGCGTGTTTCCTCTGGCCCCTAGCAGCAAGAGCACCAGCGGCGGCACCGCCGCCCTGGGCTGCCTGGTGAAGGACTACTTCCCTGAGCCTGTGACCGTGAGCTGGAACAGCGGCGCCCTGACCAGCGGCGTGCACACCTTCCCTGCCGTGCTGCAGAGCAGCGGCCTGTACAGCCTGAGCAGCGTGGTGACCGTGCCTAGCAGCAGCCTGGGCACCCAGACCTACATCTGCAACGTGAACCACAAGCCTAGCAACACCAAGGTGGACAAGAAGGTG" + +scoreFun :: Spec +scoreFun = + describe "scoreFun" $ + it "should correct count score" $ do + score conf "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGCAGCACCGGC" `shouldBe` 40.88886033302965 + score conf "GCCAGCGGCGACAAGACCCACACCTGTCCT" `shouldBe` 45.91010855311742 + score conf "CCCTGCCCCGCCCCCGAGGCCGCCGGCGGCCCTAGCGTGTTCCTGTTCCCTCCTAAGCCTAAGGACACCCTGATGATCAGCAGAACCCCCGAGGTGACCTGCGTGGTGGTGGACGTGAGCCACGAGGACCCTGAGGTGAAGTTCAATTGGTACGTGGACGGCGTGGAGGTGCACAACGCCAAGACCAAGCCTAGAGAGGAGCAGTACAACAGCACCTACAGAGTGGTGAGCGTGCTGACCGTGCTGCACCAAGACTGGCTGAACGGCAAGGAGTACAAGTGCAAGGTGAGCAACAAGGCCCTGCCCGCCCCTATCGAGAAGACCATCAGCAAGGCCAAG" + `shouldBe` 66.00968914282373 + +scoreFunDifferentCodonUsageWeight :: Spec +scoreFunDifferentCodonUsageWeight = + describe "scoreFunDifferentCodonUsageWeight" $ + it "should correct count score with g_cu=9" $ do + let conf' = CodonConfig 3 1 (CodonScoreConfig 9 0.5 1.4 40 0.001 2.6 100 1 43) + score conf' "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGCAGCACCGGC" `shouldBe` 840.8888603330297 + score conf' "GCCAGCGGCGACAAGACCCACACCTGTCCT" `shouldBe` 748.8141376168824 + score conf' "CCCTGCCCCGCCCCCGAGGCCGCCGGCGGCCCTAGCGTGTTCCTGTTCCCTCCTAAGCCTAAGGACACCCTGATGATCAGCAGAACCCCCGAGGTGACCTGCGTGGTGGTGGACGTGAGCCACGAGGACCCTGAGGTGAAGTTCAATTGGTACGTGGACGGCGTGGAGGTGCACAACGCCAAGACCAAGCCTAGAGAGGAGCAGTACAACAGCACCTACAGAGTGGTGAGCGTGCTGACCGTGCTGCACCAAGACTGGCTGAACGGCAAGGAGTACAAGTGCAAGGTGAGCAACAAGGCCCTGCCCGCCCCTATCGAGAAGACCATCAGCAAGGCCAAG" + `shouldBe` 866.0096891428237 + +scoreFunDifferentGCWeight :: Spec +scoreFunDifferentGCWeight = + describe "scoreFunDifferentGCWeight" $ + it "should correct count score with g_gc=0.2" $ do + let conf' = CodonConfig 3 1 (CodonScoreConfig 1 0.2 1.4 40 0.001 2.6 100 1 43) + score conf' "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGCAGCACCGGC" `shouldBe` 75.15554413321186 + score conf' "GCCAGCGGCGACAAGACCCACACCTGTCCT" `shouldBe` 71.08184560102933 + score conf' "CCCTGCCCCGCCCCCGAGGCCGCCGGCGGCCCTAGCGTGTTCCTGTTCCCTCCTAAGCCTAAGGACACCCTGATGATCAGCAGAACCCCCGAGGTGACCTGCGTGGTGGTGGACGTGAGCCACGAGGACCCTGAGGTGAAGTTCAATTGGTACGTGGACGGCGTGGAGGTGCACAACGCCAAGACCAAGCCTAGAGAGGAGCAGTACAACAGCACCTACAGAGTGGTGAGCGTGCTGACCGTGCTGCACCAAGACTGGCTGAACGGCAAGGAGTACAAGTGCAAGGTGAGCAACAAGGCCCTGCCCGCCCCTATCGAGAAGACCATCAGCAAGGCCAAG" + `shouldBe` 85.20387565712949 + +scoreFunDifferentGCFactor :: Spec +scoreFunDifferentGCFactor = + describe "scoreFunDifferentGCFactor" $ + it "should correct count score with f_gc=1.7" $ do + let conf' = CodonConfig 3 1 (CodonScoreConfig 1 0.5 1.7 40 0.001 2.6 100 1 43) + score conf' "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGCAGCACCGGC" `shouldBe` -59.640364905662636 + score conf' "GCCAGCGGCGACAAGACCCACACCTGTCCT" `shouldBe` -20.530452529069947 + score conf' "CCCTGCCCCGCCCCCGAGGCCGCCGGCGGCCCTAGCGTGTTCCTGTTCCCTCCTAAGCCTAAGGACACCCTGATGATCAGCAGAACCCCCGAGGTGACCTGCGTGGTGGTGGACGTGAGCCACGAGGACCCTGAGGTGAAGTTCAATTGGTACGTGGACGGCGTGGAGGTGCACAACGCCAAGACCAAGCCTAGAGAGGAGCAGTACAACAGCACCTACAGAGTGGTGAGCGTGCTGACCGTGCTGCACCAAGACTGGCTGAACGGCAAGGAGTACAAGTGCAAGGTGAGCAACAAGGCCCTGCCCGCCCCTATCGAGAAGACCATCAGCAAGGCCAAG" + `shouldBe` 20.011809536823847 + +scoreFunDifferentGCWindow :: Spec +scoreFunDifferentGCWindow = + describe "scoreFunDifferentGCWindow" $ + it "should correct count score with w_gc=10" $ do + let conf' = CodonConfig 3 1 (CodonScoreConfig 1 0.5 1.4 10 0.001 2.6 100 1 43) + score conf' "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGCAGCACCGGC" `shouldBe` 19.575380904184584 + score conf' "GCCAGCGGCGACAAGACCCACACCTGTCCT" `shouldBe` 61.4633338693899 + score conf' "CCCTGCCCCGCCCCCGAGGCCGCCGGCGGCCCTAGCGTGTTCCTGTTCCCTCCTAAGCCTAAGGACACCCTGATGATCAGCAGAACCCCCGAGGTGACCTGCGTGGTGGTGGACGTGAGCCACGAGGACCCTGAGGTGAAGTTCAATTGGTACGTGGACGGCGTGGAGGTGCACAACGCCAAGACCAAGCCTAGAGAGGAGCAGTACAACAGCACCTACAGAGTGGTGAGCGTGCTGACCGTGCTGCACCAAGACTGGCTGAACGGCAAGGAGTACAAGTGCAAGGTGAGCAACAAGGCCCTGCCCGCCCCTATCGAGAAGACCATCAGCAAGGCCAAG" + `shouldBe` 71.60033023641928 + +scoreFunDifferentFoldingWeight :: Spec +scoreFunDifferentFoldingWeight = + describe "scoreFunDifferentFoldingWeight" $ + it "should correct count score with weight of folding = 0.03" $ do + let conf' = CodonConfig 3 1 (CodonScoreConfig 1 0.5 1.4 40 0.03 2.6 100 1 43) + score conf' "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGCAGCACCGGC" `shouldBe` -45.11113966697035 + score conf' "GCCAGCGGCGACAAGACCCACACCTGTCCT" `shouldBe` 43.91010855311742 + score conf' "CCCTGCCCCGCCCCCGAGGCCGCCGGCGGCCCTAGCGTGTTCCTGTTCCCTCCTAAGCCTAAGGACACCCTGATGATCAGCAGAACCCCCGAGGTGACCTGCGTGGTGGTGGACGTGAGCCACGAGGACCCTGAGGTGAAGTTCAATTGGTACGTGGACGGCGTGGAGGTGCACAACGCCAAGACCAAGCCTAGAGAGGAGCAGTACAACAGCACCTACAGAGTGGTGAGCGTGCTGACCGTGCTGCACCAAGACTGGCTGAACGGCAAGGAGTACAAGTGCAAGGTGAGCAACAAGGCCCTGCCCGCCCCTATCGAGAAGACCATCAGCAAGGCCAAG" + `shouldBe` -15.990310857176269 + +scoreFunDifferentFoldingFactor :: Spec +scoreFunDifferentFoldingFactor = + describe "scoreFunDifferentFoldingFactor" $ + it "should correct count score with folding factor = 4.999" $ do + let conf' = CodonConfig 3 1 (CodonScoreConfig 1 0.5 1.4 40 0.001 4.999 100 1 43) + score conf' "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGCAGCACCGGC" `shouldBe` -4644.11113966697 + score conf' "GCCAGCGGCGACAAGACCCACACCTGTCCT" `shouldBe` 42.91010855311742 + score conf' "CCCTGCCCCGCCCCCGAGGCCGCCGGCGGCCCTAGCGTGTTCCTGTTCCCTCCTAAGCCTAAGGACACCCTGATGATCAGCAGAACCCCCGAGGTGACCTGCGTGGTGGTGGACGTGAGCCACGAGGACCCTGAGGTGAAGTTCAATTGGTACGTGGACGGCGTGGAGGTGCACAACGCCAAGACCAAGCCTAGAGAGGAGCAGTACAACAGCACCTACAGAGTGGTGAGCGTGCTGACCGTGCTGCACCAAGACTGGCTGAACGGCAAGGAGTACAAGTGCAAGGTGAGCAACAAGGCCCTGCCCGCCCCTATCGAGAAGACCATCAGCAAGGCCAAG" + `shouldBe` -4200.990310857176 + +scoreFunDifferentFoldingWindow :: Spec +scoreFunDifferentFoldingWindow = + describe "scoreFunDifferentFoldingWindow" $ + it "should correct count score with folding window = 23" $ do + let conf' = CodonConfig 3 1 (CodonScoreConfig 1 0.5 1.4 40 0.001 2.6 23 1 43) + score conf' "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGCAGCACCGGC" `shouldBe` 42.88886033302965 + score conf' "GCCAGCGGCGACAAGACCCACACCTGTCCT" `shouldBe` 45.91010855311742 + score conf' "CCCTGCCCCGCCCCCGAGGCCGCCGGCGGCCCTAGCGTGTTCCTGTTCCCTCCTAAGCCTAAGGACACCCTGATGATCAGCAGAACCCCCGAGGTGACCTGCGTGGTGGTGGACGTGAGCCACGAGGACCCTGAGGTGAAGTTCAATTGGTACGTGGACGGCGTGGAGGTGCACAACGCCAAGACCAAGCCTAGAGAGGAGCAGTACAACAGCACCTACAGAGTGGTGAGCGTGCTGACCGTGCTGCACCAAGACTGGCTGAACGGCAAGGAGTACAAGTGCAAGGTGAGCAACAAGGCCCTGCCCGCCCCTATCGAGAAGACCATCAGCAAGGCCAAG" + `shouldBe` 68.00968914282373 + +scoreFunWithForbiddenSeq :: Spec +scoreFunWithForbiddenSeq = + describe "scoreFunWithForbiddenSeq" $ + it "should correct count score with forbidden sequence" $ do + score conf "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGCAGCACCATACTCCCCCGGC" `shouldBe` -3.4521030544284343 + score conf "GCCAGCGGCGACCGGCGGGTCAAGACCCACACCTGTCCT" `shouldBe` -17.346165749671798 + score conf "ACAGCCAGCGAATAAACCCCGAGGCCGCCGGCGGCCCTAGCGTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAATGACTACAGCCAAG" `shouldBe` 3.4278532842107197 + +scoreFunDifferentForbiddenSeqWeight :: Spec +scoreFunDifferentForbiddenSeqWeight = + describe "scoreFunDifferentForbiddenSeqWeight" $ + it "should correct count score with forbidden sequence weight = 4" $ do + let conf' = CodonConfig 3 1 (CodonScoreConfig 1 0.5 1.4 40 0.001 2.6 100 4 43) + score conf' "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGCAGCACCATACTCCCCCGGC" `shouldBe` -153.45210305442845 + score conf' "GCCAGCGGCGACCGGCGGGTCAAGACCCACACCTGTCCT" `shouldBe` -167.3461657496718 + score conf' "ACAGCCAGCGAATAAACCCCGAGGCCGCCGGCGGCCCTAGCGTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAATGACTACAGCCAAG" `shouldBe` -146.57214671578927 + score conf' "ACAGCCAGCGAATAAACCCCGACGATCGGGCCGCCGGCGAATGACAGCCAAG" `shouldBe` -137.875795016417 + +scoreFunDifferentGCDesired :: Spec +scoreFunDifferentGCDesired = + describe "scoreFunDifferentGCDesired" $ + it "should correct count score with gc desired = 60" $ do + let conf' = CodonConfig 3 1 (CodonScoreConfig 1 0.5 1.4 40 0.001 2.6 100 1 60) + score conf' "ATGGAGACCGACACCCTGCTGCTGTGGGTGCTGCTGCTGTGGGTGCCTGGCAGCACCGGC" `shouldBe` 80.83499151084192 + score conf' "GCCAGCGGCGACAAGACCCACACCTGTCCT" `shouldBe` 80.74362252733329 + score conf' "CCCTGCCCCGCCCCCGAGGCCGCCGGCGGCCCTAGCGTGTTCCTGTTCCCTCCTAAGCCTAAGGACACCCTGATGATCAGCAGAACCCCCGAGGTGACCTGCGTGGTGGTGGACGTGAGCCACGAGGACCCTGAGGTGAAGTTCAATTGGTACGTGGACGGCGTGGAGGTGCACAACGCCAAGACCAAGCCTAGAGAGGAGCAGTACAACAGCACCTACAGAGTGGTGAGCGTGCTGACCGTGCTGCACCAAGACTGGCTGAACGGCAAGGAGTACAAGTGCAAGGTGAGCAACAAGGCCCTGCCCGCCCCTATCGAGAAGACCATCAGCAAGGCCAAG" + `shouldBe` 96.19662511761598