haskell - BioHaskell:读取 FASTA 文件

Question

使用 BioHaskell , 如何读取包含氨基酸序列的 FASTA 文件？

我希望能够:

• 获取String序列的列表
• 从 FASTA 注释(假设是唯一的)到序列字符串中获取 Map String String(来自 Data.Map)
• 在 BioHaskell 中实现的算法中使用序列。

注意:这个问题有意不显示研究成果，因为它会立即以问答式的方式回答。

Answer 1

提取原始序列字符串

从现在开始我们假设文件 aa.fa 包含一些氨基酸 FASTA 序列。让我们从一个提取序列列表的简单示例开始。

import Bio.Sequence.Fasta (readFasta)
import Bio.Sequence.SeqData (seqdata)
import qualified Data.ByteString.Lazy.Char8 as LB

main = do
    sequences <- readFasta "aa.fa"
    let listOfSequences = map (LB.unpack . seqdata) sequences :: [String]
    -- Just for show, we will print one sequence per line here
    -- This will basically execute putStrLn for each sequence
    mapM_ putStrLn listOfSequences

readFasta返回 IO [Sequence Unknown]。基本上这意味着没有关于序列是否包含氨基酸或核苷酸的信息。

请注意，我们使用 LB.unpack而不是此处的 show，因为 show 在生成的 String 的开头和结尾添加了双引号 (") >。使用 LB.unpack 是可行的，因为在当前的 BioHaskell 0.5.3 版本中，SeqData 只是定义为惰性 ByteString。

我们可以使用 castToAmino 来解决这个问题或 castToNuc :

转换为 AA/核苷酸序列

let aaSequences = map castToAmino sequences :: [Sequence Amino]

请注意，这些函数目前(BioHaskell 版本 0.5.3)不执行任何有效性检查。您可以在 BioHaskell 算法中使用 [Sequence Amino] 或 [Sequence Nuc]。

通过FASTA header 查找序列

我们现在假设我们的 aa.fa 包含一个序列

>abc123
MGLIFARATNA...

现在，我们将从 FASTA 文件构建一个 Map String String(我们将在本例中使用 Data.Map.Strict)。我们可以使用这个映射来查找序列。

查找将产生一个Maybe String。此示例中的预期行为是在找到序列时打印序列，或者如果在 Map 中未找到任何内容则不打印任何内容。

由于 Data.Maybe 是一个 Monad，我们可以使用 Data.Foldable.mapM_为了这个任务。

import Bio.Sequence.Fasta (readFasta)
import Bio.Sequence.SeqData (Sequence, seqdata, seqheader)
import qualified Data.ByteString.Lazy.Char8 as LB
import Data.Foldable (mapM_)
import qualified Data.Map.Strict as Map

-- | Convert a Sequence to a String tuple (sequence label, sequence)
sequenceToMapTuple :: Sequence a -> (String, String)
sequenceToMapTuple s = (LB.unpack $ seqheader s, LB.unpack $ seqdata s)

main = do
    sequences <- readFasta "aa.fa"
    -- Build the sequence map (by header)
    let sequenceMap = Map.fromList $ map sequenceToMapTuple sequences
    -- Lookup the sequence for the "abc123" header
    mapM_ print $ Map.lookup "abc123" sequenceMap

编辑:感谢@GabrielGonzalez 的建议，最后一个示例现在使用 Data.Foldable.mapM_ 而不是 Data.Maybe.fromJust