liftovervcf without sorting the variants #1306

lindenb · 2019-04-10T10:49:06Z

Description

I'm trying to liftover gnomad/genome to hg38 but it takes to much time and memory. I'd rather have the resulting VCF file without sorting and then possibly sort later with bcftools.

I added an option DISABLE_SORT to liftover vcf:

now the sorter miight be null if DISABLE_SORT=true
I moved the variable out to a class member accept
I added a test that just write with/without DISABLE_SORT

now vcfliftover can write to stdout, variants are not sorted:

$ java -jar build/libs/picard.jar  LiftoverVcf INPUT=testdata/picard/vcf/LiftOver/testLiftoverBiallelicIndels.vcf OUTPUT=/dev/stdout CHAIN=testdata/picard/vcf/LiftOver/test.over.chain REJECT=/dev/null  REFERENCE_SEQUENCE=testdata/picard/vcf/LiftOver/dummy.reference.fasta DISABLE_SORT=true 2> /dev/null 
##fileformat=VCFv4.2
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
##INFO=<ID=ReverseComplementedAlleles,Number=0,Type=Flag,Description="The REF and the ALT alleles have been reverse complemented in liftover since the mapping from the previous reference to the current one was on the negative strand.">
##INFO=<ID=SwappedAlleles,Number=0,Type=Flag,Description="The REF and the ALT alleles have been swapped in liftover due to changes in the reference. It is possible that not all INFO annotations reflect this swap, and in the genotypes, only the GT, PL, and AD fields have been modified. You should check the TAGS_TO_REVERSE parameter that was used during the LiftOver to be sure.">
##contig=<ID=chr1,length=540>
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	Sample1
chr1	539	.	A	AAGGG	15676.17	PASS	END=539;ReverseComplementedAlleles	GT	0/0
chr1	421	.	CA	C	724.43	PASS	END=422;ReverseComplementedAlleles	GT	0/1
chr1	469	.	A	T	100	PASS	END=469;ReverseComplementedAlleles	GT	0/1
chr1	469	.	A	T	100	PASS	END=469;ReverseComplementedAlleles	GT	./.
chr1	421	.	CA	C	100	PASS	END=422;ReverseComplementedAlleles	GT	1/.

Checklist (never delete this)

Never delete this, it is our record that procedure was followed. If you find that for whatever reason one of the checklist points doesn't apply to your PR, you can leave it unchecked but please add an explanation below.

Content

Added or modified tests to cover changes and any new functionality
Edited the README / documentation (if applicable)
All tests passing on Travis

Review

Final thumbs-up from reviewer
Rebase, squash and reword as applicable

For more detailed guidelines, see https:/broadinstitute/picard/wiki/Guidelines-for-pull-requests

yfarjoun

Nice addition @lindenb! thanks!

Lifting gnomad is indeed a pain....it's interesting that you think that sorting it with a different tool will be better...

I've put in a few comments, but overall it looks great! thanks again.

src/main/java/picard/vcf/LiftoverVcf.java

yfarjoun · 2019-04-10T13:33:03Z

src/test/java/picard/util/LiftoverVcfTest.java

+
+ //try to open with / without index
+ try (final VCFFileReader liftReader = new VCFFileReader(liftOutputFile, !disableSort)) {
+ // nothing


at least count the variants and make sure that you got the right number....

yfarjoun · 2019-04-10T13:38:23Z

src/main/java/picard/vcf/LiftoverVcf.java

@@ -250,7 +253,10 @@
 );

 private VariantContextWriter rejects;
+ /** the output VariantContextWriter */
+ private VariantContextWriter accept;


For consistency, can you rename this accepts (with s) or change this and rejects to acceptedRecords and rejectedRecords ?

yfarjoun · 2019-04-10T13:39:04Z

src/main/java/picard/vcf/LiftoverVcf.java

