2

我有一个由我创建的 ruby​​ 脚本创建的 fasta 文件；但是，它有许多重复的序列，它们在不同的序列 ID 下，我想展平文件，以便我可以从这个...

输入 fasta 文件 (example.fa) 中的数据

>isotig00009_f3.4.1 
ITLKPCGVPFSCCIPDQASGVANTQCGYGVRSPEQQNTFHTKIYTTGCADMFTMWINRYLYYIAGIAGVIVLVELFGFCFAHSLINDIKRQKARWAHR
>isotig00037_f3.1.1 
KLSLIVVNHSMVASKFERVILAYTSIIIEVSPAKRRRNGKILRKNTIRFRWQTFRILSAFTVVTFSKMMTQKI
>isotig00045_f1.15.3
YKINKRP
>isotig00046_f3.15.3
YKINKRP
>isotig00047_f3.15.3
YKINKRP
>isotig00048_f1.15.3
YKINKRP
>isotig00049_f1.15.3
YKINKRP
>isotig00050_f2.15.3
YKINKRP
>isotig00051_f1.15.3
YKINKRP

对此（在输出文件 - output.fa 中）...

>isotig00009_f3.4.1 
ITLKPCGVPFSCCIPDQASGVANTQCGYGVRSPEQQNTFHTKIYTTGCADMFTMWINRYLYYIAGIAGVIVLVELFGFCFAHSLINDIKRQKARWAHR
>isotig00037_f3.1.1 
KLSLIVVNHSMVASKFERVILAYTSIIIEVSPAKRRRNGKILRKNTIRFRWQTFRILSAFTVVTFSKMMTQKI
>isotig00045_f1.15.3 : isotig00046_f3.15.3 : isotig00047_f3.15.3 : isotig00048_f1.15.3 : isotig00049_f1.15.3 : isotig00050_f2.15.3 : isotig00051_f1.15.3
YKINKRP

我制作了一个小脚本，取出序列的一个副本，但尽管尝试了很长时间，但我似乎无法为每个序列添加任何序列 ID。我尝试过但不起作用的一件事（见下面的评论部分），是尝试提取唯一序列，然后提取 sequence_IDs ......

这是我正在使用的脚本：

#!/usr/bin/env ruby

filename = "./example.fa"
text = File.read(filename)

def seq_uniq(input, output)
    parser = /^>.*\n(.*)/i
    seq_id_parser = /^(>.*)\n(.*)/i
    file = File.new("#{output}", "w")
    input.scan(parser).uniq.each do |seq|
        file.puts seq
#        input.scan(seq_id_parser) do |seq_id, seq_actual|
#           if seq_actual == seq
#               file.puts seq_id
#           end
#        end
    end
    file.close
end


seq_uniq(text, "./output.fa")

如果有人能指出我正确的方向，我将不胜感激。由于这需要嵌入到生成 html 网页的大型 ruby​​ 脚本中，如果您只使用 ruby​​（或

编辑：

为了澄清，我想展平 fasta 文件以放置所有 Seq. 在单行上标识相同序列的 ID（由 : 分隔），然后在下一行具有序列。

非常感谢

4

1 回答 1

1

您可以读取 fasta 文件并创建一个“sequence => list of ids with same sequence”的散列。请注意，此解决方案假定输入文件每行只有一个 id，如果有多个，则需要改进解析器：

def uniq_sequences(input)
  ids_by_sequence = {}
  input.split(/\n/).each_slice(2) do |id, sequence|
    id = id.gsub(/^>\s*/,'') # Remove the leading ">".
    ids_by_sequence[sequence] = [] unless ids_by_sequence[sequence]
    ids_by_sequence[sequence] << id # Store all ids with the same sequence.
  end
  # Return a string of aggregated ids per same sequence.
  ids_by_sequence.map do |sequence, ids|
    ">#{ids.join(' : ')}\n#{sequence}"
  end.join("\n")
end

puts uniq_sequences(File.read('./example.fa'))

在您的示例输入中，您应该得到以下输出：

>isotig00009_f3.4.1 
ITLKPCGVPFSCCIPDQASGVANTQCGYGVRSPEQQNTFHTKIYTTGCADMFTMWINRYLYYIAGIAGVIVLVELFGFCFAHSLINDIKRQKARWAHR
>isotig00037_f3.1.1 
KLSLIVVNHSMVASKFERVILAYTSIIIEVSPAKRRRNGKILRKNTIRFRWQTFRILSAFTVVTFSKMMTQKI
>isotig00045_f1.15.3 : isotig00046_f3.15.3 : isotig00047_f3.15.3 : isotig00048_f1.15.3 : isotig00049_f1.15.3 : isotig00050_f2.15.3 : isotig00051_f1.15.3
YKINKRP
于 2013-08-26T23:30:41.960 回答