diff --git a/Biobase/Taxonomy.hs b/Biobase/Taxonomy.hs new file mode 100644 --- /dev/null +++ b/Biobase/Taxonomy.hs @@ -0,0 +1,13 @@ +-- | Functions for parsing, processing and visualization of taxonomy data. +-- +module Biobase.Taxonomy ( -- * Datatypes + -- Datatypes used to represent taxonomy data + module Biobase.Taxonomy.Types, + module Biobase.Taxonomy.Import, + module Biobase.Taxonomy.Utils, + module Biobase.Taxonomy.Visualization + ) where +import Biobase.Taxonomy.Types +import Biobase.Taxonomy.Import +import Biobase.Taxonomy.Utils +import Biobase.Taxonomy.Visualization diff --git a/Biobase/Taxonomy/Import.hs b/Biobase/Taxonomy/Import.hs new file mode 100644 --- /dev/null +++ b/Biobase/Taxonomy/Import.hs @@ -0,0 +1,432 @@ +-- | Functions for parsing, processing and visualization of taxonomy data. +-- +-- === Usage example: +-- * Read in taxonomy data +-- +-- > eitherTaxtree <- readNamedTaxonomy "/path/to/NCBI_taxonomydump_directory" +-- +-- * Process data +-- +-- > let subtree = extractTaxonomySubTreebyLevel [562] (fromRight eitherTaxTree) (Just 4) +-- +-- * Visualize result +-- +-- tput "/path/to/dotdirectory" subtree +module Biobase.Taxonomy.Import ( -- * Datatypes + -- Datatypes used to represent taxonomy data + module Biobase.Taxonomy.Types, + -- * Parsing + -- Functions prefixed with "read" read from filepaths, functions with parse from Haskell Strings. + readTaxonomy, + readNamedTaxonomy, + parseTaxonomy, + parseNCBITaxCitations, + readNCBITaxCitations, + parseNCBITaxDelNodes, + readNCBITaxDelNodes, + parseNCBITaxDivisions, + readNCBITaxDivisions, + parseNCBITaxGenCodes, + readNCBITaxGenCodes, + parseNCBITaxMergedNodes, + readNCBITaxMergedNodes, + parseNCBITaxNames, + readNCBITaxNames, + parseNCBITaxNodes, + readNCBITaxNodes, + parseNCBISimpleTaxons, + readNCBISimpleTaxons, + readNCBITaxonomyDatabase + ) where +import Prelude +import System.IO +import Biobase.Taxonomy.Types +import Text.Parsec.Prim (runP) +import Text.ParserCombinators.Parsec +import Control.Monad +import Data.List +import Data.Maybe +import qualified Data.Either.Unwrap as E +import Data.Graph.Inductive.Graph +import Data.Graph.Inductive.Tree +import qualified Data.ByteString.Char8 as B +import qualified Data.Text.Lazy as T +-------------------------------------------------------- + +--------------------------------------- +-- Parsing functions + +-- | NCBI taxonomy dump nodes and names in the input directory path are parsed and a SimpleTaxon tree with scientific names for each node is generated. +readNamedTaxonomy :: String -> IO (Either ParseError (Gr SimpleTaxon Double)) +readNamedTaxonomy directoryPath = do + nodeNames <- readNCBITaxNames (directoryPath ++ "names.dmp") + if E.isLeft nodeNames + then return (Left (E.fromLeft nodeNames)) + else do + let rightNodeNames = E.fromRight nodeNames + let filteredNodeNames = filter isScientificName rightNodeNames + let namedTaxonomyGraph = genParserNamedTaxonomyGraph filteredNodeNames + parseFromFileEncISO88591 namedTaxonomyGraph (directoryPath ++ "nodes.dmp") + +isScientificName :: TaxName -> Bool +isScientificName name = nameClass name == scientificNameT + where scientificNameT = B.pack "scientific name" + +-- | NCBI taxonomy dump nodes and names in the input directory path are parsed and a SimpleTaxon tree is generated. +readTaxonomy :: String -> IO (Either ParseError (Gr SimpleTaxon Double)) +readTaxonomy = parseFromFileEncISO88591 genParserTaxonomyGraph + +-- | NCBI taxonomy dump nodes and names in the input directory path are parsed and a SimpleTaxon tree is generated. +parseTaxonomy :: String -> Either ParseError (Gr SimpleTaxon Double) +parseTaxonomy = parse genParserTaxonomyGraph "parseTaxonomy" + +genParserTaxonomyGraph :: GenParser Char st (Gr SimpleTaxon Double) +genParserTaxonomyGraph = do + nodesEdges <- many1 (try genParserGraphNodeEdge) + optional eof + let (nodesList,edgesList) = unzip nodesEdges + --let taxedges = filter (\(a,b,_) -> a /= b) edgesList + let taxedges = filter notLoopEdge edgesList + --let taxnodes = concat nodesList + --return (mkGraph taxnodes taxedges) + let currentGraph = mkGraph nodesList taxedges + return currentGraph + + +notLoopEdge :: (Int,Int,a) -> Bool +notLoopEdge (a,b,_) = a /= b + +--genParserNodeEdges :: [TaxName] -> GenParser Char st [(Int,SimpleTaxon),(Int,Int,Double)] +--genParserNodeEdges = do +-- nodesEdges <- (many1 (try genParserGraphNodeEdge)) +-- optional eof +-- return (nodesList,edgesList) + + + --let taxedges = filter notLoopEdge edgesList + --let taxnamednodes = map (setNodeScientificName filteredNodeNames) nodesList + --let currentGraph = mkGraph taxnamednodes taxedges + --return currentGraph + +genParserNamedTaxonomyGraph :: [TaxName] -> GenParser Char st (Gr SimpleTaxon Double) +genParserNamedTaxonomyGraph filteredNodeNames = do + nodesEdges <- (many1 (try genParserGraphNodeEdge)) + optional eof + let (nodesList,edgesList) = unzip nodesEdges + let taxedges = filter notLoopEdge edgesList + let taxnamednodes = map (setNodeScientificName filteredNodeNames) nodesList + let currentGraph = mkGraph taxnamednodes taxedges + return currentGraph + +setNodeScientificName :: [TaxName] -> (t, SimpleTaxon) -> (t, SimpleTaxon) +setNodeScientificName inputTaxNames (inputNode,inputTaxon) = outputNode + where maybeRetrievedName = find (isTaxNameIdSimpleTaxid inputTaxon) inputTaxNames + retrievedName = maybe (T.pack "no name") nameTxt maybeRetrievedName + outputNode = (inputNode,inputTaxon{simpleScientificName = retrievedName}) + +isTaxNameIdSimpleTaxid :: SimpleTaxon -> TaxName -> Bool +isTaxNameIdSimpleTaxid inputTaxon inputTaxName = nameTaxId inputTaxName == simpleTaxId inputTaxon + + +genParserGraphNodeEdge :: GenParser Char st ((Int,SimpleTaxon),(Int,Int,Double)) +genParserGraphNodeEdge = do + _simpleTaxId <- many1 digit + string "\t|\t" + _simpleParentTaxId <- many1 digit + string "\t|\t" + _simpleRank <- many1 (noneOf "\t") + many1 (noneOf "\n") + char '\n' + let _simpleTaxIdInt = readInt _simpleTaxId + let _simpleParentTaxIdInt = readInt _simpleParentTaxId + return ((_simpleTaxIdInt,SimpleTaxon _simpleTaxIdInt T.empty _simpleParentTaxIdInt (readRank _simpleRank)),(_simpleTaxIdInt,_simpleParentTaxIdInt,1 :: Double)) + +-- | parse NCBITaxCitations from input string +parseNCBITaxCitations :: String -> Either ParseError [TaxCitation] +parseNCBITaxCitations = parse genParserNCBITaxCitations "parseTaxCitations" + +-- | parse NCBITaxCitations from input filePath +readNCBITaxCitations :: String -> IO (Either ParseError [TaxCitation]) +readNCBITaxCitations = parseFromFileEncISO88591 genParserNCBITaxCitations + +-- | parse NCBITaxDelNodes from input string +parseNCBITaxDelNodes :: String -> Either ParseError [TaxDelNode] +parseNCBITaxDelNodes = parse genParserNCBITaxDelNodes "parseTaxDelNodes" + +-- | parse NCBITaxDelNodes from input filePath +readNCBITaxDelNodes :: String -> IO (Either ParseError [TaxDelNode]) +readNCBITaxDelNodes = parseFromFile genParserNCBITaxDelNodes + +-- | parse NCBITaxDivisons from input string +parseNCBITaxDivisions :: String -> Either ParseError [TaxDivision] +parseNCBITaxDivisions = parse genParserNCBITaxDivisons "parseTaxDivisons" + +-- | parse NCBITaxDivisons from input filePath +readNCBITaxDivisions :: String -> IO (Either ParseError [TaxDivision]) +readNCBITaxDivisions = parseFromFile genParserNCBITaxDivisons + +-- | parse NCBITaxGenCodes from input string +parseNCBITaxGenCodes :: String -> Either ParseError [TaxGenCode] +parseNCBITaxGenCodes = parse genParserNCBITaxGenCodes "parseTaxGenCodes" + +-- | parse NCBITaxGenCodes from input filePath +readNCBITaxGenCodes :: String -> IO (Either ParseError [TaxGenCode]) +readNCBITaxGenCodes = parseFromFile genParserNCBITaxGenCodes + +-- | parse NCBITaxMergedNodes from input string +parseNCBITaxMergedNodes :: String -> Either ParseError [TaxMergedNode] +parseNCBITaxMergedNodes = parse genParserNCBITaxMergedNodes "parseTaxMergedNodes" + +-- | parse NCBITaxMergedNodes from input filePath +readNCBITaxMergedNodes :: String -> IO (Either ParseError [TaxMergedNode]) +readNCBITaxMergedNodes = parseFromFile genParserNCBITaxMergedNodes + +-- | parse NCBITaxNames from input string +parseNCBITaxNames :: String -> Either ParseError [TaxName] +parseNCBITaxNames = parse genParserNCBITaxNames "parseTaxNames" + +-- | parse NCBITaxNames from input filePath +readNCBITaxNames :: String -> IO (Either ParseError [TaxName]) +readNCBITaxNames = parseFromFile genParserNCBITaxNames + +-- | parse NCBITaxNames from input string +parseNCBITaxNodes :: String -> Either ParseError TaxNode +parseNCBITaxNodes = parse genParserNCBITaxNode "parseTaxNode" + +-- | parse NCBITaxCitations from input filePath +readNCBITaxNodes :: String -> IO (Either ParseError [TaxNode]) +readNCBITaxNodes = parseFromFile genParserNCBITaxNodes + +-- | parse NCBISimpleTaxNames from input string +parseNCBISimpleTaxons :: String -> Either ParseError SimpleTaxon +parseNCBISimpleTaxons = parse genParserNCBISimpleTaxon "parseSimpleTaxon" + +-- | parse NCBITaxCitations from input filePath +readNCBISimpleTaxons :: String -> IO (Either ParseError [SimpleTaxon]) +readNCBISimpleTaxons = parseFromFile genParserNCBISimpleTaxons + +-- | Parse the input as NCBITax datatype +readNCBITaxonomyDatabase :: String -> IO (Either [String] NCBITaxDump) +readNCBITaxonomyDatabase folder = do + citations <- readNCBITaxCitations (folder ++ "citations.dmp") + let citationsError = extractParseError citations + taxdelNodes <- readNCBITaxDelNodes (folder ++ "delnodes.dmp") + let delNodesError = extractParseError taxdelNodes + divisons <- readNCBITaxDivisions (folder ++ "division.dmp") + let divisonsError = extractParseError divisons + genCodes <- readNCBITaxGenCodes (folder ++ "gencode.dmp") + let genCodesError = extractParseError genCodes + mergedNodes <- readNCBITaxMergedNodes (folder ++ "merged.dmp") + let mergedNodesError = extractParseError mergedNodes + names <- readNCBITaxNames (folder ++ "names.dmp") + let namesError = extractParseError names + taxnodes <- readNCBITaxNodes (folder ++ "nodes.dmp") + let nodesError = extractParseError taxnodes + let parseErrors = [citationsError, delNodesError, divisonsError, genCodesError, mergedNodesError, namesError, nodesError] + return (checkParsing parseErrors citations taxdelNodes divisons genCodes mergedNodes names taxnodes) + +genParserNCBITaxCitations :: GenParser Char st [TaxCitation] +genParserNCBITaxCitations = many1 genParserNCBITaxCitation + +genParserNCBITaxDelNodes :: GenParser Char st [TaxDelNode] +genParserNCBITaxDelNodes = many1 genParserNCBITaxDelNode + +genParserNCBITaxDivisons :: GenParser Char st [TaxDivision] +genParserNCBITaxDivisons = many1 genParserNCBITaxDivision + +genParserNCBITaxGenCodes :: GenParser Char st [TaxGenCode] +genParserNCBITaxGenCodes = many1 genParserNCBITaxGenCode + + +genParserNCBITaxMergedNodes :: GenParser Char st [TaxMergedNode] +genParserNCBITaxMergedNodes = many1 genParserNCBITaxMergedNode + + +genParserNCBITaxNames :: GenParser Char st [TaxName] +genParserNCBITaxNames = many1 genParserNCBITaxName + +genParserNCBITaxNodes :: GenParser Char st [TaxNode] +genParserNCBITaxNodes = many1 genParserNCBITaxNode + +genParserNCBISimpleTaxons :: GenParser Char st [SimpleTaxon] +genParserNCBISimpleTaxons = many1 genParserNCBISimpleTaxon + + +genParserNCBITaxCitation :: GenParser Char st TaxCitation +genParserNCBITaxCitation = do + _citId <- many1 digit + string "\t|\t" + _citKey <- many (noneOf "\t") + string "\t|\t" + _pubmedId <- optionMaybe (many1 digit) + string "\t|\t" + _medlineId <- optionMaybe (many1 digit) + tab + char '|' + _url <- genParserTaxURL + char '|' + tab + _text <- (many (noneOf "\t")) + string "\t|\t" + _taxIdList <- (many genParserTaxIdList) + string "\t|\n" + return $ TaxCitation (readInt _citId) (B.pack _citKey) (liftM readInt _pubmedId) (liftM readInt _medlineId) _url (B.pack _text) _taxIdList + +genParserNCBITaxDelNode :: GenParser Char st TaxDelNode +genParserNCBITaxDelNode = do + taxdelNode <- many1 digit + space + char '|' + char '\n' + return $ TaxDelNode (readInt taxdelNode) + +genParserNCBITaxDivision :: GenParser Char st TaxDivision +genParserNCBITaxDivision = do + _divisionId <- many1 digit + string "\t|\t" + _divisionCDE <- many1 upper + string "\t|\t" + _divisionName <- many1 (noneOf "\t") + string "\t|\t" + _comments <- many1 (noneOf "\t") + string "\t|\n" + return $ TaxDivision (readInt _divisionId) (B.pack _divisionCDE) (B.pack _divisionName) (B.pack _comments) + +genParserNCBITaxGenCode :: GenParser Char st TaxGenCode +genParserNCBITaxGenCode = do + _geneticCodeId <- many1 digit + string "\t|\t" + _abbreviation <- (many1 (noneOf "\t")) + string "\t|\t" + _genCodeName <- many1 (noneOf "\t") + string "\t|\t" + _cde <- many1 (noneOf "\t") + string "\t|\t" + _starts <- many1 (noneOf "\t") + string "\t|\n" + return $ TaxGenCode (readInt _geneticCodeId) (B.pack _abbreviation) (B.pack _genCodeName) (B.pack _cde) (B.pack _starts) + +genParserNCBITaxMergedNode :: GenParser Char st TaxMergedNode +genParserNCBITaxMergedNode = do + _oldTaxId <- many1 digit + string "\t|\t" + _newTaxId <- many1 digit + string "\t|\n" + return $ TaxMergedNode (readInt _oldTaxId) (readInt _newTaxId) + +genParserNCBITaxName :: GenParser Char st TaxName +genParserNCBITaxName = do + _taxId <- many1 digit + string "\t|\t" + _nameTxt <- many1 (noneOf "\t\n") + string "\t|\t" + _uniqueName <- many (noneOf "\t\n") + string "\t|\t" + _nameClass <- many1 (noneOf "\t\n") + tab + char '|' + newline + return $! TaxName (readInt _taxId) (T.pack _nameTxt) (B.pack _uniqueName) (B.pack _nameClass) + +genParserNCBISimpleTaxon :: GenParser Char st SimpleTaxon +genParserNCBISimpleTaxon = do + _simpleTaxId <- many1 digit + string "\t|\t" + _simpleParentTaxId <- many1 digit + string "\t|\t" + _simpleRank <- many1 (noneOf "\t") + many1 (noneOf "\n") + char '\n' + return $! SimpleTaxon (readInt _simpleTaxId) T.empty (readInt _simpleParentTaxId) (readRank _simpleRank) + +genParserNCBITaxNode :: GenParser Char st TaxNode +genParserNCBITaxNode = do + _taxId <- many1 digit + string "\t|\t" + _parentTaxId <- many1 digit + string "\t|\t" + _rank <- many1 (noneOf "\t") + string "\t|\t" + _emblCode <- (many (noneOf "\t")) + string "\t|\t" + _divisionId <- many1 digit + string "\t|\t" + _inheritedDivFlag <- many1 digit + string "\t|\t" + _geneticCodeId <- many1 digit + string "\t|\t" + _inheritedGCFlag <- many1 digit + string "\t|\t" + _mitochondrialGeneticCodeId <- many1 digit + string "\t|\t" + _inheritedMGCFlag <- many1 digit + string "\t|\t" + _genBankHiddenFlag <- many1 digit + string "\t|\t" + _hiddenSubtreeRootFlag <- many1 digit + string "\t|\t" + _comments <- many (noneOf "\t") + tab + char '|' + char '\n' + return $ TaxNode (readInt _taxId) (readInt _parentTaxId) (readRank _rank) (B.pack _emblCode) (read _divisionId :: Int) (readBool _inheritedDivFlag) (read _geneticCodeId ::Int) (readBool _inheritedGCFlag) (read _mitochondrialGeneticCodeId ::Int) (readBool _inheritedMGCFlag) (readBool _genBankHiddenFlag) (readBool _hiddenSubtreeRootFlag) (B.pack _comments) + +--------------------------------------- +-- Auxiliary functions +readInt :: String -> Int +readInt = read + +readBool :: String -> Bool +readBool "0" = False +readBool "1" = True +readBool _ = False + +readRank :: String -> Rank +readRank a = read a :: Rank + +genParserTaxIdList :: GenParser Char st Int +genParserTaxIdList = do + optional (char ' ') + _taxId <- many1 digit + optional (char ' ') + return (readInt _taxId) + +genParserTaxURL :: GenParser Char st B.ByteString +genParserTaxURL = do + tab + url1 <- many (noneOf "\t") + tab + url2 <- many (noneOf "|") + return (B.pack (url1 ++ url2)) + --return (concatenateURLParts url1 url2) + +concatenateURLParts :: Maybe String -> Maybe String -> Maybe String +concatenateURLParts url1 url2 + | isJust url1 && isJust url2 = maybeStringConcat url1 url2 + | isJust url1 && isNothing url2 = url1 + | otherwise = Nothing + +maybeStringConcat :: Maybe String -> Maybe String -> Maybe String +maybeStringConcat = liftM2 (++) + +readEncodedFile :: TextEncoding -> FilePath -> IO String +readEncodedFile encoding name = do + handle <- openFile name ReadMode + hSetEncoding handle encoding + hGetContents handle + +parseFromFileEncISO88591 :: Parser a -> String -> IO (Either ParseError a) +parseFromFileEncISO88591 parser fname = do + input <- readEncodedFile latin1 fname + return (runP parser () fname input) + +-- | check a list of parsing results for presence of Left aka Parse error +checkParsing :: [String] -> Either ParseError [TaxCitation] -> Either ParseError [TaxDelNode] -> Either ParseError [TaxDivision] -> Either ParseError [TaxGenCode] -> Either ParseError [TaxMergedNode] -> Either ParseError [TaxName] -> Either ParseError [TaxNode]-> Either [String] NCBITaxDump +checkParsing parseErrors citations taxdelNodes divisons genCodes mergedNodes names taxnodes + | join parseErrors == "" = Right (NCBITaxDump (E.fromRight citations) (E.fromRight taxdelNodes) (E.fromRight divisons) (E.fromRight genCodes) (E.fromRight mergedNodes) (E.fromRight names) (E.fromRight taxnodes)) + | otherwise = Left parseErrors + +extractParseError :: Either ParseError a -> String +extractParseError _parse + | E.isLeft _parse = show (E.fromLeft _parse) + | otherwise = "" diff --git a/Biobase/Taxonomy/Types.hs b/Biobase/Taxonomy/Types.hs new file mode 100644 --- /dev/null +++ b/Biobase/Taxonomy/Types.hs @@ -0,0 +1,303 @@ +-- | This module contains data structures for +-- taxonomy data + +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE DuplicateRecordFields #-} + +module Biobase.Taxonomy.Types where +import Prelude +import qualified Data.ByteString.Char8 as B +import qualified Data.Aeson as A +import qualified Data.Vector as V +--import Data.Graph.Inductive +import Data.Graph.Inductive.Graph +import Data.Graph.Inductive.Tree +import qualified Data.Text as T +import qualified Data.Text.Lazy as TL +--import qualified Data.Text.Encoding + +-- | SimpleTaxon only contains the most relevant fields of a taxonomy entry. +-- For all annotaded fields use the Taxon datatype and its associated functions +data SimpleTaxon = SimpleTaxon + { + -- node id in GenBank + simpleTaxId :: Int, + simpleScientificName :: TL.Text, + -- parent node id in GenBank taxonomy database + simpleParentTaxId :: Int, + -- rank of this node (superkingdom, kingdom, ...) + simpleRank :: Rank + } + deriving (Show, Read, Eq) + +-- | Datastructure for tree comparisons +data CompareTaxon = CompareTaxon + { + compareScientificName :: TL.Text, + compareRank :: Rank, + -- number indicating in which trees, + inTree :: [Int] + } + deriving (Show, Read, Eq) + +-- | Data structure for Entrez taxonomy fetch result +data Taxon = Taxon + { taxonTaxId :: Int + , taxonScientificName :: B.ByteString + , taxonParentTaxId :: Int + , taxonRank :: Rank + , division :: B.ByteString + , geneticCode :: TaxGenCode + , mitoGeneticCode :: TaxGenCode + , lineage :: B.ByteString + , lineageEx :: [LineageTaxon] + , createDate :: B.ByteString + , updateDate :: B.ByteString + , pubDate :: B.ByteString + } deriving (Show, Eq) + + +data TaxonName = TaxonName + { classCDE :: B.ByteString + , dispName :: B.ByteString + } deriving (Show, Eq) + +-- | Lineage Taxons denote all parent Taxonomy nodes of a node retrieved by Entrez fetch +data Lineage = Lineage + { lineageStartTaxId :: Int + , lineageStartScienticName :: B.ByteString + , lineageStartRank :: Rank + , lineageTaxons :: [LineageTaxon] + } + deriving (Show, Eq) + + +-- | Lineage Taxons denote all parent Taxonomy nodes of a node retrieved by Entrez fetch +data LineageTaxon = LineageTaxon + { lineageTaxId :: Int + , lineageScienticName :: B.ByteString + , lineageRank :: Rank} + deriving (Show, Eq) + +-- | NCBI Taxonomy database dump hierachichal data structure +-- as defined in ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump_readme.txt +data NCBITaxDump = NCBITaxDump + { + taxCitations :: [TaxCitation], + taxDelNodes :: [TaxDelNode], + taxDivisions :: [TaxDivision], + taxGenCodes :: [TaxGenCode], + taxMergedNodes :: [TaxMergedNode], + taxNames :: [TaxName], + taxNodes :: [TaxNode] + } + deriving (Show, Read, Eq) + +-- | Datastructure for entries of Taxonomy database dump citations file +data TaxCitation = TaxCitation + { + -- the unique id of citation + citId :: Int, + -- citation key + citKey :: B.ByteString, + -- unique id in PubMed database (0 if not in PubMed) + pubmedId :: Maybe Int, + -- unique id in MedLine database (0 if not in MedLine) + medlineId :: Maybe Int, + -- URL associated with citation + url :: B.ByteString, + -- any text (usually article name and authors) + -- The following characters are escaped in this text by a backslash: + -- newline (appear as "\n"), + -- tab character ("\t"), + -- double quotes ('\"'), + -- backslash character ("\\"). + text :: B.ByteString, + -- list of node ids separated by a single space + taxIdList :: [Int] + } + deriving (Show, Read, Eq) + +-- | Datastructure for entries of Taxonomy database dump deleted nodes file +data TaxDelNode = TaxDelNode + { + -- deleted node id + delTaxId :: Int + } + deriving (Show, Read, Eq) + +-- | Datastructure for entries of Taxonomy database dump division file +data TaxDivision = TaxDivision + { + -- taxonomy database division id + divisionId :: Int, + -- GenBank division code (three characters) + divisionCDE :: B.ByteString, + -- e.g. BCT, PLN, VRT, MAM, PRI... + divisonName :: B.ByteString, + divisionComments :: B.ByteString + } + deriving (Show, Read, Eq) + +-- | Datastructure for entries of Taxonomy database dump gencode file +data TaxGenCode = TaxGenCode + { + -- GenBank genetic code id + geneticCodeId :: Int, + -- genetic code name abbreviation + abbreviation :: B.ByteString, + -- genetic code name + geneCodeName :: B.ByteString, + -- translation table for this genetic code + cde :: B.ByteString, + -- start codons for this genetic code + starts :: B.ByteString + } + deriving (Show, Read, Eq) + +-- | Datastructure for entries of Taxonomy database dump mergednodes file +data TaxMergedNode = TaxMergedNode + { + -- id of nodes which has been merged + oldTaxId :: Int, + -- id of nodes which is result of merging + newTaxId :: Int + } + deriving (Show, Read, Eq) + +-- | Datastructure for entries of Taxonomy database dump names file +data TaxName = TaxName + { + -- the id of node associated with this name + nameTaxId :: Int, + -- name itself + nameTxt :: TL.Text, + -- the unique variant of this name if name not unique + uniqueName :: B.ByteString, + -- (synonym, common name, ...) + nameClass :: B.ByteString + } + deriving (Show, Read, Eq) + +-- | Taxonomic ranks: NCBI uses the uncommon Speciessubgroup +data Rank = Norank | Form | Variety | Infraspecies | Subspecies | Speciessubgroup | Species | Speciesgroup | Superspecies | Series | Section | Subgenus | Genus | Subtribe | Tribe | Supertribe | Subfamily | Family | Superfamily | Parvorder | Infraorder | Suborder | Order | Superorder | Magnorder | Cohort | Legion | Parvclass | Infraclass | Subclass | Class | Superclass | Microphylum | Infraphylum | Subphylum | Phylum | Superphylum | Infrakingdom | Subkingdom | Kingdom | Superkingdom | Domain deriving (Eq, Ord, Show, Bounded, Enum) + +readsRank :: String -> [(Rank, String)] +instance Read Rank where + readsPrec _ = readsRank + +readsRank input -- = [(Domain x)| x <- reads input ] + | input == "domain" = [(Domain,"")] + | input == "superkingdom" = [(Superkingdom,"")] + | input == "kingdom" = [(Kingdom,"")] + | input == "subkingdom" = [(Subkingdom,"")] + | input == "infrakingdom" = [(Infrakingdom,"")] + | input == "superphylum" = [(Superphylum,"")] + | input == "phylum" = [(Phylum,"")] + | input == "subphylum" = [(Subphylum,"")] + | input == "infraphylum" = [(Infraphylum,"")] + | input == "microphylum" = [(Microphylum,"")] + | input == "superclass" = [(Superclass,"")] + | input == "class" = [(Class,"")] + | input == "subclass" = [(Subclass,"")] + | input == "infraclass" = [(Infraclass,"")] + | input == "parvclass " = [(Parvclass ,"")] + | input == "legion" = [(Legion,"")] + | input == "cohort" = [(Cohort,"")] + | input == "magnorder " = [(Magnorder ,"")] + | input == "superorder" = [(Superorder,"")] + | input == "order" = [(Order,"")] + | input == "suborder" = [(Suborder,"")] + | input == "infraorder" = [(Infraorder,"")] + | input == "parvorder" = [(Parvorder,"")] + | input == "superfamily" = [(Superfamily,"")] + | input == "family" = [(Family,"")] + | input == "subfamily" = [(Subfamily,"")] + | input == "supertribe" = [(Supertribe,"")] + | input == "tribe" = [(Tribe,"")] + | input == "subtribe" = [(Subtribe,"")] + | input == "genus" = [(Genus,"")] + | input == "subgenus" = [(Subgenus,"")] + | input == "section" = [(Section,"")] + | input == "series" = [(Series,"")] + | input == "superspecies" = [(Superspecies,"")] + | input == "species group" = [(Speciesgroup,"")] + | input == "species" = [(Species,"")] + | input == "species subgroup" = [(Speciessubgroup,"")] + | input == "subspecies" = [(Subspecies,"")] + | input == "infraspecies" = [(Infraspecies,"")] + | input == "varietas" = [(Variety,"")] + | input == "forma" = [(Form,"")] + | input == "no rank" = [(Norank,"")] + | otherwise = [(Norank,"")] + +-- | Datastructure for entries of Taxonomy database dump nodes file +data TaxNode = TaxNode + { + -- node id in GenBank + taxId :: Int, + -- parent node id in GenBank taxonomy database + parentTaxId :: Int, + -- rank of this node (superkingdom, kingdom, ...) + rank :: Rank, + -- locus-name prefix; not unique + emblCode :: B.ByteString, + -- see division.dmp file + nodeDivisionId :: Int, + -- 1 if node inherits division from parent + inheritedDivFlag :: Bool, + -- see gencode.dmp file + nodeGeneticCodeId :: Int, + -- 1 if node inherits genetic code from parent + inheritedGCFlag :: Bool, + -- see gencode.dmp file + mitochondrialGeneticCodeId :: Int, + -- 1 if node inherits mitochondrial gencode from parent + inheritedMGCFlag :: Bool, + -- 1 if name is suppressed in GenBank entry lineage + genBankHiddenFlag :: Bool, + -- 1 if this subtree has no sequence data yet + hiddenSubtreeRootFlag :: Bool, + -- free-text comments and citations + nodeComments :: B.ByteString + } + deriving (Show, Read, Eq) + +-- | Simple Gene2Accession table +data SimpleGene2Accession = SimpleGene2Accession + { simpleTaxIdEntry :: Int, + simpleGenomicNucleotideAccessionVersion :: B.ByteString + } deriving (Show, Eq, Read) + +-- | Datastructure for Gene2Accession table +data Gene2Accession = Gene2Accession + { taxIdEntry :: Int, + geneID :: Int, + status :: B.ByteString, + rnaNucleotideAccessionVersion :: B.ByteString, + rnaNucleotideGi :: B.ByteString, + proteinAccessionVersion :: B.ByteString, + proteinGi :: B.ByteString, + genomicNucleotideAccessionVersion :: B.ByteString, + genomicNucleotideGi :: B.ByteString, + startPositionOnTheGenomicAccession :: B.ByteString, + endPositionOnTheGenomicAccession :: B.ByteString, + orientation :: B.ByteString, + assembly :: B.ByteString, + maturePeptideAccessionVersion :: B.ByteString, + maturePeptideGi :: B.ByteString + } deriving (Show, Eq, Read) + +instance A.ToJSON (Gr SimpleTaxon Double) where + toJSON inputGraph = simpleTaxonJSONValue inputGraph 1 + +simpleTaxonJSONValue :: Gr SimpleTaxon Double -> Node -> A.Value +simpleTaxonJSONValue inputGraph node = jsonValue + where jsonValue = A.object [currentScientificName,T.pack "children" A..= children] + childNodes = suc inputGraph node + currentLabel = lab inputGraph node + currentScientificName = T.pack "name" A..= maybe (T.pack "notFound") (TL.toStrict . simpleScientificName) currentLabel + children = A.Array (V.fromList (map (simpleTaxonJSONValue inputGraph) childNodes)) + --jsonValue = A.object [currentScientificName,currentId,currentRank,(T.pack "children") A..= children] + --currentId = (T.pack "id") A..= (maybe (T.pack "notFound") (\a -> T.pack (show (simpleTaxId a))) currentLabel) + --currentRank = (T.pack "rank") A..= (maybe (T.pack "notFound") (\a -> T.pack (show (simpleRank a))) currentLabel) diff --git a/Biobase/Taxonomy/Utils.hs b/Biobase/Taxonomy/Utils.hs new file mode 100644 --- /dev/null +++ b/Biobase/Taxonomy/Utils.hs @@ -0,0 +1,142 @@ +-- | Functions for processing of taxonomy data. +-- +module Biobase.Taxonomy.Utils ( -- * Datatypes + -- Datatypes used to represent taxonomy data + module Biobase.Taxonomy.Types, + -- * Processing + compareSubTrees, + extractTaxonomySubTreebyLevel, + extractTaxonomySubTreebyLevelNew, + extractTaxonomySubTreebyRank, + safeNodePath, + getParentbyRank, + ) where +import Prelude +import Biobase.Taxonomy.Types +import Data.List +import qualified Data.Vector as V +import Data.Maybe +import Data.Graph.Inductive.Graph +import Data.Graph.Inductive.Query.SP (sp) +import Data.Graph.Inductive.Query.BFS (level) +import Data.Graph.Inductive.Tree +import Data.Graph.Inductive.Basic + +--------------------------------------- +-- Processing functions + +-- | Extract a subtree correpsonding to input node paths to root. Only nodes in level number distance to root are included. Used in Ids2TreeCompare tool. +compareSubTrees :: [Gr SimpleTaxon Double] -> (Int,Gr CompareTaxon Double) +compareSubTrees graphs = (length graphs,resultGraph) + where treesLabNodes = map labNodes graphs + treesLabEdges = map labEdges graphs + mergedNodes = nub (concat treesLabNodes) + mergedEdges = nub (concat treesLabEdges) + --annotate node in which of the compared trees they are present + comparedNodes = annotateTaxonsDifference treesLabNodes mergedNodes + resultGraph = mkGraph comparedNodes mergedEdges :: Gr CompareTaxon Double + +annotateTaxonsDifference :: [[LNode SimpleTaxon]] -> [LNode SimpleTaxon] -> [LNode CompareTaxon] +annotateTaxonsDifference treesNodes mergedtreeNodes = comparedNodes + where comparedNodes = map (annotateTaxonDifference indexedTreesNodes) mergedtreeNodes + indexedTreesNodes = zip [0..(length treesNodes)] treesNodes + + +annotateTaxonDifference :: [(Int,[LNode SimpleTaxon])] -> LNode SimpleTaxon -> LNode CompareTaxon +annotateTaxonDifference indexedTreesNodes mergedtreeNode = comparedNode + where comparedNode = (simpleTaxId (snd mergedtreeNode),CompareTaxon (simpleScientificName (snd mergedtreeNode)) (simpleRank (snd mergedtreeNode)) currentInTree) + currentInTree = concatMap (\(i,treeNodes) -> [i | mergedtreeNode `elem` treeNodes]) indexedTreesNodes + +-- | Extract a subtree corresponding to input node paths to root. Only nodes in level number distance to root are included. Used in Ids2Tree tool. +extractTaxonomySubTreebyLevel :: [Node] -> Gr SimpleTaxon Double -> Maybe Int -> Gr SimpleTaxon Double +extractTaxonomySubTreebyLevel inputNodes graph levelNumber = taxonomySubTree + where paths = nub (concatMap (getPath (1 :: Node) graph) inputNodes) + contexts = map (context graph) paths + lnodes = map labNode' contexts + ledges = nub (concatMap (out graph . fst) lnodes) + unfilteredTaxonomySubTree = mkGraph lnodes ledges :: Gr SimpleTaxon Double + filteredLNodes = filterNodesByLevel levelNumber lnodes unfilteredTaxonomySubTree + filteredledges = nub (concatMap (out graph . fst) filteredLNodes) + taxonomySubTree = mkGraph filteredLNodes filteredledges :: Gr SimpleTaxon Double + +-- | Extract a subtree corresponding to input node paths to root. Only nodes in level number distance to root are included. Used in Ids2Tree tool. +extractTaxonomySubTreebyLevelNew :: [Node] -> Gr SimpleTaxon Double -> Maybe Int -> Gr SimpleTaxon Double +extractTaxonomySubTreebyLevelNew inputNodes graph levelNumber = taxonomySubTree + where inputNodeVector = V.fromList inputNodes + paths = V.concatMap (getVectorPath (1 :: Node) graph) inputNodeVector + contexts = V.map (context graph) paths + vlnodes = V.map labNode' contexts + ledges = concatMap (out graph . fst) lnodes + lnodes = V.toList vlnodes + --ledges = V.toList vledges + unfilteredTaxonomySubTree = mkGraph lnodes ledges :: Gr SimpleTaxon Double + filteredLNodes = filterNodesByLevel levelNumber lnodes unfilteredTaxonomySubTree + --filteredLNodesVector = V.fromList filteredLNodes + filteredledges = concatMap (out graph . fst) filteredLNodes + --filteredledges = V.toList filteredledgesVector + taxonomySubTree = mkGraph filteredLNodes filteredledges :: Gr SimpleTaxon Double + +-- | Extract a subtree corresponding to input node paths to root. If a Rank is provided, all node that are less or equal are omitted +extractTaxonomySubTreebyRank :: [Node] -> Gr SimpleTaxon Double -> Maybe Rank -> Gr SimpleTaxon Double +extractTaxonomySubTreebyRank inputNodes graph highestRank = taxonomySubTree + where paths = nub (concatMap (getPath (1 :: Node) graph) inputNodes) + contexts = map (context graph) paths + lnodes = map labNode' contexts + filteredLNodes = filterNodesByRank highestRank lnodes + filteredledges = nub (concatMap (out graph . fst) filteredLNodes) + taxonomySubTree = mkGraph filteredLNodes filteredledges :: Gr SimpleTaxon Double + +getVectorPath :: Node -> Gr SimpleTaxon Double -> Node -> V.Vector Node +getVectorPath root graph node = maybe V.empty V.fromList (sp node root graph) + +getPath :: Node -> Gr SimpleTaxon Double -> Node -> Path +getPath root graph node = maybe [] id (sp node root graph) + +-- | Extract parent node with specified Rank +getParentbyRank :: Node -> Gr SimpleTaxon Double -> Maybe Rank -> Maybe (Node, SimpleTaxon) +getParentbyRank inputNode graph requestedRank = filteredLNode + where path = maybe [] id (sp (inputNode :: Node) (1 :: Node) graph) + nodeContext = map (context graph) path + lnode = map labNode' nodeContext + filteredLNode = findNodeByRank requestedRank lnode + +-- | Filter nodes by distance from root +filterNodesByLevel :: Maybe Int -> [(Node, SimpleTaxon)] -> Gr SimpleTaxon Double -> [(Node, SimpleTaxon)] +filterNodesByLevel levelNumber inputNodes graph + | isJust levelNumber = filteredNodes + | otherwise = inputNodes + --distances of all nodes to root + where nodedistances = level (1::Node) (undir graph) + sortedNodeDistances = sortBy sortByNodeID nodedistances + sortedInputNodes = sortBy sortByNodeID inputNodes + zippedNodeDistancesInputNodes = zip sortedNodeDistances sortedInputNodes + zippedFilteredNodes = filter (\((_,distance),(_,_)) -> distance <= fromJust levelNumber) zippedNodeDistancesInputNodes + filteredNodes = map snd zippedFilteredNodes + +sortByNodeID :: (Node,a) -> (Node,a) -> Ordering +sortByNodeID (n1, _) (n2, _) + | n1 < n2 = GT + | n1 > n2 = LT + | n1 == n2 = EQ + | otherwise = EQ + +-- | Find only taxons of a specific rank in a list of input taxons +findNodeByRank :: Maybe Rank -> [(t, SimpleTaxon)] -> Maybe (t, SimpleTaxon) +findNodeByRank requestedRank inputNodes + | isJust requestedRank = filteredNodes + | otherwise = Nothing + where filteredNodes = find (\(_,t) -> simpleRank t == fromJust requestedRank) inputNodes + +-- | Filter a list of input taxons for a minimal provided rank +filterNodesByRank :: Maybe Rank -> [(t, SimpleTaxon)] -> [(t, SimpleTaxon)] +filterNodesByRank highestRank inputNodes + | isJust highestRank = filteredNodes + | otherwise = inputNodes + where filteredNodes = filter (\(_,t) -> simpleRank t >= fromJust highestRank) inputNodes ++ noRankNodes + noRankNodes = filter (\(_,t) -> simpleRank t == Norank) inputNodes + +-- | Returns path between 2 maybe nodes. Used in TreeDistance tool. +safeNodePath :: Maybe Node -> Gr SimpleTaxon Double -> Maybe Node -> Either String Path +safeNodePath nodeid1 graphOutput nodeid2 + | isJust nodeid1 && isJust nodeid2 = Right (maybe [] id (sp (fromJust nodeid1) (fromJust nodeid2) (undir graphOutput))) + | otherwise = Left "Both taxonomy ids must be provided for distance computation" diff --git a/Biobase/Taxonomy/Visualization.hs b/Biobase/Taxonomy/Visualization.hs new file mode 100644 --- /dev/null +++ b/Biobase/Taxonomy/Visualization.hs @@ -0,0 +1,100 @@ +-- | Functions for visualization of taxonomy data. +module Biobase.Taxonomy.Visualization ( -- * Datatypes + -- Datatypes used to represent taxonomy data + module Biobase.Taxonomy.Types, + -- * Visualization + drawTaxonomyComparison, + drawTaxonomy, + writeTree, + writeDotTree, + writeJsonTree + ) where +import Prelude +import Biobase.Taxonomy.Types +import Data.Graph.Inductive.Tree +import Data.Graph.Inductive.Basic +import qualified Data.GraphViz as GV +import qualified Data.GraphViz.Printing as GVP +import qualified Data.GraphViz.Attributes.Colors as GVAC +import qualified Data.GraphViz.Attributes.Complete as GVA +import qualified Data.ByteString.Lazy.Char8 as L +import qualified Data.Aeson as AE +import qualified Data.Text.Lazy as T + +--------------------------------------- +-- Visualisation functions + +-- | Draw graph in dot format. Used in Ids2Tree tool. +drawTaxonomy :: Bool -> Gr SimpleTaxon Double -> String +drawTaxonomy withRank inputGraph = do + let nodeFormating = if withRank then nodeFormatWithRank else nodeFormatWithoutRank + let params = GV.nonClusteredParams {GV.isDirected = True + , GV.globalAttributes = [GV.GraphAttrs [GVA.Size (GVA.GSize (20 :: Double) (Just (20 :: Double)) False)]] + , GV.isDotCluster = const True + --, GV.fmtNode = \ (_,l) -> [GV.textLabel (TL.pack (show (simpleRank l) ++ "\n" ++ T.unpack (simpleScientificName l)))] + , GV.fmtNode = nodeFormating + , GV.fmtEdge = const [] + } + let dotFormat = GV.graphToDot params inputGraph + let dottext = GVP.renderDot $ GVP.toDot dotFormat + T.unpack dottext + +nodeFormatWithRank :: (t, SimpleTaxon) -> [GVA.Attribute] +nodeFormatWithRank (_,l) = [GV.textLabel (T.concat [T.pack (show (simpleRank l)), T.pack ("\n") , simpleScientificName l])] + +nodeFormatWithoutRank :: (t, SimpleTaxon) -> [GVA.Attribute] +nodeFormatWithoutRank (_,l) = [GV.textLabel (simpleScientificName l)] + +-- | Draw tree comparison graph in dot format. Used in Ids2TreeCompare tool. +drawTaxonomyComparison :: Bool -> (Int,Gr CompareTaxon Double) -> String +drawTaxonomyComparison withRank (treeNumber,inputGraph) = do + let cList = makeColorList treeNumber + let nodeFormating = if withRank then (compareNodeFormatWithRank cList) else (compareNodeFormatWithoutRank cList) + let params = GV.nonClusteredParams {GV.isDirected = True + , GV.globalAttributes = [] + , GV.isDotCluster = const True + --, GV.fmtNode = \ (_,l) -> [GV.textLabel (TL.pack (show (compareRank l) ++ "\n" ++ B.unpack (compareScientificName l))), GV.style GV.wedged, GVA.Color (selectColors (inTree l) cList)] + , GV.fmtNode = nodeFormating + , GV.fmtEdge = const [] + } + let dotFormat = GV.graphToDot params (grev inputGraph) + let dottext = GVP.renderDot $ GVP.toDot dotFormat + T.unpack dottext + +compareNodeFormatWithRank :: [GVA.Color] -> (t, CompareTaxon) -> [GVA.Attribute] +compareNodeFormatWithRank cList (_,l) = [GV.textLabel (T.concat [T.pack (show (compareRank l) ++ "\n"),compareScientificName l]), GV.style GV.wedged, GVA.Color (selectColors (inTree l) cList)] + +compareNodeFormatWithoutRank :: [GVA.Color] -> (t, CompareTaxon) -> [GVA.Attribute] +compareNodeFormatWithoutRank cList (_,l) = [GV.textLabel (compareScientificName l), GV.style GV.wedged, GVA.Color (selectColors (inTree l) cList)] + +-- | Colors from color list are selected according to in which of the compared trees the node is contained. +selectColors :: [Int] -> [GVA.Color] -> GVAC.ColorList +selectColors inTrees currentColorList = GVAC.toColorList (map (\i -> currentColorList !! i) inTrees) + +-- | A color list is sampled from the spectrum according to how many trees are compared. +makeColorList :: Int -> [GVA.Color] +makeColorList treeNumber = cList + where cList = map (\i -> GVAC.HSV ((fromIntegral i/fromIntegral neededColors) * 0.708) 0.5 1.0) [0..neededColors] + neededColors = treeNumber - 1 + +-- | Write tree representation either as dot or json to provided file path +writeTree :: String -> String -> Bool -> Gr SimpleTaxon Double -> IO () +writeTree requestedFormat outputDirectoryPath withRank inputGraph = do + case requestedFormat of + "dot" -> writeDotTree outputDirectoryPath withRank inputGraph + "json"-> writeJsonTree outputDirectoryPath inputGraph + _ -> writeDotTree outputDirectoryPath withRank inputGraph + +-- | Write tree representation as dot to provided file path. +-- Graphviz tools like dot can be applied to the written .dot file to generate e.g. svg-format images. +writeDotTree :: String -> Bool -> Gr SimpleTaxon Double -> IO () +writeDotTree outputDirectoryPath withRank inputGraph = do + let diagram = drawTaxonomy withRank (grev inputGraph) + writeFile (outputDirectoryPath ++ "taxonomy.dot") diagram + +-- | Write tree representation as json to provided file path. +-- You can visualize the result for example with 3Djs. +writeJsonTree :: String -> Gr SimpleTaxon Double -> IO () +writeJsonTree outputDirectoryPath inputGraph = do + let jsonOutput = AE.encode (grev inputGraph) + L.writeFile (outputDirectoryPath ++ "taxonomy.json") jsonOutput diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,17 @@ +-*-change-log-*- + +2.0.0 [Florian Eggenhofer](mailto:egg@cs.uni-freiburg.de) 18. November 2019 + + * Changed to Biobase repository layout + * Added new Lineage datatype + * Using bytestring instead of string + +1.0.3 [Florian Eggenhofer](mailto:egg@cs.uni-freiburg.de) 10. August 2017 + + * Compatibility with newest versions of fgl + * Improved travis file + +1.0.2 [Florian Eggenhofer](mailto:egg@cs.uni-freiburg.de) 23. October 2016 + + * Changed datastructures for cutting the tree from list to vectors + * Improved travis file diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,15 @@ +Taxonomy [![Hackage](https://img.shields.io/hackage/v/Taxonomy.svg)](https://hackage.haskell.org/package/Taxonomy) [![Build Status](https://travis-ci.org/eggzilla/Taxonomy.svg?branch=master)](https://travis-ci.org/eggzilla/Taxonomy) +============= + +Haskell cabal Taxonomy libary contains tools, parsers, datastructures and visualisation +for the NCBI (National Center for Biotechnology Information) Taxonomy datasources. + +It can utilize information from the Entrez REST interface (via [EntrezHTTP](https://github.com/eggzilla/EntrezHTTP), +as well as from the files of the Taxonomy database dump. + +Input data is parsed into a FGL based datastructure, which enables a wealth of processing +steps like node distances, retrieval of parent nodes or extraction of +subtrees. + +Trees can be visualised via dot-format (graphviz) or +via json-format (d3js). diff --git a/Taxonomy.cabal b/Taxonomy.cabal --- a/Taxonomy.cabal +++ b/Taxonomy.cabal @@ -1,52 +1,53 @@ name: Taxonomy -version: 1.0.3 +version: 2.0.0 synopsis: Libary for parsing, processing and vizualization of taxonomy data description: Haskell cabal Taxonomy libary contains tools, parsers, datastructures and visualisation for the NCBI (National Center for Biotechnology Information) Taxonomy datasources. . It can utilize information from the REST interface via , - as well as from the files of the Taxonomy database . + as well as from the files of the Taxonomy database . . Input data is parsed into a FGL based datastructure, which enables a wealth of processing - steps like node distances, retrieval of parent nodes or extraction of - subtrees. + steps like node distances, retrieval of parent nodes or extraction of + subtrees. . Trees can be visualised via dot-format () . - <> + <> . or via json-format (). - . - The package contains tools based on this package. - + . + The package contains tools based on this package. + license: GPL-3 license-file: LICENSE author: Florian Eggenhofer -maintainer: florian.eggenhofer@univie.ac.at --- copyright: +maintainer: egg@informatik.uni-freiburg.de +-- copyright: category: Bioinformatics build-type: Simple -cabal-version: >=1.8 -Tested-With: GHC == 7.6.3, GHC == 7.8.4, GHC == 7.10.3, GHC == 8.0.2, GHC == 8.2.1 +cabal-version: >=1.10.0 +Tested-With: GHC == 8.4.4, GHC == 8.6.5 +Extra-Source-Files: + README.md ChangeLog.md + + source-repository head type: git location: https://github.com/eggzilla/Taxonomy source-repository this type: git - location: https://github.com/eggzilla/Taxonomy/tree/v1.0.3 - tag: v1.0.3 + location: https://github.com/eggzilla/Taxonomy/tree/v2.0.0 + tag: v2.0.0 -library - -- Modules exported by the library. - exposed-modules: Bio.TaxonomyData, Bio.Taxonomy - - -- compiler-options: +Library ghc-options: -Wall -fno-warn-unused-do-bind - - -- Other library packages from which modules are imported. + default-language: Haskell2010 build-depends: base >=4.5 && <5, parsec, either-unwrap, fgl>=5.5.4.0, text, graphviz>=2999.12.0.4, bytestring, aeson, vector - - -- Directories containing source files. - hs-source-dirs: src - + Hs-source-dirs: . + Exposed-modules: Biobase.Taxonomy.Types + Biobase.Taxonomy.Import + Biobase.Taxonomy.Visualization + Biobase.Taxonomy.Utils + Biobase.Taxonomy diff --git a/src/Bio/Taxonomy.hs b/src/Bio/Taxonomy.hs deleted file mode 100644 --- a/src/Bio/Taxonomy.hs +++ /dev/null @@ -1,635 +0,0 @@ --- | Functions for parsing, processing and visualization of taxonomy data. --- --- === Usage example: --- * Read in taxonomy data --- --- > eitherTaxtree <- readNamedTaxonomy "/path/to/NCBI_taxonomydump_directory" --- --- * Process data --- --- > let subtree = extractTaxonomySubTreebyLevel [562] (fromRight eitherTaxTree) (Just 4) --- --- * Visualize result --- --- tput "/path/to/dotdirectory" subtree -module Bio.Taxonomy ( -- * Datatypes - -- Datatypes used to represent taxonomy data - module Bio.TaxonomyData, - -- * Parsing - -- Functions prefixed with "read" read from filepaths, functions with parse from Haskell Strings. - readTaxonomy, - readNamedTaxonomy, - parseTaxonomy, - parseNCBITaxCitations, - readNCBITaxCitations, - parseNCBITaxDelNodes, - readNCBITaxDelNodes, - parseNCBITaxDivisions, - readNCBITaxDivisions, - parseNCBITaxGenCodes, - readNCBITaxGenCodes, - parseNCBITaxMergedNodes, - readNCBITaxMergedNodes, - parseNCBITaxNames, - readNCBITaxNames, - parseNCBITaxNodes, - readNCBITaxNodes, - parseNCBISimpleTaxons, - readNCBISimpleTaxons, - readNCBITaxonomyDatabase, - -- * Processing - compareSubTrees, - extractTaxonomySubTreebyLevel, - extractTaxonomySubTreebyLevelNew, - extractTaxonomySubTreebyRank, - safeNodePath, - getParentbyRank, - -- * Visualization - drawTaxonomyComparison, - drawTaxonomy, - writeTree, - writeDotTree, - writeJsonTree - ) where -import Prelude -import System.IO -import Bio.TaxonomyData -import Text.Parsec.Prim (runP) -import Text.ParserCombinators.Parsec -import Control.Monad -import Data.List -import qualified Data.Vector as V -import Data.Maybe -import qualified Data.Either.Unwrap as E -import Data.Graph.Inductive.Graph -import Data.Graph.Inductive.Query.SP (sp) -import Data.Graph.Inductive.Query.BFS (level) -import Data.Graph.Inductive.Tree -import Data.Graph.Inductive.Basic -import qualified Data.GraphViz as GV -import qualified Data.GraphViz.Printing as GVP -import qualified Data.GraphViz.Attributes.Colors as GVAC -import qualified Data.GraphViz.Attributes.Complete as GVA -import qualified Data.ByteString.Lazy.Char8 as L -import qualified Data.Aeson as AE -import qualified Data.Text.Lazy as T --------------------------------------------------------- - ---------------------------------------- --- Parsing functions - --- | NCBI taxonomy dump nodes and names in the input directory path are parsed and a SimpleTaxon tree with scientific names for each node is generated. -readNamedTaxonomy :: String -> IO (Either ParseError (Gr SimpleTaxon Double)) -readNamedTaxonomy directoryPath = do - nodeNames <- readNCBITaxNames (directoryPath ++ "names.dmp") - if E.isLeft nodeNames - then return (Left (E.fromLeft nodeNames)) - else do - let rightNodeNames = V.fromList (E.fromRight nodeNames) - let filteredNodeNames = V.filter isScientificName rightNodeNames - parseFromFileEncISO88591 (genParserNamedTaxonomyGraph filteredNodeNames) (directoryPath ++ "nodes.dmp") - -isScientificName :: TaxName -> Bool -isScientificName name = nameClass name == scientificNameT - where scientificNameT = T.pack "scientific name" - --- | NCBI taxonomy dump nodes and names in the input directory path are parsed and a SimpleTaxon tree is generated. -readTaxonomy :: String -> IO (Either ParseError (Gr SimpleTaxon Double)) -readTaxonomy = parseFromFileEncISO88591 genParserTaxonomyGraph - --- | NCBI taxonomy dump nodes and names in the input directory path are parsed and a SimpleTaxon tree is generated. -parseTaxonomy :: String -> Either ParseError (Gr SimpleTaxon Double) -parseTaxonomy = parse genParserTaxonomyGraph "parseTaxonomy" - -genParserTaxonomyGraph :: GenParser Char st (Gr SimpleTaxon Double) -genParserTaxonomyGraph = do - nodesEdges <- many1 (try genParserGraphNodeEdge) - optional eof - let (nodesList,edgesList) = unzip nodesEdges - --let taxedges = filter (\(a,b,_) -> a /= b) edgesList - let taxedges = filter notLoopEdge edgesList - --let taxnodes = concat nodesList - --return (mkGraph taxnodes taxedges) - return $! mkGraph nodesList taxedges - - -notLoopEdge :: (Int,Int,a) -> Bool -notLoopEdge (a,b,_) = a /= b - -genParserNamedTaxonomyGraph :: V.Vector TaxName -> GenParser Char st (Gr SimpleTaxon Double) -genParserNamedTaxonomyGraph filteredNodeNames = do - nodesEdges <- (many1 (try genParserGraphNodeEdge)) - optional eof - let (nodesList,edgesList) = V.unzip (V.fromList nodesEdges) - let taxedges = V.filter notLoopEdge edgesList - let taxnamednodes = V.map (setNodeScientificName filteredNodeNames) nodesList - return $! mkGraph (V.toList taxnamednodes) (V.toList taxedges) - -setNodeScientificName :: V.Vector TaxName -> (t, SimpleTaxon) -> (t, SimpleTaxon) -setNodeScientificName inputTaxNames (inputNode,inputTaxon) = outputNode - where maybeRetrievedName = V.find (isTaxNameIdSimpleTaxid inputTaxon) inputTaxNames - retrievedName = maybe (T.pack "no name") nameTxt maybeRetrievedName - outputNode = (inputNode,inputTaxon{simpleScientificName = retrievedName}) - -isTaxNameIdSimpleTaxid :: SimpleTaxon -> TaxName -> Bool -isTaxNameIdSimpleTaxid inputTaxon inputTaxName = nameTaxId inputTaxName == simpleTaxId inputTaxon - - -genParserGraphNodeEdge :: GenParser Char st ((Int,SimpleTaxon),(Int,Int,Double)) -genParserGraphNodeEdge = do - _simpleTaxId <- many1 digit - string "\t|\t" - _simpleParentTaxId <- many1 digit - string "\t|\t" - _simpleRank <- many1 (noneOf "\t") - many1 (noneOf "\n") - char '\n' - let _simpleTaxIdInt = readInt _simpleTaxId - let _simpleParentTaxIdInt = readInt _simpleParentTaxId - return ((_simpleTaxIdInt,SimpleTaxon _simpleTaxIdInt T.empty _simpleParentTaxIdInt (readRank _simpleRank)),(_simpleTaxIdInt,_simpleParentTaxIdInt,1 :: Double)) - --- | parse NCBITaxCitations from input string -parseNCBITaxCitations :: String -> Either ParseError [TaxCitation] -parseNCBITaxCitations = parse genParserNCBITaxCitations "parseTaxCitations" - --- | parse NCBITaxCitations from input filePath -readNCBITaxCitations :: String -> IO (Either ParseError [TaxCitation]) -readNCBITaxCitations = parseFromFileEncISO88591 genParserNCBITaxCitations - --- | parse NCBITaxDelNodes from input string -parseNCBITaxDelNodes :: String -> Either ParseError [TaxDelNode] -parseNCBITaxDelNodes = parse genParserNCBITaxDelNodes "parseTaxDelNodes" - --- | parse NCBITaxDelNodes from input filePath -readNCBITaxDelNodes :: String -> IO (Either ParseError [TaxDelNode]) -readNCBITaxDelNodes = parseFromFile genParserNCBITaxDelNodes - --- | parse NCBITaxDivisons from input string -parseNCBITaxDivisions :: String -> Either ParseError [TaxDivision] -parseNCBITaxDivisions = parse genParserNCBITaxDivisons "parseTaxDivisons" - --- | parse NCBITaxDivisons from input filePath -readNCBITaxDivisions :: String -> IO (Either ParseError [TaxDivision]) -readNCBITaxDivisions = parseFromFile genParserNCBITaxDivisons - --- | parse NCBITaxGenCodes from input string -parseNCBITaxGenCodes :: String -> Either ParseError [TaxGenCode] -parseNCBITaxGenCodes = parse genParserNCBITaxGenCodes "parseTaxGenCodes" - --- | parse NCBITaxGenCodes from input filePath -readNCBITaxGenCodes :: String -> IO (Either ParseError [TaxGenCode]) -readNCBITaxGenCodes = parseFromFile genParserNCBITaxGenCodes - --- | parse NCBITaxMergedNodes from input string -parseNCBITaxMergedNodes :: String -> Either ParseError [TaxMergedNode] -parseNCBITaxMergedNodes = parse genParserNCBITaxMergedNodes "parseTaxMergedNodes" - --- | parse NCBITaxMergedNodes from input filePath -readNCBITaxMergedNodes :: String -> IO (Either ParseError [TaxMergedNode]) -readNCBITaxMergedNodes = parseFromFile genParserNCBITaxMergedNodes - --- | parse NCBITaxNames from input string -parseNCBITaxNames :: String -> Either ParseError [TaxName] -parseNCBITaxNames = parse genParserNCBITaxNames "parseTaxNames" - --- | parse NCBITaxNames from input filePath -readNCBITaxNames :: String -> IO (Either ParseError [TaxName]) -readNCBITaxNames = parseFromFile genParserNCBITaxNames - --- | parse NCBITaxNames from input string -parseNCBITaxNodes :: String -> Either ParseError TaxNode -parseNCBITaxNodes = parse genParserNCBITaxNode "parseTaxNode" - --- | parse NCBITaxCitations from input filePath -readNCBITaxNodes :: String -> IO (Either ParseError [TaxNode]) -readNCBITaxNodes = parseFromFile genParserNCBITaxNodes - --- | parse NCBISimpleTaxNames from input string -parseNCBISimpleTaxons :: String -> Either ParseError SimpleTaxon -parseNCBISimpleTaxons = parse genParserNCBISimpleTaxon "parseSimpleTaxon" - --- | parse NCBITaxCitations from input filePath -readNCBISimpleTaxons :: String -> IO (Either ParseError [SimpleTaxon]) -readNCBISimpleTaxons = parseFromFile genParserNCBISimpleTaxons - --- | Parse the input as NCBITax datatype -readNCBITaxonomyDatabase :: String -> IO (Either [String] NCBITaxDump) -readNCBITaxonomyDatabase folder = do - citations <- readNCBITaxCitations (folder ++ "citations.dmp") - let citationsError = extractParseError citations - taxdelNodes <- readNCBITaxDelNodes (folder ++ "delnodes.dmp") - let delNodesError = extractParseError taxdelNodes - divisons <- readNCBITaxDivisions (folder ++ "division.dmp") - let divisonsError = extractParseError divisons - genCodes <- readNCBITaxGenCodes (folder ++ "gencode.dmp") - let genCodesError = extractParseError genCodes - mergedNodes <- readNCBITaxMergedNodes (folder ++ "merged.dmp") - let mergedNodesError = extractParseError mergedNodes - names <- readNCBITaxNames (folder ++ "names.dmp") - let namesError = extractParseError names - taxnodes <- readNCBITaxNodes (folder ++ "nodes.dmp") - let nodesError = extractParseError taxnodes - let parseErrors = [citationsError, delNodesError, divisonsError, genCodesError, mergedNodesError, namesError, nodesError] - return (checkParsing parseErrors citations taxdelNodes divisons genCodes mergedNodes names taxnodes) - -genParserNCBITaxCitations :: GenParser Char st [TaxCitation] -genParserNCBITaxCitations = many1 genParserNCBITaxCitation - -genParserNCBITaxDelNodes :: GenParser Char st [TaxDelNode] -genParserNCBITaxDelNodes = many1 genParserNCBITaxDelNode - -genParserNCBITaxDivisons :: GenParser Char st [TaxDivision] -genParserNCBITaxDivisons = many1 genParserNCBITaxDivision - -genParserNCBITaxGenCodes :: GenParser Char st [TaxGenCode] -genParserNCBITaxGenCodes = many1 genParserNCBITaxGenCode - - -genParserNCBITaxMergedNodes :: GenParser Char st [TaxMergedNode] -genParserNCBITaxMergedNodes = many1 genParserNCBITaxMergedNode - - -genParserNCBITaxNames :: GenParser Char st [TaxName] -genParserNCBITaxNames = many1 genParserNCBITaxName - -genParserNCBITaxNodes :: GenParser Char st [TaxNode] -genParserNCBITaxNodes = many1 genParserNCBITaxNode - -genParserNCBISimpleTaxons :: GenParser Char st [SimpleTaxon] -genParserNCBISimpleTaxons = many1 genParserNCBISimpleTaxon - - -genParserNCBITaxCitation :: GenParser Char st TaxCitation -genParserNCBITaxCitation = do - _citId <- many1 digit - string "\t|\t" - _citKey <- optionMaybe (many1 (noneOf "\t")) - string "\t|\t" - _pubmedId <- optionMaybe (many1 digit) - string "\t|\t" - _medlineId <- optionMaybe (many1 digit) - tab - char '|' - _url <- genParserTaxURL - char '|' - tab - _text <- optionMaybe (many1 (noneOf "\t")) - string "\t|\t" - _taxIdList <- optionMaybe (many1 genParserTaxIdList) - string "\t|\n" - return $ TaxCitation (readInt _citId) _citKey (liftM readInt _pubmedId) (liftM readInt _medlineId) _url _text _taxIdList - -genParserNCBITaxDelNode :: GenParser Char st TaxDelNode -genParserNCBITaxDelNode = do - taxdelNode <- many1 digit - space - char '|' - char '\n' - return $ TaxDelNode (readInt taxdelNode) - -genParserNCBITaxDivision :: GenParser Char st TaxDivision -genParserNCBITaxDivision = do - _divisionId <- many1 digit - string "\t|\t" - _divisionCDE <- many1 upper - string "\t|\t" - _divisionName <- many1 (noneOf "\t") - string "\t|\t" - _comments <- optionMaybe (many1 (noneOf "\t")) - string "\t|\n" - return $ TaxDivision (readInt _divisionId) _divisionCDE _divisionName _comments - -genParserNCBITaxGenCode :: GenParser Char st TaxGenCode -genParserNCBITaxGenCode = do - _geneticCodeId <- many1 digit - string "\t|\t" - _abbreviation <- optionMaybe (many1 (noneOf "\t")) - string "\t|\t" - _genCodeName <- many1 (noneOf "\t") - string "\t|\t" - _cde <- many1 (noneOf "\t") - string "\t|\t" - _starts <- many1 (noneOf "\t") - string "\t|\n" - return $ TaxGenCode (readInt _geneticCodeId) _abbreviation _genCodeName _cde _starts - -genParserNCBITaxMergedNode :: GenParser Char st TaxMergedNode -genParserNCBITaxMergedNode = do - _oldTaxId <- many1 digit - string "\t|\t" - _newTaxId <- many1 digit - string "\t|\n" - return $ TaxMergedNode (readInt _oldTaxId) (readInt _newTaxId) - -genParserNCBITaxName :: GenParser Char st TaxName -genParserNCBITaxName = do - _taxId <- many1 digit - string "\t|\t" - _nameTxt <- many1 (noneOf "\t\n") - string "\t|\t" - _uniqueName <- optionMaybe (many1 (noneOf "\t\n")) - string "\t|\t" - _nameClass <- many1 (noneOf "\t\n") - tab - char '|' - newline - return $! TaxName (readInt _taxId) (T.pack _nameTxt) (maybe T.empty T.pack _uniqueName) (T.pack _nameClass) - -genParserNCBISimpleTaxon :: GenParser Char st SimpleTaxon -genParserNCBISimpleTaxon = do - _simpleTaxId <- many1 digit - string "\t|\t" - _simpleParentTaxId <- many1 digit - string "\t|\t" - _simpleRank <- many1 (noneOf "\t") - many1 (noneOf "\n") - char '\n' - return $! SimpleTaxon (readInt _simpleTaxId) T.empty (readInt _simpleParentTaxId) (readRank _simpleRank) - -genParserNCBITaxNode :: GenParser Char st TaxNode -genParserNCBITaxNode = do - _taxId <- many1 digit - string "\t|\t" - _parentTaxId <- many1 digit - string "\t|\t" - _rank <- many1 (noneOf "\t") - string "\t|\t" - _emblCode <- optionMaybe (many1 (noneOf "\t")) - string "\t|\t" - _divisionId <- many1 digit - string "\t|\t" - _inheritedDivFlag <- many1 digit - string "\t|\t" - _geneticCodeId <- many1 digit - string "\t|\t" - _inheritedGCFlag <- many1 digit - string "\t|\t" - _mitochondrialGeneticCodeId <- many1 digit - string "\t|\t" - _inheritedMGCFlag <- many1 digit - string "\t|\t" - _genBankHiddenFlag <- many1 digit - string "\t|\t" - _hiddenSubtreeRootFlag <- many1 digit - string "\t|\t" - _comments <- optionMaybe (many1 (noneOf "\t")) - tab - char '|' - char '\n' - return $ TaxNode (readInt _taxId) (readInt _parentTaxId) (readRank _rank) _emblCode _divisionId (readBool _inheritedDivFlag) _geneticCodeId (readBool _inheritedGCFlag) _mitochondrialGeneticCodeId (readBool _inheritedMGCFlag) (readBool _genBankHiddenFlag) (readBool _hiddenSubtreeRootFlag) _comments - ---------------------------------------- --- Processing functions - --- | Extract a subtree correpsonding to input node paths to root. Only nodes in level number distance to root are included. Used in Ids2TreeCompare tool. -compareSubTrees :: [Gr SimpleTaxon Double] -> (Int,Gr CompareTaxon Double) -compareSubTrees graphs = (length graphs,resultGraph) - where treesLabNodes = map labNodes graphs - treesLabEdges = map labEdges graphs - mergedNodes = nub (concat treesLabNodes) - mergedEdges = nub (concat treesLabEdges) - --annotate node in which of the compared trees they are present - comparedNodes = annotateTaxonsDifference treesLabNodes mergedNodes - resultGraph = mkGraph comparedNodes mergedEdges :: Gr CompareTaxon Double - -annotateTaxonsDifference :: [[LNode SimpleTaxon]] -> [LNode SimpleTaxon] -> [LNode CompareTaxon] -annotateTaxonsDifference treesNodes mergedtreeNodes = comparedNodes - where comparedNodes = map (annotateTaxonDifference indexedTreesNodes) mergedtreeNodes - indexedTreesNodes = zip [0..(length treesNodes)] treesNodes - - -annotateTaxonDifference :: [(Int,[LNode SimpleTaxon])] -> LNode SimpleTaxon -> LNode CompareTaxon -annotateTaxonDifference indexedTreesNodes mergedtreeNode = comparedNode - where comparedNode = (simpleTaxId (snd mergedtreeNode),CompareTaxon (simpleScientificName (snd mergedtreeNode)) (simpleRank (snd mergedtreeNode)) currentInTree) - currentInTree = concatMap (\(i,treeNodes) -> [i | mergedtreeNode `elem` treeNodes]) indexedTreesNodes - --- | Extract a subtree corresponding to input node paths to root. Only nodes in level number distance to root are included. Used in Ids2Tree tool. -extractTaxonomySubTreebyLevel :: [Node] -> Gr SimpleTaxon Double -> Maybe Int -> Gr SimpleTaxon Double -extractTaxonomySubTreebyLevel inputNodes graph levelNumber = taxonomySubTree - where paths = nub (concatMap (getPath (1 :: Node) graph) inputNodes) - contexts = map (context graph) paths - lnodes = map labNode' contexts - ledges = nub (concatMap (out graph . fst) lnodes) - unfilteredTaxonomySubTree = mkGraph lnodes ledges :: Gr SimpleTaxon Double - filteredLNodes = filterNodesByLevel levelNumber lnodes unfilteredTaxonomySubTree - filteredledges = nub (concatMap (out graph . fst) filteredLNodes) - taxonomySubTree = mkGraph filteredLNodes filteredledges :: Gr SimpleTaxon Double - --- | Extract a subtree corresponding to input node paths to root. Only nodes in level number distance to root are included. Used in Ids2Tree tool. -extractTaxonomySubTreebyLevelNew :: [Node] -> Gr SimpleTaxon Double -> Maybe Int -> Gr SimpleTaxon Double -extractTaxonomySubTreebyLevelNew inputNodes graph levelNumber = taxonomySubTree - where inputNodeVector = V.fromList inputNodes - paths = V.concatMap (getVectorPath (1 :: Node) graph) inputNodeVector - contexts = V.map (context graph) paths - vlnodes = V.map labNode' contexts - ledges = concatMap (out graph . fst) lnodes - lnodes = V.toList vlnodes - --ledges = V.toList vledges - unfilteredTaxonomySubTree = mkGraph lnodes ledges :: Gr SimpleTaxon Double - filteredLNodes = filterNodesByLevel levelNumber lnodes unfilteredTaxonomySubTree - --filteredLNodesVector = V.fromList filteredLNodes - filteredledges = concatMap (out graph . fst) filteredLNodes - --filteredledges = V.toList filteredledgesVector - taxonomySubTree = mkGraph filteredLNodes filteredledges :: Gr SimpleTaxon Double - --- | Extract a subtree corresponding to input node paths to root. If a Rank is provided, all node that are less or equal are omitted -extractTaxonomySubTreebyRank :: [Node] -> Gr SimpleTaxon Double -> Maybe Rank -> Gr SimpleTaxon Double -extractTaxonomySubTreebyRank inputNodes graph highestRank = taxonomySubTree - where paths = nub (concatMap (getPath (1 :: Node) graph) inputNodes) - contexts = map (context graph) paths - lnodes = map labNode' contexts - filteredLNodes = filterNodesByRank highestRank lnodes - filteredledges = nub (concatMap (out graph . fst) filteredLNodes) - taxonomySubTree = mkGraph filteredLNodes filteredledges :: Gr SimpleTaxon Double - -getVectorPath :: Node -> Gr SimpleTaxon Double -> Node -> V.Vector Node -getVectorPath root graph node = maybe V.empty V.fromList (sp node root graph) - -getPath :: Node -> Gr SimpleTaxon Double -> Node -> Path -getPath root graph node = maybe [] id (sp node root graph) - --- | Extract parent node with specified Rank -getParentbyRank :: Node -> Gr SimpleTaxon Double -> Maybe Rank -> Maybe (Node, SimpleTaxon) -getParentbyRank inputNode graph requestedRank = filteredLNode - where path = maybe [] id (sp (inputNode :: Node) (1 :: Node) graph) - nodeContext = map (context graph) path - lnode = map labNode' nodeContext - filteredLNode = findNodeByRank requestedRank lnode - --- | Filter nodes by distance from root -filterNodesByLevel :: Maybe Int -> [(Node, SimpleTaxon)] -> Gr SimpleTaxon Double -> [(Node, SimpleTaxon)] -filterNodesByLevel levelNumber inputNodes graph - | isJust levelNumber = filteredNodes - | otherwise = inputNodes - --distances of all nodes to root - where nodedistances = level (1::Node) (undir graph) - sortedNodeDistances = sortBy sortByNodeID nodedistances - sortedInputNodes = sortBy sortByNodeID inputNodes - zippedNodeDistancesInputNodes = zip sortedNodeDistances sortedInputNodes - zippedFilteredNodes = filter (\((_,distance),(_,_)) -> distance <= fromJust levelNumber) zippedNodeDistancesInputNodes - filteredNodes = map snd zippedFilteredNodes - -sortByNodeID :: (Node,a) -> (Node,a) -> Ordering -sortByNodeID (n1, _) (n2, _) - | n1 < n2 = GT - | n1 > n2 = LT - | n1 == n2 = EQ - | otherwise = EQ - --- | Find only taxons of a specific rank in a list of input taxons -findNodeByRank :: Maybe Rank -> [(t, SimpleTaxon)] -> Maybe (t, SimpleTaxon) -findNodeByRank requestedRank inputNodes - | isJust requestedRank = filteredNodes - | otherwise = Nothing - where filteredNodes = find (\(_,t) -> simpleRank t == fromJust requestedRank) inputNodes - --- | Filter a list of input taxons for a minimal provided rank -filterNodesByRank :: Maybe Rank -> [(t, SimpleTaxon)] -> [(t, SimpleTaxon)] -filterNodesByRank highestRank inputNodes - | isJust highestRank = filteredNodes - | otherwise = inputNodes - where filteredNodes = filter (\(_,t) -> simpleRank t >= fromJust highestRank) inputNodes ++ noRankNodes - noRankNodes = filter (\(_,t) -> simpleRank t == Norank) inputNodes - --- | Returns path between 2 maybe nodes. Used in TreeDistance tool. -safeNodePath :: Maybe Node -> Gr SimpleTaxon Double -> Maybe Node -> Either String Path -safeNodePath nodeid1 graphOutput nodeid2 - | isJust nodeid1 && isJust nodeid2 = Right (maybe [] id (sp (fromJust nodeid1) (fromJust nodeid2) (undir graphOutput))) - | otherwise = Left "Both taxonomy ids must be provided for distance computation" - ---------------------------------------- --- Visualisation functions - --- | Draw graph in dot format. Used in Ids2Tree tool. -drawTaxonomy :: Bool -> Gr SimpleTaxon Double -> String -drawTaxonomy withRank inputGraph = do - let nodeFormating = if withRank then nodeFormatWithRank else nodeFormatWithoutRank - let params = GV.nonClusteredParams {GV.isDirected = True - , GV.globalAttributes = [GV.GraphAttrs [GVA.Size (GVA.GSize (20 :: Double) (Just (20 :: Double)) False)]] - , GV.isDotCluster = const True - --, GV.fmtNode = \ (_,l) -> [GV.textLabel (TL.pack (show (simpleRank l) ++ "\n" ++ T.unpack (simpleScientificName l)))] - , GV.fmtNode = nodeFormating - , GV.fmtEdge = const [] - } - let dotFormat = GV.graphToDot params inputGraph - let dottext = GVP.renderDot $ GVP.toDot dotFormat - T.unpack dottext - -nodeFormatWithRank :: (t, SimpleTaxon) -> [GVA.Attribute] -nodeFormatWithRank (_,l) = [GV.textLabel (T.concat [T.pack (show (simpleRank l)), T.pack ("\n") , simpleScientificName l])] - -nodeFormatWithoutRank :: (t, SimpleTaxon) -> [GVA.Attribute] -nodeFormatWithoutRank (_,l) = [GV.textLabel (simpleScientificName l)] - --- | Draw tree comparison graph in dot format. Used in Ids2TreeCompare tool. -drawTaxonomyComparison :: Bool -> (Int,Gr CompareTaxon Double) -> String -drawTaxonomyComparison withRank (treeNumber,inputGraph) = do - let cList = makeColorList treeNumber - let nodeFormating = if withRank then (compareNodeFormatWithRank cList) else (compareNodeFormatWithoutRank cList) - let params = GV.nonClusteredParams {GV.isDirected = True - , GV.globalAttributes = [] - , GV.isDotCluster = const True - --, GV.fmtNode = \ (_,l) -> [GV.textLabel (TL.pack (show (compareRank l) ++ "\n" ++ B.unpack (compareScientificName l))), GV.style GV.wedged, GVA.Color (selectColors (inTree l) cList)] - , GV.fmtNode = nodeFormating - , GV.fmtEdge = const [] - } - let dotFormat = GV.graphToDot params (grev inputGraph) - let dottext = GVP.renderDot $ GVP.toDot dotFormat - T.unpack dottext - -compareNodeFormatWithRank :: [GVA.Color] -> (t, CompareTaxon) -> [GVA.Attribute] -compareNodeFormatWithRank cList (_,l) = [GV.textLabel (T.concat [T.pack (show (compareRank l) ++ "\n"),compareScientificName l]), GV.style GV.wedged, GVA.Color (selectColors (inTree l) cList)] - -compareNodeFormatWithoutRank :: [GVA.Color] -> (t, CompareTaxon) -> [GVA.Attribute] -compareNodeFormatWithoutRank cList (_,l) = [GV.textLabel (compareScientificName l), GV.style GV.wedged, GVA.Color (selectColors (inTree l) cList)] - --- | Colors from color list are selected according to in which of the compared trees the node is contained. -selectColors :: [Int] -> [GVA.Color] -> GVAC.ColorList -selectColors inTrees currentColorList = GVAC.toColorList (map (\i -> currentColorList !! i) inTrees) - --- | A color list is sampled from the spectrum according to how many trees are compared. -makeColorList :: Int -> [GVA.Color] -makeColorList treeNumber = cList - where cList = map (\i -> GVAC.HSV ((fromIntegral i/fromIntegral neededColors) * 0.708) 0.5 1.0) [0..neededColors] - neededColors = treeNumber - 1 - --- | Write tree representation either as dot or json to provided file path -writeTree :: String -> String -> Bool -> Gr SimpleTaxon Double -> IO () -writeTree requestedFormat outputDirectoryPath withRank inputGraph = do - case requestedFormat of - "dot" -> writeDotTree outputDirectoryPath withRank inputGraph - "json"-> writeJsonTree outputDirectoryPath inputGraph - _ -> writeDotTree outputDirectoryPath withRank inputGraph - --- | Write tree representation as dot to provided file path. --- Graphviz tools like dot can be applied to the written .dot file to generate e.g. svg-format images. -writeDotTree :: String -> Bool -> Gr SimpleTaxon Double -> IO () -writeDotTree outputDirectoryPath withRank inputGraph = do - let diagram = drawTaxonomy withRank (grev inputGraph) - writeFile (outputDirectoryPath ++ "taxonomy.dot") diagram - --- | Write tree representation as json to provided file path. --- You can visualize the result for example with 3Djs. -writeJsonTree :: String -> Gr SimpleTaxon Double -> IO () -writeJsonTree outputDirectoryPath inputGraph = do - let jsonOutput = AE.encode (grev inputGraph) - L.writeFile (outputDirectoryPath ++ "taxonomy.json") jsonOutput - ---------------------------------------- --- Auxiliary functions -readInt :: String -> Int -readInt = read - -readBool :: String -> Bool -readBool "0" = False -readBool "1" = True -readBool _ = False - -readRank :: String -> Rank -readRank a = read a :: Rank - -genParserTaxIdList :: GenParser Char st Int -genParserTaxIdList = do - optional (char ' ') - _taxId <- many1 digit - optional (char ' ') - return (readInt _taxId) - -genParserTaxURL :: GenParser Char st (Maybe String) -genParserTaxURL = do - tab - url1 <- optionMaybe (many1 (noneOf "\t")) - tab - url2 <- optionMaybe (many1 (noneOf "|")) - return (concatenateURLParts url1 url2) - -concatenateURLParts :: Maybe String -> Maybe String -> Maybe String -concatenateURLParts url1 url2 - | isJust url1 && isJust url2 = maybeStringConcat url1 url2 - | isJust url1 && isNothing url2 = url1 - | otherwise = Nothing - -maybeStringConcat :: Maybe String -> Maybe String -> Maybe String -maybeStringConcat = liftM2 (++) - -readEncodedFile :: TextEncoding -> FilePath -> IO String -readEncodedFile encoding name = do - handle <- openFile name ReadMode - hSetEncoding handle encoding - hGetContents handle - -parseFromFileEncISO88591 :: Parser a -> String -> IO (Either ParseError a) -parseFromFileEncISO88591 parser fname = do - input <- readEncodedFile latin1 fname - return (runP parser () fname input) - --- | check a list of parsing results for presence of Left aka Parse error -checkParsing :: [String] -> Either ParseError [TaxCitation] -> Either ParseError [TaxDelNode] -> Either ParseError [TaxDivision] -> Either ParseError [TaxGenCode] -> Either ParseError [TaxMergedNode] -> Either ParseError [TaxName] -> Either ParseError [TaxNode]-> Either [String] NCBITaxDump -checkParsing parseErrors citations taxdelNodes divisons genCodes mergedNodes names taxnodes - | join parseErrors == "" = Right (NCBITaxDump (E.fromRight citations) (E.fromRight taxdelNodes) (E.fromRight divisons) (E.fromRight genCodes) (E.fromRight mergedNodes) (E.fromRight names) (E.fromRight taxnodes)) - | otherwise = Left parseErrors - -extractParseError :: Either ParseError a -> String -extractParseError _parse - | E.isLeft _parse = show (E.fromLeft _parse) - | otherwise = "" diff --git a/src/Bio/TaxonomyData.hs b/src/Bio/TaxonomyData.hs deleted file mode 100644 --- a/src/Bio/TaxonomyData.hs +++ /dev/null @@ -1,292 +0,0 @@ --- | This module contains data structures for --- taxonomy data - -{-# LANGUAGE FlexibleInstances #-} - -module Bio.TaxonomyData where -import Prelude ---import qualified Data.ByteString as B -import qualified Data.Aeson as A -import qualified Data.Vector as V ---import Data.Graph.Inductive -import Data.Graph.Inductive.Graph -import Data.Graph.Inductive.Tree -import qualified Data.Text as T -import qualified Data.Text.Lazy as TL ---import qualified Data.Text.Encoding - --- | SimpleTaxon only contains the most relevant fields of a taxonomy entry. --- For all annotaded fields use the Taxon datatype and its associated functions -data SimpleTaxon = SimpleTaxon - { - -- node id in GenBank - simpleTaxId :: Int, - simpleScientificName :: TL.Text, - -- parent node id in GenBank taxonomy database - simpleParentTaxId :: Int, - -- rank of this node (superkingdom, kingdom, ...) - simpleRank :: Rank - } - deriving (Show, Read, Eq) - --- | Datastructure for tree comparisons -data CompareTaxon = CompareTaxon - { - compareScientificName :: TL.Text, - compareRank :: Rank, - -- number indicating in which trees, - inTree :: [Int] - } - deriving (Show, Read, Eq) - --- | Data structure for Entrez taxonomy fetch result -data Taxon = Taxon - { taxonTaxId :: Int - , taxonScientificName :: String - , taxonParentTaxId :: Int - , taxonRank :: Rank - , division :: String - , geneticCode :: TaxGenCode - , mitoGeneticCode :: TaxGenCode - , lineage :: String - , lineageEx :: [LineageTaxon] - , createDate :: String - , updateDate :: String - , pubDate :: String - } deriving (Show, Eq) - - -data TaxonName = TaxonName - { classCDE :: String - , dispName :: String - } deriving (Show, Eq) - --- | Lineage Taxons denote all parent Taxonomy nodes of a node retrieved by Entrez fetch -data LineageTaxon = LineageTaxon - { lineageTaxId :: Int - , lineageScienticName :: String - , lineageRank :: Rank} - deriving (Show, Eq) - --- | NCBI Taxonomy database dump hierachichal data structure --- as defined in ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump_readme.txt -data NCBITaxDump = NCBITaxDump - { - taxCitations :: [TaxCitation], - taxDelNodes :: [TaxDelNode], - taxDivisions :: [TaxDivision], - taxGenCodes :: [TaxGenCode], - taxMergedNodes :: [TaxMergedNode], - taxNames :: [TaxName], - taxNodes :: [TaxNode] - } - deriving (Show, Read, Eq) - --- | Datastructure for entries of Taxonomy database dump citations file -data TaxCitation = TaxCitation - { - -- the unique id of citation - citId :: Int, - -- citation key - citKey :: Maybe String, - -- unique id in PubMed database (0 if not in PubMed) - pubmedId :: Maybe Int, - -- unique id in MedLine database (0 if not in MedLine) - medlineId :: Maybe Int, - -- URL associated with citation - url :: Maybe String, - -- any text (usually article name and authors) - -- The following characters are escaped in this text by a backslash: - -- newline (appear as "\n"), - -- tab character ("\t"), - -- double quotes ('\"'), - -- backslash character ("\\"). - text :: Maybe String, - -- list of node ids separated by a single space - taxIdList :: Maybe [Int] - } - deriving (Show, Read, Eq) - --- | Datastructure for entries of Taxonomy database dump deleted nodes file -data TaxDelNode = TaxDelNode - { - -- deleted node id - delTaxId :: Int - } - deriving (Show, Read, Eq) - --- | Datastructure for entries of Taxonomy database dump division file -data TaxDivision = TaxDivision - { - -- taxonomy database division id - divisionId :: Int, - -- GenBank division code (three characters) - divisionCDE :: String, - -- e.g. BCT, PLN, VRT, MAM, PRI... - divisonName :: String, - divisionComments :: Maybe String - } - deriving (Show, Read, Eq) - --- | Datastructure for entries of Taxonomy database dump gencode file -data TaxGenCode = TaxGenCode - { - -- GenBank genetic code id - geneticCodeId :: Int, - -- genetic code name abbreviation - abbreviation :: Maybe String, - -- genetic code name - geneCodeName :: String, - -- translation table for this genetic code - cde :: String, - -- start codons for this genetic code - starts :: String - } - deriving (Show, Read, Eq) - --- | Datastructure for entries of Taxonomy database dump mergednodes file -data TaxMergedNode = TaxMergedNode - { - -- id of nodes which has been merged - oldTaxId :: Int, - -- id of nodes which is result of merging - newTaxId :: Int - } - deriving (Show, Read, Eq) - --- | Datastructure for entries of Taxonomy database dump names file -data TaxName = TaxName - { - -- the id of node associated with this name - nameTaxId :: Int, - -- name itself - nameTxt :: TL.Text, - -- the unique variant of this name if name not unique - uniqueName :: TL.Text, - -- (synonym, common name, ...) - nameClass :: TL.Text - } - deriving (Show, Read, Eq) - --- | Taxonomic ranks: NCBI uses the uncommon Speciessubgroup -data Rank = Norank | Form | Variety | Infraspecies | Subspecies | Speciessubgroup | Species | Speciesgroup | Superspecies | Series | Section | Subgenus | Genus | Subtribe | Tribe | Supertribe | Subfamily | Family | Superfamily | Parvorder | Infraorder | Suborder | Order | Superorder | Magnorder | Cohort | Legion | Parvclass | Infraclass | Subclass | Class | Superclass | Microphylum | Infraphylum | Subphylum | Phylum | Superphylum | Infrakingdom | Subkingdom | Kingdom | Superkingdom | Domain deriving (Eq, Ord, Show, Bounded, Enum) - -readsRank :: String -> [(Rank, String)] -instance Read Rank where - readsPrec _ = readsRank - -readsRank input -- = [(Domain x)| x <- reads input ] - | input == "domain" = [(Domain,"")] - | input == "superkingdom" = [(Superkingdom,"")] - | input == "kingdom" = [(Kingdom,"")] - | input == "subkingdom" = [(Subkingdom,"")] - | input == "infrakingdom" = [(Infrakingdom,"")] - | input == "superphylum" = [(Superphylum,"")] - | input == "phylum" = [(Phylum,"")] - | input == "subphylum" = [(Subphylum,"")] - | input == "infraphylum" = [(Infraphylum,"")] - | input == "microphylum" = [(Microphylum,"")] - | input == "superclass" = [(Superclass,"")] - | input == "class" = [(Class,"")] - | input == "subclass" = [(Subclass,"")] - | input == "infraclass" = [(Infraclass,"")] - | input == "parvclass " = [(Parvclass ,"")] - | input == "legion" = [(Legion,"")] - | input == "cohort" = [(Cohort,"")] - | input == "magnorder " = [(Magnorder ,"")] - | input == "superorder" = [(Superorder,"")] - | input == "order" = [(Order,"")] - | input == "suborder" = [(Suborder,"")] - | input == "infraorder" = [(Infraorder,"")] - | input == "parvorder" = [(Parvorder,"")] - | input == "superfamily" = [(Superfamily,"")] - | input == "family" = [(Family,"")] - | input == "subfamily" = [(Subfamily,"")] - | input == "supertribe" = [(Supertribe,"")] - | input == "tribe" = [(Tribe,"")] - | input == "subtribe" = [(Subtribe,"")] - | input == "genus" = [(Genus,"")] - | input == "subgenus" = [(Subgenus,"")] - | input == "section" = [(Section,"")] - | input == "series" = [(Series,"")] - | input == "superspecies" = [(Superspecies,"")] - | input == "species group" = [(Speciesgroup,"")] - | input == "species" = [(Species,"")] - | input == "species subgroup" = [(Speciessubgroup,"")] - | input == "subspecies" = [(Subspecies,"")] - | input == "infraspecies" = [(Infraspecies,"")] - | input == "varietas" = [(Variety,"")] - | input == "forma" = [(Form,"")] - | input == "no rank" = [(Norank,"")] - | otherwise = [(Norank,"")] - --- | Datastructure for entries of Taxonomy database dump nodes file -data TaxNode = TaxNode - { - -- node id in GenBank - taxId :: Int, - -- parent node id in GenBank taxonomy database - parentTaxId :: Int, - -- rank of this node (superkingdom, kingdom, ...) - rank :: Rank, - -- locus-name prefix; not unique - emblCode :: Maybe String, - -- see division.dmp file - nodeDivisionId :: String, - -- 1 if node inherits division from parent - inheritedDivFlag :: Bool, - -- see gencode.dmp file - nodeGeneticCodeId :: String, - -- 1 if node inherits genetic code from parent - inheritedGCFlag :: Bool, - -- see gencode.dmp file - mitochondrialGeneticCodeId :: String, - -- 1 if node inherits mitochondrial gencode from parent - inheritedMGCFlag :: Bool, - -- 1 if name is suppressed in GenBank entry lineage - genBankHiddenFlag :: Bool, - -- 1 if this subtree has no sequence data yet - hiddenSubtreeRootFlag :: Bool, - -- free-text comments and citations - nodeComments :: Maybe String - } - deriving (Show, Read, Eq) - --- | Simple Gene2Accession table -data SimpleGene2Accession = SimpleGene2Accession - { simpleTaxIdEntry :: Int, - simpleGenomicNucleotideAccessionVersion :: String - } deriving (Show, Eq, Read) - --- | Datastructure for Gene2Accession table -data Gene2Accession = Gene2Accession - { taxIdEntry :: Int, - geneID :: Int, - status :: String, - rnaNucleotideAccessionVersion :: String, - rnaNucleotideGi :: String, - proteinAccessionVersion :: String, - proteinGi :: String, - genomicNucleotideAccessionVersion :: String, - genomicNucleotideGi :: String, - startPositionOnTheGenomicAccession :: String, - endPositionOnTheGenomicAccession :: String, - orientation :: String, - assembly :: String, - maturePeptideAccessionVersion :: String, - maturePeptideGi :: String - } deriving (Show, Eq, Read) - -instance A.ToJSON (Gr SimpleTaxon Double) where - toJSON inputGraph = simpleTaxonJSONValue inputGraph 1 - -simpleTaxonJSONValue :: Gr SimpleTaxon Double -> Node -> A.Value -simpleTaxonJSONValue inputGraph node = jsonValue - where jsonValue = A.object [currentScientificName,T.pack "children" A..= children] - childNodes = suc inputGraph node - currentLabel = lab inputGraph node - currentScientificName = T.pack "name" A..= maybe (T.pack "notFound") (TL.toStrict . simpleScientificName) currentLabel - children = A.Array (V.fromList (map (simpleTaxonJSONValue inputGraph) childNodes)) - --jsonValue = A.object [currentScientificName,currentId,currentRank,(T.pack "children") A..= children] - --currentId = (T.pack "id") A..= (maybe (T.pack "notFound") (\a -> T.pack (show (simpleTaxId a))) currentLabel) - --currentRank = (T.pack "rank") A..= (maybe (T.pack "notFound") (\a -> T.pack (show (simpleRank a))) currentLabel)