Rewrite pileup module #140

alumi · 2018-08-10T10:03:21Z

Summary

BREAKING CHANGES
Completely rewrote pileup module in order to improve its performance and usability.

Changes

Don't pile up uncovered regions
- Sparsely piled-up alignments
Extract I/O from pileup algorithm
- Reader/Writer of pileup format
- Improved serialization with start/end (^/$) and mapq support
Move everything in cljam.pileup.mpileup to cljam.pileup
- Remove cljam.pileup.common which is no longer referred
Add a benchmarking script

Performance comparison

Piling up chr1:1-247249719 of large.bam with a single thread.

master

time: 2301.451698 sec, sd: 40390.179732 sec

feature/sparse-mpileup

time: 48.248668 sec, sd: 648.517263 µs

Affects

cljam.algo.pileup
cljam.io.sam.util.cigar
cljam.io.sam.util.quality

Tests

lein check 🆗
lein test :all 🆗

Notes

It could be even faster if we implement the eager version of pileup like cljam.algo.depth/depth
Upgraded cloverage to 1.0.11

codecov · 2018-08-15T15:06:40Z

Codecov Report

Merging #140 into master will decrease coverage by 0.08%.
The diff coverage is 75.16%.

@@            Coverage Diff            @@
##           master    #140      +/-   ##
=========================================
- Coverage   85.58%   85.5%   -0.09%     
=========================================
  Files          69      68       -1     
  Lines        4474    4580     +106     
  Branches      419     444      +25     
=========================================
+ Hits         3829    3916      +87     
+ Misses        226     220       -6     
- Partials      419     444      +25

Impacted Files	Coverage Δ
src/cljam/io/sam/util/quality.clj	`72.41% <0%> (-10.35%)`	⬇️
src/cljam/io/sam/util/cigar.clj	`74.02% <100%> (ø)`	⬆️
src/cljam/io/pileup.clj	`58.06% <58.06%> (ø)`
src/cljam/tools/cli.clj	`59.03% <58.06%> (+5.74%)`	⬆️
src/cljam/algo/pileup.clj	`95.74% <95.65%> (-4.26%)`	⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 125ee30...bcbc66f. Read the comment docs.

totakke

Thank you for amazing performance improvement.

However, strange Ġ seems to be written by pileup (~ is correct).

$ bash -c 'diff <(lein run pileup test-resources/bam/small.bam) <(samtools mpileup test-resources/bam/small.bam) | head'
1c1
< chr1  23000088        N       1       ^ĠT     H
---
> chr1  23000088        N       1       ^~T     H
15c15
< chr1  23000102        N       2       G^ĠG    HH
---
> chr1  23000102        N       2       G^~G    HH
30c30
< chr1  23000117        N       3       GG^Ġt   EEE

I added some other comments. Please check them.

totakke · 2018-08-16T08:33:45Z

test/cljam/algo/pileup_test.clj