+ .modifyOption(Options.ALLOW_MISSING_FIELDS_IN_HEADER, ALLOW_MISSING_FIELDS_IN_HEADER)
+ .setOutputFile(OUTPUT)
+ .setReferenceDictionary(walker.getSequenceDictionary())
+ ;


keep ; in same line.

yfarjoun · 2019-04-10T13:41:12Z

src/main/java/picard/vcf/LiftoverVcf.java

@@ -477,7 +508,11 @@ private void trackLiftedVariantContig(Map<String, Long> map, String contig) {
 private void addAndTrack(final VariantContext toAdd, final VariantContext source) {
 trackLiftedVariantContig(liftedBySourceContig, source.getContig());
 trackLiftedVariantContig(liftedByDestContig, toAdd.getContig());
- sorter.add(toAdd);
+ if (sorter != null) { //we're sorting the variants


Suggested change

if (sorter != null) { //we're sorting the variants

if (!DISABLE_SORT) { //we're sorting the variants

yfarjoun · 2019-04-10T13:43:23Z

src/test/java/picard/util/LiftoverVcfTest.java

+ // nothing
+ }
+ }
+


too many NLs

yfarjoun · 2019-04-10T13:47:07Z

src/main/java/picard/vcf/LiftoverVcf.java

- for (final VariantContext ctx : sorter) {
- out.add(ctx);
- progress.record(ctx.getContig(), ctx.getStart());
+


yfarjoun · 2019-04-10T13:49:30Z

src/main/java/picard/vcf/LiftoverVcf.java

+ log.info("Writing out sorted records to final VCF.");
+
+ for (final VariantContext ctx : sorter) {
+ this.accept.add(ctx);


indentation is off.

yfarjoun · 2019-04-10T13:50:02Z

src/main/java/picard/vcf/LiftoverVcf.java

@@ -477,7 +508,11 @@ private void trackLiftedVariantContig(Map<String, Long> map, String contig) {
 private void addAndTrack(final VariantContext toAdd, final VariantContext source) {
 trackLiftedVariantContig(liftedBySourceContig, source.getContig());
 trackLiftedVariantContig(liftedByDestContig, toAdd.getContig());
- sorter.add(toAdd);
+ if (sorter != null) { //we're sorting the variants
+ sorter.add(toAdd);


please use spaces (4) not tabs for indent.

src/test/java/picard/util/LiftoverVcfTest.java

lindenb · 2019-04-11T09:57:08Z

@yfarjoun back to you.

Thank you for the review , I hope I removed all the formatting problems

shouldn't the CREATE_INDEX be implicit?

well one can imagine that you want a sorted VCF without index. I added a test

        if (CREATE_INDEX && DISABLE_SORT) {
            log.error("CREATE_INDEX=true and DISABLE_SORT=true are mutually exclusive.");
            return 1;
         }

yfarjoun · 2019-04-11T13:13:43Z

what I meant is that your test only tests the explicit invocation, and not the implicit one...

yfarjoun

almost, a few extra newlines and one misformatted } else {

yfarjoun · 2019-04-11T13:12:52Z

src/main/java/picard/vcf/LiftoverVcf.java

@@ -287,6 +287,11 @@ protected int doWork() {
 IOUtil.assertFileIsWritable(OUTPUT);
 IOUtil.assertFileIsWritable(REJECT);

+
+ if (CREATE_INDEX && DISABLE_SORT) {


this type of input validation should be in customCommandLineValidation()

in fact I realized that CREATED_INDEX is only for BAM-related tools, for LiftoverVCF, the index is always created event if CREATED_INDEX=false

yfarjoun · 2019-04-11T13:14:08Z

src/main/java/picard/vcf/LiftoverVcf.java

@@ -181,6 +181,9 @@
 @Argument(doc = "INFO field annotations that should be deleted when swapping reference with variant alleles.", optional = true)
 public Collection<String> TAGS_TO_DROP = new ArrayList<>(LiftoverUtils.DEFAULT_TAGS_TO_DROP);

+ @Argument(doc = "Output VCF file will be written on the fly but it won't be sorted and indexed.", optional = true)


Suggested change

@Argument(doc = "Output VCF file will be written on the fly but it won't be sorted and indexed.", optional = true)

@Argument(doc = "Output VCF file will be written on the fly but it won't be sorted nor indexed.", optional = true)

“or” actually 😄

I'm not native english speaker, nor do I care that much :-D

yfarjoun · 2019-04-11T13:15:10Z

src/main/java/picard/vcf/LiftoverVcf.java

- out.writeHeader(outHeader);
+ this.acceptedRecords = new VariantContextWriterBuilder()
+ .modifyOption(Options.ALLOW_MISSING_FIELDS_IN_HEADER, ALLOW_MISSING_FIELDS_IN_HEADER)
+ .modifyOption(Options.INDEX_ON_THE_FLY,!DISABLE_SORT)


Suggested change

.modifyOption(Options.INDEX_ON_THE_FLY,!DISABLE_SORT)

.modifyOption(Options.INDEX_ON_THE_FLY, !DISABLE_SORT)

lindenb · 2019-04-11T13:57:46Z

what I meant is that your test only tests the explicit invocation, and not the implicit one...

Sorry, I'm not sure I understand 🤔 . for example in a test like LiftoverVcfTest/testWriteOriginalPosition how is the condition implicit/explicit tested ?

yfarjoun · 2019-04-12T16:17:59Z

I guess I mean that you explicitly told the program to CREATE_INDEX rather than it doing whatever it thinks

liftovervcf without sorting

2559a90

yfarjoun suggested changes Apr 10, 2019

View reviewed changes

answer to broadinstitute#1306

ca771b3

yfarjoun approved these changes Apr 11, 2019

View reviewed changes

yfarjoun merged commit a6c6d88 into broadinstitute:master Apr 22, 2019

lindenb deleted the pl_vcfliftover_nosort branch April 23, 2019 09:14

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

liftovervcf without sorting the variants #1306

liftovervcf without sorting the variants #1306

lindenb commented Apr 10, 2019

yfarjoun left a comment

yfarjoun Apr 10, 2019

yfarjoun Apr 10, 2019

yfarjoun Apr 10, 2019

yfarjoun Apr 10, 2019

yfarjoun Apr 10, 2019

yfarjoun Apr 10, 2019

yfarjoun Apr 10, 2019

yfarjoun Apr 10, 2019

lindenb commented Apr 11, 2019

yfarjoun commented Apr 11, 2019

yfarjoun left a comment

yfarjoun Apr 11, 2019

lindenb Apr 11, 2019

yfarjoun Apr 11, 2019

jmarshall Apr 11, 2019

yfarjoun Apr 11, 2019

yfarjoun Apr 11, 2019

lindenb commented Apr 11, 2019

yfarjoun commented Apr 12, 2019

	if (sorter != null) { //we're sorting the variants
	if (!DISABLE_SORT) { //we're sorting the variants

	@Argument(doc = "Output VCF file will be written on the fly but it won't be sorted and indexed.", optional = true)
	@Argument(doc = "Output VCF file will be written on the fly but it won't be sorted nor indexed.", optional = true)

	.modifyOption(Options.INDEX_ON_THE_FLY,!DISABLE_SORT)
	.modifyOption(Options.INDEX_ON_THE_FLY, !DISABLE_SORT)

liftovervcf without sorting the variants #1306

liftovervcf without sorting the variants #1306

Conversation

lindenb commented Apr 10, 2019

Description

Checklist (never delete this)

Content

Review

yfarjoun 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

lindenb commented Apr 11, 2019

yfarjoun commented Apr 11, 2019

yfarjoun 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

lindenb commented Apr 11, 2019

yfarjoun commented Apr 12, 2019