比较Clojure中的两个大文件(即;在顶帽对齐中查找未映射的读取)

Question

问题:找到一个文件中但不存在于另一个文件中的ID.每个文件大约6.5 GB.具体而言(对于生物信息学领域中的那些),一个文件是测序读取的fastq文件,另一个是来自tophat运行的sam对齐文件.我想确定fastq文件中的哪些读取不在sam对齐文件中.

我收到了java.lang.OutOfMemory: Java heap space错误.正如所建议的那样(ref1,ref2)我正在使用惰性序列.但是,我的内存仍然不足.我看过这个教程,但我还不太了解它.所以我发布了一个不太复杂的尝试解决方案,希望我只是犯了一个小错误.

我的尝试:

由于这两个文件都不适合内存,因此sam文件中的行一次读取一个块,并将块中每行的id放入一个集合中.然后使用集合中的sam id过滤一个懒惰的fastq id列表,只保留那些不在集合中的id.使用下一个sam行和剩余的fastq id重复此操作.

(defn ids-not-in-sam 
  [ids samlines chunk-size]
  (lazy-seq
    (if (seq samlines)
      (ids-not-in-sam (not-in (into #{} (qnames (take chunk-size samlines))) ids)
                      (drop chunk-size samlines) chunk-size)
      ids)))

not-in 确定哪些ID不在集合中.

(defn not-in 
  ; Return the elements x of xs which are not in the set s
  [s xs]
  (filter (complement s) xs))

qnames 从sam文件中的一行获取id字段.

(defn qnames [samlines]
  (map #(first (.split #"\t" %)) samlines))

最后,它与io放在一起(使用read-lines和write-lines来自clojure.contrib.io.

(defn write-fq-not-in-sam [fqfile samfile fout chunk-size] 
    (io/write-lines fout (ids-not-in-sam (map fq-id (read-fastq fqfile))
                                         (read-sam samfile) chunk-size))) 

我很确定我是以懒惰的方式做所有事情.但我可能会抓住某个我不注意的序列的头部.

上面的代码中是否有错误导致堆填满？更重要的是,我对这个问题的处理方法都错了吗？这是否适用于懒惰序列,我是否期望过多？

(错误可能在read-sam和read-fastq函数中,但我的帖子已经有点长了.如果需要,我可以稍后展示它们).

Answer 1

1. 对数据集进行排序(将块读入内存,排序,写入临时文件,合并已排序的文件)
2. 遍历两个集合以找到相交/缺失元素(我希望算法很清楚)