+(deftest pileup-region
+  (with-open [br (sam/bam-reader test-sorted-bam-file)]
+    (let [plp-ref1 (doall (plp/pileup br {:chr "ref" :start 1 :end 40}))
+          plp-ref2 (doall (plp/pileup br {:chr "ref2" :start 1 :end 40}))]


plp-ref2 is never used. I'm not sure whether plp-ref2 is unnecessary or tests are insufficient.

totakke · 2018-08-16T08:37:15Z

src/cljam/tools/cli.clj


 (defn pileup [args]
-  (let [{:keys [options arguments errors summary]} (parse-opts args pileup-cli-options)]
+  (let [{:keys [options arguments errors summary]


options in this line can be removed.

totakke · 2018-08-16T08:41:05Z

src/cljam/tools/cli.clj

-            [cljam.util.region :as region])
+            [cljam.algo.pileup :as plp]
+            [cljam.util.region :as region]
+            [clojure.java.io :as cio])
  (:import [java.io Closeable BufferedWriter OutputStreamWriter]))


[cljam.io.sequence :as cseq] is no longer used.

totakke · 2018-08-16T08:59:25Z

test/cljam/algo/pileup_test.clj

+            [5 [{:pos 5, :pile [{:pos 3, :end 5}]} {:pos 5, :pile [{:pos 4, :end 6}]}]]
+            [6 [nil {:pos 6, :pile [{:pos 4, :end 6}]}]]]
+           (plp/align-pileup-seqs [{:pos 3 :pile [{:pos 3 :end 5}]}
+                                    {:pos 4 :pile [{:pos 3 :end 5}]}


Wrong indent

totakke · 2018-08-16T08:59:58Z

test/cljam/algo/pileup_test.clj

+            [8 [nil {:pos 8, :pile [{:pos 8, :end 9}]}]]
+            [9 [nil {:pos 9, :pile [{:pos 8, :end 9}]}]]]
+           (plp/align-pileup-seqs [{:pos 3 :pile [{:pos 3 :end 5}]}
+                                    {:pos 4 :pile [{:pos 3 :end 5}]}


Wrong indent

totakke · 2018-08-16T09:00:07Z

test/cljam/algo/pileup_test.clj

+    (is (= [[3 [{:pos 3, :pile [{:pos 3, :end 4}]} nil {:pos 3, :pile [{:pos 3, :end 3}]} nil]]
+            [4 [{:pos 4, :pile [{:pos 3, :end 4}]} {:pos 4, :pile [{:pos 4, :end 4}]} nil nil]]]
+           (plp/align-pileup-seqs [{:pos 3 :pile [{:pos 3 :end 4}]}
+                                    {:pos 4 :pile [{:pos 3 :end 4}]}]


Wrong indent

totakke · 2018-08-16T09:33:42Z

src/cljam/tools/cli.clj

+(defn- depth
+  [f region n-threads]
+  (with-open [r (sam/reader f)
+              w (cio/writer *out*)]


What is this w used for ...?

This line causes IOException:

$ lein run pileup -s test-resources/bam/test.sorted.bam ... Exception in thread "main" java.io.IOException: Stream closed

totakke · 2018-08-16T09:43:36Z

src/cljam/tools/cli.clj

+      (if simple
+        (depth f region thread)
+        (with-open [w (cio/writer (cio/output-stream System/out))]
+          (plp/create-mpileup f ref w (some-> region parse-region)))))))


parse-region returns nil if arg is nil, so some-> is not necessary.

totakke · 2018-08-16T10:19:49Z

src/cljam/algo/pileup.clj

+            [cljam.io.sam.util.quality :as qual]
+            [cljam.io.sam.util.refs :as refs]
+            [cljam.io.pileup :as plpio])
+  (:import [java.io Closeable]


No Closeable exists in this ns.

totakke · 2018-08-16T11:34:23Z

test/cljam/io/pileup_test.clj

+                   (cio/delete-file (cio/file tmp)))))))
+      "chr1\t10\tA\t1\t.\tI\n"
+      "chr1\t10\tA\t4\t.,Tt\tIABC\n"
+      "chr1\t10\tA\t4\t^].+3TTT,-2tg$Tt\tIABC\n")))


I'd like you to add source-type-test of reader/writer in the same way as other io tests.

alumi · 2018-08-17T03:38:52Z

@totakke Thanks for your helpful review! 🙇

I added some commits to fix problems you pointed. The Ġ was caused by encoding MAPQ 255 into phred33 char without range checks, now fixed with c72c48d.

bash -c 'diff <(lein run pileup test-resources/bam/small.bam -r chr1:23000088-23000117) <(samtools mpileup test-resources/bam/small.bam -r chr1:23000088-23000117) | head' 2>/dev/null
# no diffs in the region

athos

Thank you for your great effort to improve the pileup performance. It would be really helpful to analyze a large dataset.

One thing I would like to confirm is if PileupBase should have the alignment field. A .pileup file doesn't necessarily contain enough information to restore a complete alignment back from it, so it feels a little weird that PileupBase has that field.
While I'm aware that the qname field of the alignment is used to correct overlapped reads, it seems to me more intuitive that PileupBase has the qname field directly, instead of the alignment field.

Perhaps the field can be used to generate VCF with CIGAR INFO field (or any other information that isn't stored explicitly in .pileup files), but even in that case, I think the alignment field could be an optional field that isn't declared in the PileupBase record definition.

athos · 2018-08-17T04:33:16Z

src/cljam/algo/pileup.clj

+                   (keep (partial ->locus-pile chr)))))))))
+
+(defn align-pileup-seqs
+  "Align multiple pileed-up seqs."


s/pileed-up/piled-up/

athos

Also added some trivial comments.

athos · 2018-08-17T04:36:13Z

src/cljam/algo/pileup.clj

+           (<= min-mapq (.mapq aln))))))
+
+(defn resolve-base
+  "Find a piled-up base and an indel from a alignment."


"an alignment" or "the alignment"

- Make :alignment optional - Add :qname to store a name of a query sequence

athos · 2018-08-17T05:23:54Z

src/cljam/io/pileup.clj

+;; Writer
+;; ------
+
+(defn- ^String stringify-mpileup-alignment


I think the pileup writer could be more efficient if it writes out the pileup contents directly to the BufferedWriter with java.io.Writer#append or something, instead of once stringifying the contents and writing them out to the writer.

.pileup files could be GB size in some cases, so saving the constructions of StringBuilders and a plenty of stringified contents could drastically reduce the time taken to write out a whole pileup, although I think it's also OK to leave this as a further TODO improvement.

alumi · 2018-08-17T06:35:12Z

@athos Thanks for your review!! 🙇
I added alignment field of PileupBase just for a use of the direct accessor .alignment.
So I change it to an extra entry and add qname field instead.

totakke · 2018-08-17T09:38:56Z

Thank you!

athos · 2018-08-17T09:40:56Z

Thank you for tackling a tough work, @alumi!

alumi added 7 commits July 17, 2018 15:16

Upgrade clojure 1.10 to 1.10.0-alpha6

5a148bf

Upgrade log4j to 2.11.1

589db7a

Upgrade lein-cloverage to 1.0.11

87a26f5

Add a benchmarking script for pileup

a3ff820

Rewrite mpileup module to improve the performance.

8182f43

Extract I/O from pileup

0682795

Remove namespace cljam.algo.pileup.mpileup

3bc32e7

alumi added the enhancement label Aug 10, 2018

alumi assigned athos and totakke Aug 10, 2018

alumi requested review from athos and totakke August 10, 2018 10:03

totakke requested changes Aug 16, 2018

View reviewed changes

alumi added 9 commits August 17, 2018 10:36

Add missing tests for piling up alignments on ref2

ed7c0ef

Remove an unnecessary requirement

fad8380

Remove an unnecessary binding

c22c893

Fix wrong indentations in test

9a79ef7

Remove an unnecessary import

908d055

Remove an unnecessary nil check

eec0106

Add source-type-test for pileup file format

9e9787a

Removed an unnecessary writer

5ebbc3e

Clamp quality values when encoding them into chars

c72c48d

athos reviewed Aug 17, 2018

View reviewed changes

athos requested changes Aug 17, 2018

View reviewed changes

alumi added 2 commits August 17, 2018 14:06

Fix typos

a57abc7

Change instance variables of cljam.io.pileup.PileupBase

b7a95d2

- Make :alignment optional - Add :qname to store a name of a query sequence

athos reviewed Aug 17, 2018

View reviewed changes

Directly write piles to a writer instead of building strings

bcbc66f

alumi force-pushed the feature/sparse-mpileup branch from c6ad9b3 to bcbc66f Compare August 17, 2018 06:26

athos approved these changes Aug 17, 2018

View reviewed changes

totakke approved these changes Aug 17, 2018

View reviewed changes

totakke merged commit 0a5b86c into master Aug 17, 2018

totakke deleted the feature/sparse-mpileup branch August 17, 2018 09:38

athos mentioned this pull request May 28, 2019

Improve pileup more #165

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Rewrite pileup module #140

Rewrite pileup module #140

alumi commented Aug 10, 2018 •

edited

codecov bot commented Aug 15, 2018 •

edited

totakke left a comment

totakke Aug 16, 2018

totakke Aug 16, 2018

totakke Aug 16, 2018

totakke Aug 16, 2018

totakke Aug 16, 2018

totakke Aug 16, 2018

totakke Aug 16, 2018

totakke Aug 16, 2018

totakke Aug 16, 2018

totakke Aug 16, 2018

alumi commented Aug 17, 2018

athos left a comment

athos Aug 17, 2018

athos left a comment

athos Aug 17, 2018

athos Aug 17, 2018

alumi commented Aug 17, 2018

totakke commented Aug 17, 2018

athos commented Aug 17, 2018

Rewrite pileup module #140

Rewrite pileup module #140

Conversation

alumi commented Aug 10, 2018 • edited

Summary

Changes

Performance comparison

master

feature/sparse-mpileup

Affects

Tests

Notes

codecov bot commented Aug 15, 2018 • edited

Codecov Report

totakke left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alumi commented Aug 17, 2018

athos left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

athos left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alumi commented Aug 17, 2018

totakke commented Aug 17, 2018

athos commented Aug 17, 2018

alumi commented Aug 10, 2018 •

edited

codecov bot commented Aug 15, 2018 •

edited