Msprime: RecombinationMap class

[andrewkern]

Hey all- so I ran in to an edge case today as part of the stdpopsim effort where the RecombinationMap class was choking on a map which had nonzero recombination rate in the last interval. Is there a reason this requirement is being kept? There are many systems where ends of chromosomes will be crossing over, so it would be great if we could drop this.

Andy

[jeromekelleher]

I've had a quick look into the details here, and I think the zero-end requirement is a mistake. The docs say this:

A RecombinationMap represents the changing rates of recombination along a chromosome. This is defined via two lists of numbers: positions and rates, which must be of the same length. Given an index j in these lists, the rate of recombination per base per generation is rates[j] over the interval positions[j] to positions[j + 1]. Consequently, the first position must be zero, and by convention the last rate value is also required to be zero (although it is not used).

I think what this should say is that the rate is that rate[j] covers the interval positions[j - 1] to positions[j]. So, the initial value should be zero by convention (but is ignored) and there is no requirement on the last position. Do you agree with this @petrelharp?

The recombination map code has some quirks, I was never particularly happy with it.

While we're sorting this out, we should also fix #39 and modernise the recombination map class to use numpy.

[jeromekelleher]

Looking again, I'm not sure my last comment is correct either! Here's the code that actually does the mapping. Any thoughts on what it's doing much appreciated!

[andrewkern]

heh. so it looks like this is computing something like the cumulative genetic distance to a given coordinate in physical space and then returning the genetic distance at that point (later to be discretized). I don't see anything here that would prevent us from allowing a zero value in the last interval... but it's Friday afternoon and maybe I'm missing something

[petrelharp]

@jeromekelleher was right the first time (as well as the docs): the last number is not used; to make sure it's not used, it's required to be zero. I guess that it's there because it is tidier if both positions and rates are the same length, somehow?

@andrewkern you're right; the code computes the cumulative genetic distance to x as a piecewise linear function, with positions and rates giving the nodes and the slopes.

[petrelharp]

While we're on the topic: I did previous thinking about this here: https://github.com/tskit-dev/msprime/issues/549 which resulted in these docs: https://msprime.readthedocs.io/en/stable/api.html#variable-recombination-rates

So, msprime does discrete recombination with possible breakpoints evenly spaced in genetic distance. Othe programs, e.g. SLiM, discretize recombination with evenly spaced breakpoints in physical distance. The latter is more natural, if you think about nucleotides anyhow. This means that there is no exact mapping between msprime and SLiM recombination maps. However, we still need to write one, so we can use the 'same' map in msprime and SLiM.

So, I think we need a function that takes a SLiM-style map (rates, breakpoints, and number of loci) and returns an msprime-style map that has the property that if mutations only occur at integer positions, then from the point of view of the resulting sequence the recombination map is as close as possible in some sense to the SLIM map. Said another way, if we simplify the tree sequence down to the history only at the integer positions, it should look just like it had the SLiM map. I think this is possible up to some error factor depending on machine tolerance and the minimum recombination rate.

[jeromekelleher]

@jeromekelleher was right the first time (as well as the docs): the last number is not used; to make sure it's not used, it's required to be zero. I guess that it's there because it is tidier if both positions and rates are the same length, somehow?

Thanks for confirming @petrelharp! So, this means that there's no point in allowing a non-zero rate at the end @andrewkern, as it would just invite misunderstandings and subtle bugs. I guess if you want recombination at the end you need to put in another entry just beyond the actual sequence length? It gets confusing here because msprime is working on a continuous space of [0, L) but I guess most maps are working on an integer space of [1, L]. Perhaps this is the issue?

[andrewkern]

Yes I think adding a bogus zero at the end is the simple way to go. It will be interesting to see how often this creeps up when folks start adding more and more recombination maps for the stdpopsim project.

[jeromekelleher]

OK, going to close this one for now. Can reopen if we revisit this.

[petrelharp]

We could change this to take a length n-1 vector of rates, but to allow in a deprecated way the length-n with a zero at the end.

[andrewkern]

Seems like a good idea to me

[jeromekelleher]

We could change this to take a length n-1 vector of rates, but to allow in a deprecated way the length-n with a zero at the end.

What do we gain from that though? We end up having the same issue when we read in text files --- you have to specify the end position somehow, and requiring an empty field at the end of a column is worse than requiring a zero.

[andrewkern]

Maybe I’m misinterpreting but think Peter is suggesting that the last zero would be added internally, blind to the user

[jeromekelleher]

I see where he's coming from, but I'm coming back to the text file case rather than inputting numpy arrays. How do we input two different length arrays?

position    rate
0               10
50             20
100

or

position    rate
0               10
50             20
100           0     # Required, but not used

I think the latter is better and less likely to either break things or end up in misinterpretation.

edit please squint and imaging the columns are lined up properly.

[petrelharp]

Oh, ok - it is from reading in text files. I hadn't figured out what was constraining the vectors to be the same length.

I mean, we could have the text file reading code take the zero, but remove the superfluous zero from this function. But, it is no big deal.

[jeromekelleher]

I mean, we could have the text file reading code take the zero, but remove the superfluous zero from this function. But, it is no big deal.

That's what's happening at the moment. Otherwise the last value is just being ignored in the rates list. I agree we should change this to be a length (n - 1) array when it's numpified.

[jradrion]

I believe found another issue with the handling of recombination maps: it appears that if both the last and second to last row also has a rate of 0.0, the final genomic window will not be created.

