diff --git a/Biobase/Clustal.hs b/Biobase/Clustal.hs new file mode 100644 --- /dev/null +++ b/Biobase/Clustal.hs @@ -0,0 +1,10 @@ +-- | Types and functions for Clustal +-- + +module Biobase.Clustal + ( module Biobase.Clustal.Types + , module Biobase.Clustal.Import + ) where + +import Biobase.Clustal.Import +import Biobase.Clustal.Types diff --git a/Biobase/Clustal/Import.hs b/Biobase/Clustal/Import.hs new file mode 100644 --- /dev/null +++ b/Biobase/Clustal/Import.hs @@ -0,0 +1,294 @@ +-- | Parse Clustal output +-- For more information on Clustal tools consult: +module Biobase.Clustal.Import ( + parseClustalAlignment, + readClustalAlignment, + parseStructuralClustalAlignment, + readStructuralClustalAlignment, + parseClustalSummary, + readClustalSummary, + module Biobase.Clustal.Types + ) where + +import Biobase.Clustal.Types +import Text.ParserCombinators.Parsec +import Control.Monad +import Data.List +import qualified Data.Text as T +readDouble :: String -> Double +readDouble = read + +readInt :: String -> Int +readInt = read + +-- | Parse the input as ClustalSummary datatype +genParserClustalSummary :: GenParser Char st ClustalSummary +genParserClustalSummary = do + newline + newline + newline + space + string "CLUSTAL " + version <- many1 (noneOf " ") + many1 (noneOf "\n") + newline + newline + newline + string "Sequence format is " + sequenceFormat' <- many1 (noneOf "\n") + newline + sequenceParametersList <- many1 (try genParserSequenceParameters) + string "Start of Pairwise alignments" + newline + string "Aligning..." + newline + newline + pairwiseAlignmentSummaryList <- many1 genParserPairwiseAlignmentSummary + string "Guide tree file created: [" + guideTreeFileName' <- many1 (noneOf "]") + char ']' + newline + newline + string "There are " + numberOfGroups <- many1 digit + string " groups" + newline + string "Start of Multiple Alignment" + newline + newline + string "Aligning..." + newline + groupSummaryList <- many1 genParserGroupSummary + string "Alignment Score " + alignmentScore' <- many1 digit + newline + newline + string "CLUSTAL-Alignment file created [" + alignmentFileName' <- many1 (noneOf "]") + char ']' + newline + newline + eof + return $ ClustalSummary (T.pack version) (T.pack sequenceFormat') sequenceParametersList pairwiseAlignmentSummaryList (T.pack guideTreeFileName') (readInt numberOfGroups) groupSummaryList (readInt alignmentScore') (T.pack alignmentFileName') + +genParserGroupSummary :: GenParser Char st GroupSummary +genParserGroupSummary = do + string "Group " + groupIndex <- many1 digit + string ":" + optional space + optional (string "Sequences:") + many1 space + sequenceNumber <- optionMaybe (many1 digit) + optional (many1 space) + string "Score:" <|> string "Delayed" + groupScore' <- optionMaybe (many1 digit) + newline + return $ GroupSummary (readInt groupIndex) (liftM readInt sequenceNumber) (liftM readInt groupScore') + +genParserPairwiseAlignmentSummary :: GenParser Char st PairwiseAlignmentSummary +genParserPairwiseAlignmentSummary = do + string "Sequences (" + firstSeqIndex <- many1 digit + string ":" + secondSeqIndex <- many1 digit + string ") Aligned. Score:" + many1 space + pairwiseScore <- many1 digit + newline + return $ PairwiseAlignmentSummary (readInt firstSeqIndex) (readInt secondSeqIndex) (readInt pairwiseScore) + +genParserSequenceParameters :: GenParser Char st SequenceParameters +genParserSequenceParameters = do + string "Sequence " + sequenceIndexParam <- many1 digit + string ": " + sequenceIdentifierParam <- many1 (noneOf " ") + spaces + sequenceLengthParam <- many1 digit + space + string "bp" + newline + return $ SequenceParameters (readInt sequenceIndexParam) (T.pack sequenceIdentifierParam) (readInt sequenceLengthParam) + +-- | Parse the input as ClustalAlignment datatype +genParserClustalAlignment :: GenParser Char st ClustalAlignment +genParserClustalAlignment = do + string "CLUSTAL" + many1 (noneOf "\n") + many1 (try newline) + alignmentSlices <- many1 (try genParserClustalAlignmentSlice) + optional newline + eof + return (mergealignmentSlices alignmentSlices) + +mergealignmentSlices :: [ClustalAlignmentSlice] -> ClustalAlignment +mergealignmentSlices slices = alignment + where entrySlicesList = concatMap entrySlices slices -- list of lists of entry slices + sequenceIdentifiers = nub (map entrySequenceSliceIdentifier entrySlicesList) + mergedAlignmentEntries = map (constructAlignmentEntries entrySlicesList) sequenceIdentifiers + mergedConservationTrack = concatMap conservationTrackSlice slices + alignment = ClustalAlignment mergedAlignmentEntries (T.pack mergedConservationTrack) + +constructAlignmentEntries :: [ClustalAlignmentEntrySlice] -> String -> ClustalAlignmentEntry +constructAlignmentEntries slices entryIdentifier= entry + where currentSlices = filter (\a -> entrySequenceSliceIdentifier a == entryIdentifier) slices + entrySequence = concatMap entryAlignedSliceSequence currentSlices + entry = ClustalAlignmentEntry (T.pack entryIdentifier) (T.pack entrySequence) + +genParserClustalAlignmentSlice :: GenParser Char st ClustalAlignmentSlice +genParserClustalAlignmentSlice = do + entrySlices' <- many1 genParserClustalEntrySlice + --extract length of identifier and spacer to determine offset of conservation track + let offsetLenght = length (entrySequenceSliceIdentifier (head entrySlices')) + spacerLength (head entrySlices') + conservationTrackSliceChoice <- choice [(lookAhead (string "\n")), (try (genParserConservationTrackSlice offsetLenght))] + --spacerAndConservationTrackSlice <- many1 (noneOf "\n") + --let conservationTrackSlice' = drop offsetLenght spacerAndConservationTrackSlice + --newline + let conservationTrackSlice' = if (conservationTrackSliceChoice == "\n") then "" else conservationTrackSliceChoice + optional newline + return $ ClustalAlignmentSlice entrySlices' conservationTrackSlice' + +genParserConservationTrackSlice :: Int -> GenParser Char st String +genParserConservationTrackSlice offsetLenght = do + spacerAndConservationTrackSlice <- many1 (noneOf "\n") + let conservationTrackSlice' = drop offsetLenght spacerAndConservationTrackSlice + newline + return $ conservationTrackSlice' + +genParserClustalEntrySlice :: GenParser Char st ClustalAlignmentEntrySlice +genParserClustalEntrySlice = do + sliceIdentifier <- many1 (noneOf " \n") + spacer <- many1 (char ' ') + sliceSequence <- parseNucleotideAlignmentEntry + newline + return $ ClustalAlignmentEntrySlice sliceIdentifier sliceSequence (length spacer) + +--Structural Clustal Parser functions + +-- | Parse the input as ClustalAlignment datatype as used in mlocarna +genParserStructuralClustalAlignment :: GenParser Char st StructuralClustalAlignment +genParserStructuralClustalAlignment = do + genParseMlocarnaHeader + alignmentSlices <- many1 (try genParserStructuralClustalAlignmentSlice) + secondaryStructure <- genSecondaryStructure + energy' <- genParseEnergy + eof + return (mergeStructuralAlignmentSlices (concat alignmentSlices) secondaryStructure energy') + +genSecondaryStructure :: GenParser Char st String +genSecondaryStructure = do + string "alifold" + secondaryStructure <- many1 (try genSecondaryStructureSlice) + return (concat secondaryStructure) + +genSecondaryStructureSlice :: GenParser Char st String +genSecondaryStructureSlice = do + many1 space + secondaryStructureSlice <- many1 (try (oneOf ".()")) + choice [try (string "\n"),try (string " ")] + return secondaryStructureSlice + +genParseEnergy :: GenParser Char st Double +genParseEnergy = do + string "(" + many space + energy' <- many1 (noneOf " ") + optional space + char ('=') + many1 (noneOf "\n") + newline + return (readDouble energy') + +genParseMlocarnaHeader :: GenParser Char st String +genParseMlocarnaHeader = do + string "mLocARNA" + many1 (choice [alphaNum,oneOf "-: ()."]) + newline + string "Copyright" + many1 (choice [alphaNum,char ' ']) + newline + newline + string "Compute pair probs ..." + newline + optional (try (string "Compute pairwise alignments ... ")) + optional (try newline) + string "Perform progressive alignment ..." + many1 newline + return "" + +mergeStructuralAlignmentSlices :: [StructuralClustalAlignmentEntrySlice] -> String -> Double -> StructuralClustalAlignment +mergeStructuralAlignmentSlices slices secondaryStructure energy' = alignment + where sequenceIdentifiers = nub (map structuralEntrySequenceSliceIdentifier slices) + mergedAlignmentEntries = map (constructStructuralAlignmentEntries slices) sequenceIdentifiers + alignment = StructuralClustalAlignment mergedAlignmentEntries (T.pack secondaryStructure) energy' + +constructStructuralAlignmentEntries :: [StructuralClustalAlignmentEntrySlice] -> String -> ClustalAlignmentEntry +constructStructuralAlignmentEntries slices entryIdentifier= entry + where currentSlices = filter (\a -> structuralEntrySequenceSliceIdentifier a == entryIdentifier) slices + entrySequence = concatMap structuralEntryAlignedSliceSequence currentSlices + entry = ClustalAlignmentEntry (T.pack entryIdentifier) (T.pack entrySequence) + +genParserStructuralClustalAlignmentSlice :: GenParser Char st [StructuralClustalAlignmentEntrySlice] +genParserStructuralClustalAlignmentSlice = do + entrySlices' <- many1 (try genParserStructuralClustalEntrySlice) + optional newline + return entrySlices' + +genParserStructuralClustalEntrySlice :: GenParser Char st StructuralClustalAlignmentEntrySlice +genParserStructuralClustalEntrySlice = do + sliceIdentifier <- many1 (noneOf " ") + many1 (char ' ') + sliceSequence <- parseNucleotideAlignmentEntry + newline + return $ StructuralClustalAlignmentEntrySlice (filter (/='\n') sliceIdentifier) sliceSequence + +-- exported functions + +-- | Parse Clustal alignment (.aln) from String +parseClustalAlignment :: String -> Either ParseError ClustalAlignment +parseClustalAlignment = parse genParserClustalAlignment "genParserClustalAlignment" + +-- | Parse Clustal alignment (.aln) from filehandle +readClustalAlignment :: String -> IO (Either ParseError ClustalAlignment) +readClustalAlignment = parseFromFile genParserClustalAlignment + +-- | Parse Clustal alignment (.aln) with secondary structure in dot-bracket notation from String (as produced by mlocarna) +parseStructuralClustalAlignment :: String -> Either ParseError StructuralClustalAlignment +parseStructuralClustalAlignment = parse genParserStructuralClustalAlignment "genParserStructuralClustalAlignment" + +-- | Parse Clustal alignment (.aln) with secondary structure in dot-bracket notation from filehandle (as produced by mlocarna) +readStructuralClustalAlignment :: String -> IO (Either ParseError StructuralClustalAlignment) +readStructuralClustalAlignment = parseFromFile genParserStructuralClustalAlignment + +-- | Parse Clustal summary (printed to STDOUT) from String +parseClustalSummary :: String -> Either ParseError ClustalSummary +parseClustalSummary = parse genParserClustalSummary "genParserClustalSummary" + +-- | Parse Clustal summary (printed to STDOUT) from file +readClustalSummary :: String -> IO (Either ParseError ClustalSummary) +readClustalSummary = parseFromFile genParserClustalSummary + +-- | Parse nucleotide sequence. Allowed letters according to IUPAC +--parseNucleotideSequence :: GenParser Char st String +--parseNucleotideSequence = do +-- nucleotideSequence <- many1 (oneOf "RYSWKMBDHVNATUGCryswkmbdhvnatugc") +-- return $ nucleotideSequence + +-- | Parse nucleotide alignment entry. Allowed letters according to IUPAC and commonly used gap characters +parseNucleotideAlignmentEntry :: GenParser Char st String +parseNucleotideAlignmentEntry = do + entry <- many1 (oneOf "~_-.RYSWKMBDHVNATUGCryswkmbdhvnatugc") + return $ entry + +-- | Parse protein amino acid code sequence. Allowed letters according to IUPAC +--parseProteinSequence :: GenParser Char st String +--parseProteinSequence = do +-- proteinSequence <- many1 (oneOf "ABCDEFGHIKLMNPQRSTVWXYZabcdefghiklmnpqrstvwxyz") +-- return $ proteinSequence + +-- | Parse protein amino acid code alignment entry. Allowed letters according to IUPAC and commonly used gap characters +--parseProteinAlignmentEntry :: GenParser Char st String +--parseProteinAlignmentEntry = do +-- entry <- many1 (oneOf "~_-.ABCDEFGHIKLMNPQRSTVWXYZabcdefghiklmnpqrstvwxyz") +-- return $ entry diff --git a/Biobase/Clustal/Types.hs b/Biobase/Clustal/Types.hs new file mode 100644 --- /dev/null +++ b/Biobase/Clustal/Types.hs @@ -0,0 +1,131 @@ +-- | This module contains data structures for the Clustal tools +-- For more information on Clustal tools consult: + +module Biobase.Clustal.Types where +import qualified Data.Vector as V +import qualified Data.Text as T + +-- | Data type for clustal summary, containing information about the alignment process, usually printed to STDOUT +data ClustalSummary = ClustalSummary + { + clustalw2version :: T.Text, + sequenceFormat :: T.Text, + parametersOfInputSequences :: [SequenceParameters], + pairwiseAlignmentSummaries :: [PairwiseAlignmentSummary], + guideTreeFileName :: T.Text, + groupNumber :: Int, + groupSummaries :: [GroupSummary], + alignmentScore :: Int, + alignmentFileName :: T.Text + } + deriving (Show, Eq) + +data SequenceParameters = SequenceParameters + { + inputSequenceIndex :: Int, + inputSequenceIdentifier :: T.Text, + inputSequenceLength :: Int + } + deriving (Show, Eq) + +data PairwiseAlignmentSummary = PairwiseAlignmentSummary + { + firstSequenceIndex :: Int, + secondSequenceIndex :: Int, + pairwiseAlignmentScore :: Int + } + deriving (Show, Eq) + +data GroupSummary = GroupSummary + { + alignmentGroupIndex :: Int, + numberOfAlignedSequences :: Maybe Int, + groupScore :: Maybe Int + } + deriving (Show, Eq) + +-- | Data structure for Clustal alignment format +data ClustalAlignment = ClustalAlignment + { + alignmentEntries :: [ClustalAlignmentEntry], + conservationTrack :: T.Text + } + deriving (Eq) + +instance Show ClustalAlignment where + show (ClustalAlignment _alignmentEntries _conservationTrack ) + | not (null _alignmentEntries) = header ++ alignmentString + | otherwise = header + where header = "CLUSTAL W (1.8) multiple sequence alignment\n\n\n" + longestSequenceIdLength = (maximum (map T.length (map entrySequenceIdentifier _alignmentEntries))) + 1 + totalSequenceLength = T.length (entryAlignedSequence (head _alignmentEntries)) + alignmentString = showAlignment totalSequenceLength longestSequenceIdLength 0 _alignmentEntries _conservationTrack + +showAlignment :: Int -> Int -> Int -> [ClustalAlignmentEntry] -> T.Text -> String +showAlignment totalSequenceLength longestSequenceIdLength currentWindowPosition _alignmentEntries _conservationTrack + | totalSequenceLength == 0 = [] + | currentWindowPosition < totalSequenceLength = showAlignmentBlock longestSequenceIdLength currentWindowPosition _alignmentEntries _conservationTrack ++ (showAlignment totalSequenceLength longestSequenceIdLength (currentWindowPosition + 60) _alignmentEntries _conservationTrack) + | currentWindowPosition == totalSequenceLength = [] + | otherwise = "" + +showAlignmentBlock :: Int -> Int -> [ClustalAlignmentEntry] -> T.Text -> String +showAlignmentBlock longestSequenceIdLength currentWindowPosition _alignmentEntries _conservationTrack = blockString + where blockString = entries ++ extraTrack ++ "\n" + entries = concatMap (showAlignmentLine longestSequenceIdLength currentWindowPosition) _alignmentEntries + extraTrack = concat (replicate longestSequenceIdLength " ") ++ V.toList (V.slice currentWindowPosition 60 (V.fromList (T.unpack _conservationTrack))) ++ "\n" + +showAlignmentLine :: Int -> Int -> ClustalAlignmentEntry -> String +showAlignmentLine longestSequenceIdLength currentWindowPosition _alignmentEntry = T.unpack (entrySequenceIdentifier _alignmentEntry) ++ concat (replicate (longestSequenceIdLength - T.length (entrySequenceIdentifier _alignmentEntry)) " ") ++ V.toList (V.slice currentWindowPosition 60 (V.fromList (T.unpack (entryAlignedSequence _alignmentEntry)))) ++ "\n" + +data ClustalAlignmentEntry = ClustalAlignmentEntry + { + entrySequenceIdentifier :: T.Text, + entryAlignedSequence :: T.Text + } + deriving (Show, Eq) + +data ClustalAlignmentSlice = ClustalAlignmentSlice + { + entrySlices :: [ClustalAlignmentEntrySlice], + conservationTrackSlice :: String + } + deriving (Show, Eq) + +data ClustalAlignmentEntrySlice = ClustalAlignmentEntrySlice + { + entrySequenceSliceIdentifier :: String, + entryAlignedSliceSequence :: String, + spacerLength :: Int + } + deriving (Show, Eq) + +-- | Data structure for structural Clustal alignment format +data StructuralClustalAlignment = StructuralClustalAlignment + { + structuralAlignmentEntries :: [ClustalAlignmentEntry], + secondaryStructureTrack :: T.Text, + energy :: Double + } + deriving (Eq) + +instance Show StructuralClustalAlignment where + show (StructuralClustalAlignment _alignmentEntries _secondaryStructureTrack _energy) + | not (null _alignmentEntries) = header ++ alignmentString + | otherwise = header + where header = "CLUSTAL W \n\n" + longestSequenceIdLength = (maximum (map T.length (map entrySequenceIdentifier _alignmentEntries))) + 1 + totalSequenceLength = T.length (entryAlignedSequence (head _alignmentEntries)) + alignmentString = showAlignment totalSequenceLength longestSequenceIdLength 0 _alignmentEntries _secondaryStructureTrack + +-- data StructuralClustalAlignmentSlice = StructuralClustalAlignmentSlice +-- { +-- structuralEntrySlices :: [StructuralClustalAlignmentEntrySlice] +-- } +-- deriving (Show, Eq) + +data StructuralClustalAlignmentEntrySlice = StructuralClustalAlignmentEntrySlice + { + structuralEntrySequenceSliceIdentifier :: String, + structuralEntryAlignedSliceSequence :: String + } + deriving (Show, Eq) diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,57 @@ +-*-change-log-*- + +### 1.3.0 [Florian Eggenhofer](mailto:egg@cs.uni-freiburg.de) 14. November 2019 + + * Fixed requested tick number for compilation with GHC 8.6.* + * Changed to Biobase style + +### 1.2.3 [Florian Eggenhofer](mailto:egg@cs.uni-freiburg.de) 12. March 2018 + + * Fixed parsing of additional newline in Biopythons AlignIO output without conservation track + +### 1.2.2 [Florian Eggenhofer](mailto:egg@cs.uni-freiburg.de) 07. March 2018 + + * Clustal parser can now parse alignments with missing consensus annotation + +### 1.2.1 [Florian Eggenhofer](mailto:egg@cs.uni-freiburg.de) 06. February 2017 + + * Structural alignment parser now works with multiline consensus structures + +### 1.2.0 [Florian Eggenhofer](mailto:egg@cs.uni-freiburg.de) 07. January 2017 + + * Changed datastructures for sequence identifers and sequences to Data.Text + +### 1.1.4 [Florian Eggenhofer](mailto:egg@cs.uni-freiburg.de) 30. May 2016 + + * Fixed a bug in output of clustal alignments with sequence length of 60 + +### 1.1.3 [Florian Eggenhofer](mailto:florian.eggenhofer@univie.ac.at) 4. July 2015 + + * Nucleotide sequences are now parsed by a unified function in line with IUPAC nucleotide code + +### 1.1.2 [Florian Eggenhofer](mailto:florian.eggenhofer@univie.ac.at) 3. July 2015 + + * Included parsing of optional field in mlocarna clustal output + +### 1.1.1 [Florian Eggenhofer](mailto:florian.eggenhofer@univie.ac.at) 2. July 2015 + + * Added support for cmalign clustal output . + +### 1.1.0 [Florian Eggenhofer](mailto:florian.eggenhofer@univie.ac.at) 1. July 2015 + + * Added Hspec test-suite for parsing functions + * Added Show instances for ClustalAlignment and StructuralClustalAlignment + +### 1.0.3 [Florian Eggenhofer](mailto:florian.eggenhofer@univie.ac.at) 19. April 2015 + + * Added Y (pyrimidine) and R (purine) to sequence characters + +### 1.0.2 [Florian Eggenhofer](mailto:florian.eggenhofer@univie.ac.at> 19. March 2015 + + * Linebreaks are now filtered from structural alignment sequence identifiers + +### 1.0.1 [Florian Eggenhofer](mailto:florian.eggenhofer@univie.ac.at> 27. October 2014 + + * Fixed compiler warnings and updated documentation to mention structural clustal format + * Added -Wall and -O2 compiler options + * Added support for clustal alignments with secondary structure annotation diff --git a/ClustalParser.cabal b/ClustalParser.cabal --- a/ClustalParser.cabal +++ b/ClustalParser.cabal @@ -1,29 +1,28 @@ name: ClustalParser -version: 1.2.3 +version: 1.3.0 synopsis: Libary for parsing Clustal tools output description: Currently contains parsers and datatypes for: clustalw2, clustalo, mlocarna, cmalign . Clustal tools are multiple sequence alignment tools for biological sequence like DNA, RNA and Protein. For more information on clustal Tools refer to . . - Mlocarna is a multiple sequence alignment tool for RNA sequences with secondary structure output. + Mlocarna is a multiple sequence alignment tool for RNA sequences with secondary structure output. For more information on mlocarna refer to . . Cmalign is a multiple sequence alignment program based on RNA family models and produces ,among others, clustal output. It is part of infernal . - + license: GPL-3 license-file: LICENSE author: Florian Eggenhofer maintainer: egg@informatik.uni-freiburg.de --- copyright: category: Bioinformatics build-type: Simple -cabal-version: >=1.8 -Tested-With: GHC == 7.10.3, GHC == 8.0.2, GHC == 8.2.2 +cabal-version: >=1.10.0 +Tested-With: GHC == 8.4.4, GHC == 8.6.5, GHC == 8.8.1 extra-source-files: - README.md changelog + README.md ChangeLog.md source-repository head type: git @@ -31,31 +30,35 @@ source-repository this type: git - location: https://github.com/eggzilla/ClustalParser/tree/1.2.3 - tag: 1.2.3 + location: https://github.com/eggzilla/ClustalParser/tree/1.3.0 + tag: 1.3.0 -library - -- Modules exported by the library. - exposed-modules: Bio.ClustalParser - other-modules: Bio.ClustalData +Library + Hs-Source-Dirs: . + ghc-options: -Wall -fno-warn-unused-do-bind -fsimpl-tick-factor=500 + default-language: Haskell2010 + build-depends: base >=4.5 && <5, + parsec>=3.1.9, + vector, + text + Exposed-Modules: Biobase.Clustal + Biobase.Clustal.Import + Biobase.Clustal.Types - -- Other library packages from which modules are imported. - build-depends: base >=4.5 && <5, parsec>=3.1.9, vector, text - -- compiler options - ghc-options: -Wall -fno-warn-unused-do-bind - -- Directories containing source files. - hs-source-dirs: src - executable ClustalParserTest main-is: ClustalParserTest.hs - build-depends: base >= 4 && <= 5, cmdargs, ClustalParser, either-unwrap - -- compiler options - ghc-options: -Wall + build-depends: base >= 4 && <= 5, cmdargs, ClustalParser, either-unwrap, text, vector, parsec + -- compiler options + ghc-options: -Wall -fno-warn-unused-do-bind -fsimpl-tick-factor=500 + default-language: Haskell2010 + other-modules: Biobase.Clustal.Import + Biobase.Clustal.Types test-suite hspec build-depends: base, parsec, ClustalParser, hspec >= 2.0, text, hspec-discover hs-source-dirs: test-suite + default-language: Haskell2010 main-is: Spec.hs type: exitcode-stdio-1.0 other-modules: ClustalParserSpec diff --git a/ClustalParserTest.hs b/ClustalParserTest.hs --- a/ClustalParserTest.hs +++ b/ClustalParserTest.hs @@ -1,16 +1,15 @@ --runghc -package-db --ghc-arg=.cabal-sandbox/x86_64-linux-ghc-8.0.1-packages.conf.d ClustalParserTest.hs struct-multiline.mlocarna module Main where - + import System.Environment (getArgs) -import Bio.ClustalParser +import Biobase.Clustal.Import --import Data.Either.Unwrap main :: IO () main = do args <- getArgs - let input_file = (head args) + let input_file = (head args) -- read Clustal outputfile parsedinput <- readStructuralClustalAlignment input_file print parsedinput --print $ structuralAlignmentEntries (fromRight parsedinput) - diff --git a/changelog b/changelog deleted file mode 100644 --- a/changelog +++ /dev/null @@ -1,30 +0,0 @@ --*-change-log-*- -1.2.3 Florian Eggenhofer 12. March 2018 - * Fixed parsing of additional newline in Biopythons AlignIO output without conservation track -1.2.2 Florian Eggenhofer 07. March 2018 - * Clustal parser can now parse alignments with missing consensus - annotation -1.2.1 Florian Eggenhofer 06. February 2017 - * Structural alignment parser now works with multiline consensus structures -1.2.0 Florian Eggenhofer 07. January 2017 - * Changed datastructures for sequence identifers and sequences to Data.Text -1.1.4 Florian Eggenhofer 30. May 2016 - * Fixed a bug in output of clustal alignments with sequence length of 60 -1.1.3 Florian Eggenhofer 4. July 2015 - * Nucleotide sequences are now parsed by a unified function in line - with IUPAC nucleotide code -1.1.2 Florian Eggenhofer 3. July 2015 - * Included parsing of optional field in mlocarna clustal output -1.1.1 Florian Eggenhofer 2. July 2015 - * Added support for cmalign clustal output . -1.1.0 Florian Eggenhofer 1. July 2015 - * Added Hspec test-suite for parsing functions - * Added Show instances for ClustalAlignment and StructuralClustalAlignment -1.0.3 Florian Eggenhofer 19. April 2015 - * Added Y (pyrimidine) and R (purine) to sequence characters -1.0.2 Florian Eggenhofer 19. March 2015 - * Linebreaks are now filtered from structural alignment sequence identifiers -1.0.1 Florian Eggenhofer 27. October 2014 - * Fixed compiler warnings and updated documentation to mention structural clustal format - * Added -Wall and -O2 compiler options - * Added support for clustal alignments with secondary structure annotation diff --git a/src/Bio/ClustalData.hs b/src/Bio/ClustalData.hs deleted file mode 100644 --- a/src/Bio/ClustalData.hs +++ /dev/null @@ -1,133 +0,0 @@ --- | This module contains data structures for the Clustal tools --- For more information on Clustal tools consult: - -module Bio.ClustalData where -import qualified Data.Vector as V -import qualified Data.Text as T - --- | Data type for clustal summary, containing information about the alignment process, usually printed to STDOUT -data ClustalSummary = ClustalSummary - { - clustalw2version :: T.Text, - sequenceFormat :: T.Text, - parametersOfInputSequences :: [SequenceParameters], - pairwiseAlignmentSummaries :: [PairwiseAlignmentSummary], - guideTreeFileName :: T.Text, - groupNumber :: Int, - groupSummaries :: [GroupSummary], - alignmentScore :: Int, - alignmentFileName :: T.Text - } - deriving (Show, Eq) - -data SequenceParameters = SequenceParameters - { - inputSequenceIndex :: Int, - inputSequenceIdentifier :: T.Text, - inputSequenceLength :: Int - } - deriving (Show, Eq) - -data PairwiseAlignmentSummary = PairwiseAlignmentSummary - { - firstSequenceIndex :: Int, - secondSequenceIndex :: Int, - pairwiseAlignmentScore :: Int - } - deriving (Show, Eq) - -data GroupSummary = GroupSummary - { - alignmentGroupIndex :: Int, - numberOfAlignedSequences :: Maybe Int, - groupScore :: Maybe Int - } - deriving (Show, Eq) - --- | Data structure for Clustal alignment format -data ClustalAlignment = ClustalAlignment - { - alignmentEntries :: [ClustalAlignmentEntry], - conservationTrack :: T.Text - } - deriving (Eq) - -instance Show ClustalAlignment where - show (ClustalAlignment _alignmentEntries _conservationTrack ) - | not (null _alignmentEntries) = header ++ alignmentString - | otherwise = header - where header = "CLUSTAL W (1.8) multiple sequence alignment\n\n\n" - longestSequenceIdLength = (maximum (map T.length (map entrySequenceIdentifier _alignmentEntries))) + 1 - totalSequenceLength = T.length (entryAlignedSequence (head _alignmentEntries)) - alignmentString = showAlignment totalSequenceLength longestSequenceIdLength 0 _alignmentEntries _conservationTrack - -showAlignment :: Int -> Int -> Int -> [ClustalAlignmentEntry] -> T.Text -> String -showAlignment totalSequenceLength longestSequenceIdLength currentWindowPosition _alignmentEntries _conservationTrack - | totalSequenceLength == 0 = [] - | currentWindowPosition < totalSequenceLength = showAlignmentBlock longestSequenceIdLength currentWindowPosition _alignmentEntries _conservationTrack ++ (showAlignment totalSequenceLength longestSequenceIdLength (currentWindowPosition + 60) _alignmentEntries _conservationTrack) - | currentWindowPosition == totalSequenceLength = [] - | otherwise = "" - -showAlignmentBlock :: Int -> Int -> [ClustalAlignmentEntry] -> T.Text -> String -showAlignmentBlock longestSequenceIdLength currentWindowPosition _alignmentEntries _conservationTrack = blockString - where blockString = entries ++ extraTrack ++ "\n" - entries = concatMap (showAlignmentLine longestSequenceIdLength currentWindowPosition) _alignmentEntries - extraTrack = concat (replicate longestSequenceIdLength " ") ++ V.toList (V.slice currentWindowPosition 60 (V.fromList (T.unpack _conservationTrack))) ++ "\n" - -showAlignmentLine :: Int -> Int -> ClustalAlignmentEntry -> String -showAlignmentLine longestSequenceIdLength currentWindowPosition _alignmentEntry = T.unpack (entrySequenceIdentifier _alignmentEntry) ++ concat (replicate (longestSequenceIdLength - T.length (entrySequenceIdentifier _alignmentEntry)) " ") ++ V.toList (V.slice currentWindowPosition 60 (V.fromList (T.unpack (entryAlignedSequence _alignmentEntry)))) ++ "\n" - -data ClustalAlignmentEntry = ClustalAlignmentEntry - { - entrySequenceIdentifier :: T.Text, - entryAlignedSequence :: T.Text - } - deriving (Show, Eq) - -data ClustalAlignmentSlice = ClustalAlignmentSlice - { - entrySlices :: [ClustalAlignmentEntrySlice], - conservationTrackSlice :: String - } - deriving (Show, Eq) - -data ClustalAlignmentEntrySlice = ClustalAlignmentEntrySlice - { - entrySequenceSliceIdentifier :: String, - entryAlignedSliceSequence :: String, - spacerLength :: Int - } - deriving (Show, Eq) - --- | Data structure for structural Clustal alignment format -data StructuralClustalAlignment = StructuralClustalAlignment - { - structuralAlignmentEntries :: [ClustalAlignmentEntry], - secondaryStructureTrack :: T.Text, - energy :: Double - } - deriving (Eq) - -instance Show StructuralClustalAlignment where - show (StructuralClustalAlignment _alignmentEntries _secondaryStructureTrack _energy) - | not (null _alignmentEntries) = header ++ alignmentString - | otherwise = header - where header = "CLUSTAL W \n\n" - longestSequenceIdLength = (maximum (map T.length (map entrySequenceIdentifier _alignmentEntries))) + 1 - totalSequenceLength = T.length (entryAlignedSequence (head _alignmentEntries)) - alignmentString = showAlignment totalSequenceLength longestSequenceIdLength 0 _alignmentEntries _secondaryStructureTrack - --- data StructuralClustalAlignmentSlice = StructuralClustalAlignmentSlice --- { --- structuralEntrySlices :: [StructuralClustalAlignmentEntrySlice] --- } --- deriving (Show, Eq) - -data StructuralClustalAlignmentEntrySlice = StructuralClustalAlignmentEntrySlice - { - structuralEntrySequenceSliceIdentifier :: String, - structuralEntryAlignedSliceSequence :: String - } - deriving (Show, Eq) - - diff --git a/src/Bio/ClustalParser.hs b/src/Bio/ClustalParser.hs deleted file mode 100644 --- a/src/Bio/ClustalParser.hs +++ /dev/null @@ -1,294 +0,0 @@ --- | Parse Clustal output --- For more information on Clustal tools consult: -module Bio.ClustalParser ( - parseClustalAlignment, - readClustalAlignment, - parseStructuralClustalAlignment, - readStructuralClustalAlignment, - parseClustalSummary, - readClustalSummary, - module Bio.ClustalData - ) where - -import Bio.ClustalData -import Text.ParserCombinators.Parsec -import Control.Monad -import Data.List -import qualified Data.Text as T -readDouble :: String -> Double -readDouble = read - -readInt :: String -> Int -readInt = read - --- | Parse the input as ClustalSummary datatype -genParserClustalSummary :: GenParser Char st ClustalSummary -genParserClustalSummary = do - newline - newline - newline - space - string "CLUSTAL " - version <- many1 (noneOf " ") - many1 (noneOf "\n") - newline - newline - newline - string "Sequence format is " - sequenceFormat' <- many1 (noneOf "\n") - newline - sequenceParametersList <- many1 (try genParserSequenceParameters) - string "Start of Pairwise alignments" - newline - string "Aligning..." - newline - newline - pairwiseAlignmentSummaryList <- many1 genParserPairwiseAlignmentSummary - string "Guide tree file created: [" - guideTreeFileName' <- many1 (noneOf "]") - char ']' - newline - newline - string "There are " - numberOfGroups <- many1 digit - string " groups" - newline - string "Start of Multiple Alignment" - newline - newline - string "Aligning..." - newline - groupSummaryList <- many1 genParserGroupSummary - string "Alignment Score " - alignmentScore' <- many1 digit - newline - newline - string "CLUSTAL-Alignment file created [" - alignmentFileName' <- many1 (noneOf "]") - char ']' - newline - newline - eof - return $ ClustalSummary (T.pack version) (T.pack sequenceFormat') sequenceParametersList pairwiseAlignmentSummaryList (T.pack guideTreeFileName') (readInt numberOfGroups) groupSummaryList (readInt alignmentScore') (T.pack alignmentFileName') - -genParserGroupSummary :: GenParser Char st GroupSummary -genParserGroupSummary = do - string "Group " - groupIndex <- many1 digit - string ":" - optional space - optional (string "Sequences:") - many1 space - sequenceNumber <- optionMaybe (many1 digit) - optional (many1 space) - string "Score:" <|> string "Delayed" - groupScore' <- optionMaybe (many1 digit) - newline - return $ GroupSummary (readInt groupIndex) (liftM readInt sequenceNumber) (liftM readInt groupScore') - -genParserPairwiseAlignmentSummary :: GenParser Char st PairwiseAlignmentSummary -genParserPairwiseAlignmentSummary = do - string "Sequences (" - firstSeqIndex <- many1 digit - string ":" - secondSeqIndex <- many1 digit - string ") Aligned. Score:" - many1 space - pairwiseScore <- many1 digit - newline - return $ PairwiseAlignmentSummary (readInt firstSeqIndex) (readInt secondSeqIndex) (readInt pairwiseScore) - -genParserSequenceParameters :: GenParser Char st SequenceParameters -genParserSequenceParameters = do - string "Sequence " - sequenceIndexParam <- many1 digit - string ": " - sequenceIdentifierParam <- many1 (noneOf " ") - spaces - sequenceLengthParam <- many1 digit - space - string "bp" - newline - return $ SequenceParameters (readInt sequenceIndexParam) (T.pack sequenceIdentifierParam) (readInt sequenceLengthParam) - --- | Parse the input as ClustalAlignment datatype -genParserClustalAlignment :: GenParser Char st ClustalAlignment -genParserClustalAlignment = do - string "CLUSTAL" - many1 (noneOf "\n") - many1 (try newline) - alignmentSlices <- many1 (try genParserClustalAlignmentSlice) - optional newline - eof - return (mergealignmentSlices alignmentSlices) - -mergealignmentSlices :: [ClustalAlignmentSlice] -> ClustalAlignment -mergealignmentSlices slices = alignment - where entrySlicesList = concatMap entrySlices slices -- list of lists of entry slices - sequenceIdentifiers = nub (map entrySequenceSliceIdentifier entrySlicesList) - mergedAlignmentEntries = map (constructAlignmentEntries entrySlicesList) sequenceIdentifiers - mergedConservationTrack = concatMap conservationTrackSlice slices - alignment = ClustalAlignment mergedAlignmentEntries (T.pack mergedConservationTrack) - -constructAlignmentEntries :: [ClustalAlignmentEntrySlice] -> String -> ClustalAlignmentEntry -constructAlignmentEntries slices entryIdentifier= entry - where currentSlices = filter (\a -> entrySequenceSliceIdentifier a == entryIdentifier) slices - entrySequence = concatMap entryAlignedSliceSequence currentSlices - entry = ClustalAlignmentEntry (T.pack entryIdentifier) (T.pack entrySequence) - -genParserClustalAlignmentSlice :: GenParser Char st ClustalAlignmentSlice -genParserClustalAlignmentSlice = do - entrySlices' <- many1 genParserClustalEntrySlice - --extract length of identifier and spacer to determine offset of conservation track - let offsetLenght = length (entrySequenceSliceIdentifier (head entrySlices')) + spacerLength (head entrySlices') - conservationTrackSliceChoice <- choice [(lookAhead (string "\n")), (try (genParserConservationTrackSlice offsetLenght))] - --spacerAndConservationTrackSlice <- many1 (noneOf "\n") - --let conservationTrackSlice' = drop offsetLenght spacerAndConservationTrackSlice - --newline - let conservationTrackSlice' = if (conservationTrackSliceChoice == "\n") then "" else conservationTrackSliceChoice - optional newline - return $ ClustalAlignmentSlice entrySlices' conservationTrackSlice' - -genParserConservationTrackSlice :: Int -> GenParser Char st String -genParserConservationTrackSlice offsetLenght = do - spacerAndConservationTrackSlice <- many1 (noneOf "\n") - let conservationTrackSlice' = drop offsetLenght spacerAndConservationTrackSlice - newline - return $ conservationTrackSlice' - -genParserClustalEntrySlice :: GenParser Char st ClustalAlignmentEntrySlice -genParserClustalEntrySlice = do - sliceIdentifier <- many1 (noneOf " \n") - spacer <- many1 (char ' ') - sliceSequence <- parseNucleotideAlignmentEntry - newline - return $ ClustalAlignmentEntrySlice sliceIdentifier sliceSequence (length spacer) - ---Structural Clustal Parser functions - --- | Parse the input as ClustalAlignment datatype as used in mlocarna -genParserStructuralClustalAlignment :: GenParser Char st StructuralClustalAlignment -genParserStructuralClustalAlignment = do - genParseMlocarnaHeader - alignmentSlices <- many1 (try genParserStructuralClustalAlignmentSlice) - secondaryStructure <- genSecondaryStructure - energy' <- genParseEnergy - eof - return (mergeStructuralAlignmentSlices (concat alignmentSlices) secondaryStructure energy') - -genSecondaryStructure :: GenParser Char st String -genSecondaryStructure = do - string "alifold" - secondaryStructure <- many1 (try genSecondaryStructureSlice) - return (concat secondaryStructure) - -genSecondaryStructureSlice :: GenParser Char st String -genSecondaryStructureSlice = do - many1 space - secondaryStructureSlice <- many1 (try (oneOf ".()")) - choice [try (string "\n"),try (string " ")] - return secondaryStructureSlice - -genParseEnergy :: GenParser Char st Double -genParseEnergy = do - string "(" - many space - energy' <- many1 (noneOf " ") - optional space - char ('=') - many1 (noneOf "\n") - newline - return (readDouble energy') - -genParseMlocarnaHeader :: GenParser Char st String -genParseMlocarnaHeader = do - string "mLocARNA" - many1 (choice [alphaNum,oneOf "-: ()."]) - newline - string "Copyright" - many1 (choice [alphaNum,char ' ']) - newline - newline - string "Compute pair probs ..." - newline - optional (try (string "Compute pairwise alignments ... ")) - optional (try newline) - string "Perform progressive alignment ..." - many1 newline - return "" - -mergeStructuralAlignmentSlices :: [StructuralClustalAlignmentEntrySlice] -> String -> Double -> StructuralClustalAlignment -mergeStructuralAlignmentSlices slices secondaryStructure energy' = alignment - where sequenceIdentifiers = nub (map structuralEntrySequenceSliceIdentifier slices) - mergedAlignmentEntries = map (constructStructuralAlignmentEntries slices) sequenceIdentifiers - alignment = StructuralClustalAlignment mergedAlignmentEntries (T.pack secondaryStructure) energy' - -constructStructuralAlignmentEntries :: [StructuralClustalAlignmentEntrySlice] -> String -> ClustalAlignmentEntry -constructStructuralAlignmentEntries slices entryIdentifier= entry - where currentSlices = filter (\a -> structuralEntrySequenceSliceIdentifier a == entryIdentifier) slices - entrySequence = concatMap structuralEntryAlignedSliceSequence currentSlices - entry = ClustalAlignmentEntry (T.pack entryIdentifier) (T.pack entrySequence) - -genParserStructuralClustalAlignmentSlice :: GenParser Char st [StructuralClustalAlignmentEntrySlice] -genParserStructuralClustalAlignmentSlice = do - entrySlices' <- many1 (try genParserStructuralClustalEntrySlice) - optional newline - return entrySlices' - -genParserStructuralClustalEntrySlice :: GenParser Char st StructuralClustalAlignmentEntrySlice -genParserStructuralClustalEntrySlice = do - sliceIdentifier <- many1 (noneOf " ") - many1 (char ' ') - sliceSequence <- parseNucleotideAlignmentEntry - newline - return $ StructuralClustalAlignmentEntrySlice (filter (/='\n') sliceIdentifier) sliceSequence - --- exported functions - --- | Parse Clustal alignment (.aln) from String -parseClustalAlignment :: String -> Either ParseError ClustalAlignment -parseClustalAlignment = parse genParserClustalAlignment "genParserClustalAlignment" - --- | Parse Clustal alignment (.aln) from filehandle -readClustalAlignment :: String -> IO (Either ParseError ClustalAlignment) -readClustalAlignment = parseFromFile genParserClustalAlignment - --- | Parse Clustal alignment (.aln) with secondary structure in dot-bracket notation from String (as produced by mlocarna) -parseStructuralClustalAlignment :: String -> Either ParseError StructuralClustalAlignment -parseStructuralClustalAlignment = parse genParserStructuralClustalAlignment "genParserStructuralClustalAlignment" - --- | Parse Clustal alignment (.aln) with secondary structure in dot-bracket notation from filehandle (as produced by mlocarna) -readStructuralClustalAlignment :: String -> IO (Either ParseError StructuralClustalAlignment) -readStructuralClustalAlignment = parseFromFile genParserStructuralClustalAlignment - --- | Parse Clustal summary (printed to STDOUT) from String -parseClustalSummary :: String -> Either ParseError ClustalSummary -parseClustalSummary = parse genParserClustalSummary "genParserClustalSummary" - --- | Parse Clustal summary (printed to STDOUT) from file -readClustalSummary :: String -> IO (Either ParseError ClustalSummary) -readClustalSummary = parseFromFile genParserClustalSummary - --- | Parse nucleotide sequence. Allowed letters according to IUPAC ---parseNucleotideSequence :: GenParser Char st String ---parseNucleotideSequence = do --- nucleotideSequence <- many1 (oneOf "RYSWKMBDHVNATUGCryswkmbdhvnatugc") --- return $ nucleotideSequence - --- | Parse nucleotide alignment entry. Allowed letters according to IUPAC and commonly used gap characters -parseNucleotideAlignmentEntry :: GenParser Char st String -parseNucleotideAlignmentEntry = do - entry <- many1 (oneOf "~_-.RYSWKMBDHVNATUGCryswkmbdhvnatugc") - return $ entry - --- | Parse protein amino acid code sequence. Allowed letters according to IUPAC ---parseProteinSequence :: GenParser Char st String ---parseProteinSequence = do --- proteinSequence <- many1 (oneOf "ABCDEFGHIKLMNPQRSTVWXYZabcdefghiklmnpqrstvwxyz") --- return $ proteinSequence - --- | Parse protein amino acid code alignment entry. Allowed letters according to IUPAC and commonly used gap characters ---parseProteinAlignmentEntry :: GenParser Char st String ---parseProteinAlignmentEntry = do --- entry <- many1 (oneOf "~_-.ABCDEFGHIKLMNPQRSTVWXYZabcdefghiklmnpqrstvwxyz") --- return $ entry diff --git a/test-suite/ClustalParserSpec.hs b/test-suite/ClustalParserSpec.hs --- a/test-suite/ClustalParserSpec.hs +++ b/test-suite/ClustalParserSpec.hs @@ -1,6 +1,6 @@ module ClustalParserSpec (spec) where -import Bio.ClustalParser +import Biobase.Clustal.Import import Test.Hspec import Test.Hspec.QuickCheck import Text.Parsec.Error @@ -31,8 +31,8 @@ context "Parsing invalid input" $ do it "Returns ParseError" $ do (parseStructuralClustalAlignment "invalid input") `shouldBe` Left structuralClustalParseError - - + + clustalParseError :: ParseError clustalParseError = (addErrorMessage (Expect "\"CLUSTAL\"") (newErrorMessage clustalErrorMessage clustalErrorSourcePosition)) @@ -59,7 +59,7 @@ structuralClustalExample :: String structuralClustalExample = "mLocARNA --- multiple Local (and global) Alignment of RNA --- LocARNA 1.8.0\nCopyright Sebastian Will\n\nCompute pair probs ...\nPerform progressive alignment ...\n\n\n\nAB001721.1/2735-2851 CCCGGUGACUAUAGAGAGAGGGCCACACCCGUUCCCAUCCCGAACACGGAAGUUAAGCCUCUCAUCGCUGAUGGUACUAUGUGGUUCGCUGCAUGGGAGAGUAGGACGUUGCCGGGU\ngi|451991584:1-117 CCCGGUGACUAUAGAGAGAGGGCCACACCCGUUCCCAUCCCGAACACGGAAGUUAAGCCUCUCAUCGCUGAUGGUACUAUGUGGUUCGCUGCAUGGGAGAGUAGGACGUUGCCGGGU\nalifold (((((((((....(.(((.(((.....))).))).)...(((....)))..(((..(((((((((.((.(.(((.(((....))))))).))))))))).))..)))))))))))). (-38.30 = -38.30 + 0.00)\n" - + structuralClustalResult :: StructuralClustalAlignment structuralClustalResult = StructuralClustalAlignment [ClustalAlignmentEntry (T.pack "AB001721.1/2735-2851") (T.pack "CCCGGUGACUAUAGAGAGAGGGCCACACCCGUUCCCAUCCCGAACACGGAAGUUAAGCCUCUCAUCGCUGAUGGUACUAUGUGGUUCGCUGCAUGGGAGAGUAGGACGUUGCCGGGU"),ClustalAlignmentEntry (T.pack "gi|451991584:1-117") (T.pack "CCCGGUGACUAUAGAGAGAGGGCCACACCCGUUCCCAUCCCGAACACGGAAGUUAAGCCUCUCAUCGCUGAUGGUACUAUGUGGUUCGCUGCAUGGGAGAGUAGGACGUUGCCGGGU")] (T.pack "(((((((((....(.(((.(((.....))).))).)...(((....)))..(((..(((((((((.((.(.(((.(((....))))))).))))))))).))..)))))))))))).") (-38.30)