lightning-auriga / interpolate-genetic-position Goto Github PK

This program now needs to correctly emit output for input query bedfiles with non-contiguous regions, e.g.

chr22	17083846	17084000	thing1
chr22	17084050	17084145	thing2

Assuming an input genetic map (bedgraph format, base 0) like this:

chr22	17076254	17081023	0	0
chr22	17081023	17083846	0.137252	0
chr22	17083846	17084017	4.00811	0.000387462396
chr22	17084017	17084145	1.91398e-06	0.001072849206
chr22	17084145	17086809	0.925629	0.00107284945098944
chr22	17086809	17087057	7.31739e-12	0.00353872510698944

The bolt format genetic map would currently look something like (bolt format, base 1):

chr	position	COMBINED_rate(cM/Mb)	Genetic_Map(cM)
chr22	17083847	4.00811	0.000387462
chr22	17084051	1.91398e-06	0.00107285

This is incorrect. There needs to be an additional breakpoint in the output, corresponding to base-0 position 17084017, with rate 1.91398e-06. Otherwise, the first block incorrectly implies linearly interpolation with the first rate, 4.00811.

input file operations on variant query files have already been refactored into a derived class, whose base class can be mocked. this same modification needs to be performed for:

• input file operations on variant query files (interpolate-genetic-position/input_variant_file)
• input file operations on genetic maps (interpolate-genetic-position/input_genetic_map_file)
• output file operations on (modified) variant query files

it runs, but by the time you get to the mitochondrial variants of the test file, nonsense gpos numbers are emitted. by design, values on a chromosome absent from the genetic map should be uniformly 0.

while the gtest infrastructure is now present (in branch), there need to be unit tests for primary functionality.

I think it would be cool for the program to accept vcf input and emit vcf output. That way, you could potentially add it to a pipe chain with bcftools.

vcf operations would happen via htslib, and the genetic position would presumably be annotated as a new INFO field.

Part of the point of --region-step-interval was not to necessarily increment every query boundary, but only ones where the bed+4 column 4 fit some criterion. It seems like the criterion that's needed now is that the fourth column of two successive query regions are not the same, after trimming the suffix -.*$.

need to modify the .gitattributes file to suppress m4 files, in particular. this should also get propagated upstream to the cookiecutter.