For example the maximum position simulated with the following map is 499869:
Chromosome StartPos Rate(cM/Mb)
2L 0 0.0000
2L 100000 4.0000
2L 200000 8.0000
2L 300000 6.0000
2L 400000 2.0000
2L 500000 0.0000

while the maximum position simulated with the following map is 399855:
Chromosome StartPos Rate(cM/Mb)
2L 0 0.0000
2L 100000 4.0000
2L 200000 8.0000
2L 300000 6.0000
2L 400000 0.0000
2L 500000 0.0000

[jradrion]

Also, with respect to the issue of window start positions and requiring the last window to have rate = 0.0, my suggestion would be to require each row to have both the start and end positions. This should remove any ambiguity.

[petrelharp]

I've edited your issue, could you check I wrote what you meant to say, @jradrion? And, could you please post the code you used, and the output that says this is a problem?

[jradrion]

Your edit is fine. The issue seems to occur only when both the last and second to last row have their rates set to zero.

The bit of code I used is as follows :

def main():
     parser = argparse.ArgumentParser()
     parser.add_argument('--map',dest='map',help='Input recombination map')
     parser.add_argument('--step',dest='step',help='Window step size: used for determining simulation parameters', type=int, default=1000)
     parser.add_argument('--subProjectName',dest='subProjectName',help='Name/prefix for simuation/experiment')
     parser.add_argument('--projectDir',dest='outDir',help='Directory for all project output. NOTE: the same projectDir must be used for all functions of vcfRho')
     parser.add_argument('--nCPU',dest='nCPU',help='Number of CPUs to use', type=int, default=1)
     args = parser.parse_args()



     mapSimDir=os.path.join(args.outDir,"mapSims")
     if not os.path.exists(mapSimDir):
         os.makedirs(mapSimDir)


     tsFILE=os.path.join(mapSimDir,os.path.basename(args.map).replace(".txt",".ts"))
     if not os.path.exists(tsFILE):
         print("Simulating recombination landscape...")
         recomb_map = msp.RecombinationMap.read_hapmap(args.map)
         ts=msp.simulate(
                 sample_size=50,
                 Ne=10**5,
                 mutation_rate=5*10**-9,
                 recombination_map=recomb_map)
         ts.dump(tsFILE)

     tsvcfFILE=os.path.join(mapSimDir,os.path.basename(args.map).replace(".txt",".vcf"))
     if not os.path.exists(tsvcfFILE):
         print("Writing ts to VCF...")
         ts=msp.load(tsFILE)
         with open(tsvcfFILE, "w") as fOUT:
             ts.write_vcf(fOUT,ploidy=2,contig_id="2L")


     ## Read in VCF
     inVCF=tsvcfFILE
     print("""Importing VCF: "%s"...""" %(inVCF))
     h5FILE=inVCF.replace(".vcf",".hdf5")


     ## Convert to hdf5 file if this has not been done previously
     flag=0
     if not os.path.exists(h5FILE):
         print("Converting VCF to HDF5...")
         allel.vcf_to_hdf5(inVCF,h5FILE,fields="*")
         flag=1


     ## Read in the hdf5
     print("""Importing HDF5: "%s"...\n""" %(h5FILE))
     callset=h5py.File(h5FILE, mode="r")
     var=allel.VariantChunkedTable(callset["variants"],names=["CHROM","POS"], index="POS")
     chroms=var["CHROM"]
     pos=var["POS"]
     print(pos.max())

if __name__ == "__main__":
     main()

This will print the last position in the VCF, which I assumed should fall somewhere within the last window of the recombination map. The maps I used are those from my first post.

[jradrion]

OK, hopefully this will be easier to debug with the following python script. It doesn't look like the issue is in the read_hapmap function, as this function produces all the windows expected. I have also attached two example recombination maps to test.

import sys
import msprime as msp
import allel
import h5py

##usage: python hapMapDebug.py testMap.txt

recomb_map = msp.RecombinationMap.read_hapmap(sys.argv[1])
print("Recombination map windows:", recomb_map.get_positions())

tsFILE="./tmp.ts"
vcfFILE="./tmp.vcf"
h5FILE="./tmp.hdf5"

ts=msp.simulate(
        sample_size=50,
        Ne=10**5,
        mutation_rate=5*10**-9,
        recombination_map=recomb_map)
ts.dump(tsFILE)

ts=msp.load(tsFILE)
with open(vcfFILE, "w") as fOUT:
    ts.write_vcf(fOUT,ploidy=2,contig_id="2L")


allel.vcf_to_hdf5(vcfFILE,h5FILE,fields="*",overwrite=True)


callset=h5py.File(h5FILE, mode="r")
var=allel.VariantChunkedTable(callset["variants"],names=["CHROM","POS"], index="POS")
chroms=var["CHROM"]
pos=var["POS"]
print("Max SNP position:", pos.max())

testMapA.txt
testMapB.txt

[jradrion]

Stripped further...

import sys
import msprime as msp

##usage: python hapMapDebug.py recombinationMap.txt

recomb_map = msp.RecombinationMap.read_hapmap(sys.argv[1])
print("Recombination_map.get_positions():", recomb_map.get_positions())

ts=msp.simulate(
        sample_size=50,
        Ne=10**5,
        mutation_rate=5*10**-9,
        recombination_map=recomb_map)

i=ts.get_num_sites()
print("Last SNP simulated:", ts.site(i-1))

[petrelharp]

Ok, this is a bug, and separate from this issue. I'm moving this to a new issue.

[jeromekelleher]

Is this issue still relevant @andrewkern, or can we close it?

[jeromekelleher]

Closing this due to inactivity