Prepare input genetic map file

[XiaoleiZ]

Hi,

This is an awesome work! I am applying FastSMC in UKBB. For the input - genetic map file, what would you recommend one to use/generate? The described format looks like different from standard plink map output and public reference genetic map file (no chromosome id but with mutation rate).

Thanks in advance,
Xiaolei

[fcooper8472]

Hi @XiaoleiZ, thanks for getting in touch!

You are correct that it is not a standard map file (although this will likely change in the future).

For now, the genetic map file needs to provide physical positions (in base pairs) and genetic positions (in centimorgans). FastSMC expects a tab separated file containing at least 3 columns, with physical positions in the first column and genetic positions in the third column. The second column may contain any values. The file may contain a header row, or not, and it may be optionally compressed using gzip.

I have updated the README to make this specification clearer:
https://github.com/PalamaraLab/FastSMC#genetic-map-map

Please let me know if you have any other questions.

[XiaoleiZ]

Thanks I tried to generate the file as described.

However, I came across this 'std::bad_alloc' error. Google tells me that it is likely to be out of memory. But there was enough memory for the program. I will paste my command and the error message below. Any thoughts on how to solve it? Many thanks.

`
Exited with exit code 134.

Resource usage summary:

CPU time : 4269.78 sec.
Max Memory : 2702 MB
Average Memory : 1487.34 MB
Total Requested Memory : 80000.00 MB
Delta Memory : 77298.00 MB
Max Swap : 41 MB
Max Processes : 4
Max Threads : 5
Run time : 4317 sec.
Turnaround time : 4325 sec.

The output (if any) follows:

---

   ASMC OPTIONS       

---

Input will have prefix : /nfs/research/birney/users/xiaolei/phase_ukbb/phased_chr21
Decoding quantities file : /hps/software/users/birney/xiaolei/fast_smc_sh/FILES/DECODING_QUANTITIES/30-100-2000.decodingQuantities.gz
Output will have prefix : /nfs/research/birney/users/xiaolei/ibd_ukbb/chr21.1.1.gasmc.bibd
Binary output ? true
Time threshold to define IBD in generations : 50
Use batches ? true
Batch size : 32
Running job 1 of 1
Recall level 3
skipCSFSdistance is 0
compress ? false
useAncestral ? false
outputIbdSegmentLength ? true
doPerPairPosteriorMean ? true
doPerPairMAP ? true
noConditionalAgeEstimates ? true
Use hashing as a preprocessing step ? true

---

 hashing OPTIONS     

---

Minimum match length (in cM) : 1.5
Skipping words with (seeds/samples) less than 0
Minimum minor allele frequency : 0
Allowed gaps 1
Dynamic hash seed cutoff : 0
Read 488376 samples.
Read data for 976752 haploid samples and 10615 markers, 9 of which are monomorphic. This job will focus on 976752 haploid samples.
Using precomputed decoding info from /hps/software/users/birney/xiaolei/fast_smc_sh/FILES/DECODING_QUANTITIES/30-100-2000.decodingQuantities.gz
Will decode using AVX instruction set.

488376 sample identifiers read
*** loading word buffer 1 / 10
*** loading word buffer 2 / 10
*** loading word buffer 3 / 10
*** loading word buffer 4 / 10
*** loading word buffer 5 / 10
*** loading word buffer 6 / 10
*** loading word buffer 7 / 10
*** loading word buffer 8 / 10
*** loading word buffer 9 / 10
terminate called after throwing an instance of 'std::bad_alloc'
what(): std::bad_alloc
run_fastsmc.sh: line 19: 2094015 Aborted (core dumped) $FASTSMC/FastSMC_exe --inFileRoot $INPUT/phased_chr21 --outFileRoot $OUTPUT/chr21 --decodingQuantFile $DECODING_Q --mode array --time 50 --min_m 1.5 --segmentLength --hashing --perPairPosteriorMeans --perPairMAP --noConditionalAgeEstimates --bin`

[fcooper8472]

Hi @XiaoleiZ,

The std::bad_alloc here almost certainly indicates that FastSMC tried to allocate more memory than your system had available.
The information you have provided indicates FastSMC tried to allocate ~80GB:

Total Requested Memory : 80000.00 MB

This is in line with your description, if you're running FastSMC on the whole UK Biobank. From the program parameters it looks as though you are currently using a single job:

Running job 1 of 1

You may want to run 100s if not 1000s of jobs for each chromosome or chromosome arm to reduce the resources required for each job.

You can find an example here of running multiple jobs:
https://github.com/PalamaraLab/FastSMC/blob/master/cpp_example/FastSMC_example_multiple_jobs.sh

You will not want to run every job in parallel, of course, as that would not reduce your total memory footprint, but hopefully that example script demonstrates how to run FastSMC with multiple jobs.

Please let me know if you run in to any more problems.

[XiaoleiZ]

Thanks for the help. It works now.

Another question, I am confused for the output: what's the usage for "Column 11. (optional) Average mean posterior age estimate of the IBD segment"? I understand one shall use Column 12 for age estimate of IBD. But these two look very different in their scales. In the output files, for IBDs with the same length, mean posterior age estimate (Column 11) show different while MAP age estimates (Column 12) are always the same. What's the interpretation for mean posterior age estimate here?

Thanks very much in advance.

[jnaitsaada]

Hi @XiaoleiZ !

For a given pair of individuals and a given genomic site, FastSMC computes the posterior of the TMRCA in discretised time intervals, which are provided by the DecodingQuantities file (e.g probability of TMRCA between 0 and 5 generations, between 5 and 10, between 10 and 20 etc...).

For a given site, FastSMC computes two different age estimates:

1. the MAP age estimate at this specific site (i.e the center of the time interval that has highest probability)
2. the mean posterior age estimate at this specific site (i.e the weighted average time interval, where each time interval is weighted by its probability)

For a given IBD segments (made of multiple sites), FastSMC outputs a single number by computing averages across sites. So FastSMC outputs:

1. the average MAP age estimate along the segment (column 12)
2. the average mean posterior age estimate along the segment (column 11)

The MAP is a discrete random variable since it can only be one of the possible time intervals (e.g 2.5 for the first time interval, 7.5 for the second one, 15 for the third one etc...), while the mean posterior is a continuous random variable as it is a weighted average. This is why the MAP age estimates seem to be "always the same" while the mean posterior estimate show "different" values for different IBD segments.

There is no clear answer to which estimate you need to use, this is entirely up to you. In the FastSMC paper, we decided to consider the MAP age estimate. We believe for recent time scale (~ <200 generations) the MAP is a better estimate, while the mean posterior is better for deep time scale (~ >200 generations).

A couple other notes are that because the age estimate is bounded by the smallest possible MAP specified in the decoding quantities, for very long segments (e.g. >10 cM) you might want to use the length-based estimator we used in the paper. Also note that the age estimates are conditioned on the segment actually being IBD (i.e. a common ancestor below the specified time cutoff). You might therefore also consider the value of the IBD score when processing the output (i.e. don't trust segments with low IBD score, depending on the application). If you want to get age estimates that are not conditioned on the segment being IBD within the specified number of generations you can try using the --noConditionalAgeEstimates flag.