bedcov header to have optional headers as well#2140
Conversation
jkbonfield
left a comment
There was a problem hiding this comment.
I tested a bunch of things and it looks to work with valid input.
Invalid data with differing numbers of columns in the records and a pre-existing header do not correct things so the output is just as broken as the input. I think that's fine and attempting to do anything more is a total can of worms.
| .R Print a comment/header describing columns. When a header is available in | ||
| bed file, it is used in output, when it is not available, mandatory field's | ||
| header and empty tab separated header for optional fields are used in output. |
There was a problem hiding this comment.
This is a bit of a complex sentence.
I'd suggest breaking it up.
When a header starting in "#chrom" is available in the input bed file, it is copied to the output. When it is not available a header is created with field names matching the fields listed in the GA4GH BED specification.
The \fB-c\fR and \fB-j\fR options can add further per-file columns named
.RI in1.sam "_cov, " in1.sam "_depth and " in1.sam "_count"
| cmd("echo \"#chrom\tchromStart\tchromEnd\t$$opts{path}/bedcov/bedcov.bam_cov\" > $$opts{tmp}/bedcovH1.expected"); | ||
| cmd ("cat $$opts{path}/bedcov/bedcov.expected >> $$opts{tmp}/bedcovH1.expected"); | ||
| cmd("$$opts{bin}/samtools bedcov -H $$opts{path}/bedcov/bedcov.bed $$opts{path}/bedcov/bedcov.bam > $$opts{tmp}/out.H1"); | ||
| $out = cmd("diff $$opts{tmp}/bedcovH1.expected $$opts{tmp}/out.H1"); | ||
| if ( $out ne "" ) { | ||
| failed($opts,msg=>"coverage with header",reason=>"output does $out not match to expected\n"); | ||
| } else { | ||
| passed($opts,msg=>"coverage with header success\n"); | ||
| } | ||
| #with custom header | ||
| cmd("echo \"#chrom\tchromStart\tchromEnd\tT1\nchr1\t12209228\t12209246\t10\" > $$opts{tmp}/bedcovH2.bed"); | ||
| cmd("echo \"#chrom\tchromStart\tchromEnd\tT1\t$$opts{path}/bedcov/bedcov.bam_cov\nchr1\t12209228\t12209246\t10\t24\" > $$opts{tmp}/bedcovH2.expected"); | ||
| cmd("$$opts{bin}/samtools bedcov -H $$opts{tmp}/bedcovH2.bed $$opts{path}/bedcov/bedcov.bam > $$opts{tmp}/out.H2"); | ||
| $out = cmd("diff $$opts{tmp}/bedcovH2.expected $$opts{tmp}/out.H2"); | ||
| if ( $out ne "" ) { | ||
| failed($opts,msg=>"coverage with custom header",reason=>"output does not match to expected\n"); | ||
| } else { | ||
| passed($opts,msg=>"coverage with custom header success\n"); | ||
| } | ||
| #with empty source header | ||
| cmd("echo \"#chrom\tchromStart\tchromEnd\t\nchr1\t12209228\t12209246\t10\" > $$opts{tmp}/bedcovH3.bed"); | ||
| cmd("echo \"#chrom\tchromStart\tchromEnd\t\t$$opts{path}/bedcov/bedcov.bam_cov\nchr1\t12209228\t12209246\t10\t24\" > $$opts{tmp}/bedcovH3.expected"); | ||
| cmd("$$opts{bin}/samtools bedcov -H $$opts{tmp}/bedcovH3.bed $$opts{path}/bedcov/bedcov.bam > $$opts{tmp}/out.H3"); | ||
| $out = cmd("diff $$opts{tmp}/bedcovH3.expected $$opts{tmp}/out.H3"); | ||
| if ( $out ne "" ) { | ||
| failed($opts,msg=>"coverage with src empty header",reason=>"output does not match to expected\n"); | ||
| } else { | ||
| passed($opts,msg=>"coverage with src empty header success\n"); | ||
| } | ||
| #empty added header | ||
| #chrom\tchromStart\tchromEnd\tname\tscore\tstrand\tthickStart\t\ | ||
| # thickEnd\titemRgb\tblockCount\tblockSizes\tblockStarts | ||
| cmd("echo \"chr1\t12209228\t12209246\t4\t5\t6\t7\t8\t9\t10\t11\t12\t13\t14\" > $$opts{tmp}/bedcovH4.bed"); | ||
| cmd("echo \"#chrom\tchromStart\tchromEnd\tname\tscore\tstrand\tthickStart\tthickEnd\titemRgb\tblockCount\tblockSizes\tblockStarts\t\t\t$$opts{path}/bedcov/bedcov.bam_cov\" > $$opts{tmp}/bedcovH4.expected"); | ||
| cmd("echo \"chr1\t12209228\t12209246\t4\t5\t6\t7\t8\t9\t10\t11\t12\t13\t14\t24\" >> $$opts{tmp}/bedcovH4.expected"); | ||
| cmd("$$opts{bin}/samtools bedcov -H $$opts{tmp}/bedcovH4.bed $$opts{path}/bedcov/bedcov.bam > $$opts{tmp}/out.H4"); | ||
| $out = cmd("diff $$opts{tmp}/bedcovH4.expected $$opts{tmp}/out.H4"); | ||
| if ( $out ne "" ) { | ||
| failed($opts,msg=>"coverage with empty header",reason=>"output does not match to expected\n"); | ||
| } else { | ||
| passed($opts,msg=>"coverage with empty header success\n"); | ||
| } |
There was a problem hiding this comment.
Why are all these expected files created with echo? The other existing test data just has these files checked in. It's much easier to read than a script, and permits us to experiment on the command line outside of running the test harness.
There was a problem hiding this comment.
The header for result fields are as "< file path >_cov" / "< file path >_count"...
This file path is random path and hence dynamic creation of expected files.
There was a problem hiding this comment.
Fair enough. Either we have custom code for filtering out differences when checking or we have custom code for generating the validation data. Happy with your solution as the header is specifically something we want to be checking. Thanks for explaining it.
| "blockSizes", "blockStarts"}; | ||
| for (i = 0; i < fields; ++i) { | ||
| fprintf(fp, "%s%s", (i ? "\t" : "#"), | ||
| (i < sizeof(bedcols)/sizeof(bedcols[i]) ? bedcols[i] : "")); |
There was a problem hiding this comment.
It's good that you check for overflow here, but maybe instead of "" we should add . or ? so the separate words are still visible on a terminal and automatic parsing with things like sscanf or awk print $5 etc counts the header fields correctly.
It's a very minor point though as such data doesn't match the specificaton anyway.
There was a problem hiding this comment.
Using '.' instead of empty headers.
66b635c to
e926f63
Compare
Adds header for optional columns from bed file along with mandatory column headers.
If a header is present in bed file, it is used as such.
(A comment line starting with #chrom is considered as header in bed file.)
If header is not present in bed file, a header is made with text as mentioned in format (https://genome.ucsc.edu/FAQ/FAQformat#format1). There is difference from earlier version where "chrom tab start tab end" was used and now they are "chrom tab chromStart tab chromEnd".
If there are more fields than the format specified, empty header with tab separation is used.
Result column headers appended after this header.
fixes #2126