Calculation of DRC statistics

[kapoormanav]

Hi,

I have run the FastSMC with T = 50. I want to calculate the DRC statistics for my dataset. Is there a script to calculate this statistics from the FastSMC output?

Thanks a lot!

[jnaitsaada]

Hi kapoormanav,
Sorry, we don’t provide scripts for these analyses because they require special care to details that are specific to each analyzed dataset (e.g. population size, inflation, ...). We have added a few suggestions in the ASMC manual, which I paste below (we will also add a comment on IBD in the FastSMC docs, thanks for pointing out):

---

The DRC statistic can be obtained by summing entries of the sumOverPairs matrices that correspond to the desired time interval, and averaging across SNPs within a window. Please refer to the ASMC and FastSMC papers for details.

A few notes on analyzing the DRC statistic: it is important to remember that using a Gamma distribution as a null model for the DRC statistic is an approximation (see paper), and further work is needed to derive an exact null model. In particular the DRC cannot be larger than 1, so approximate p-values obtained under a Gamma model for very large DRC values will be artificially small. We found the Gamma approximation to be reasonable under some conditions, such as a large recent effective population size, which makes the DRC small in neutral regions, but this approximation will not work well in all populations (e.g. isolated populations) or for large time intervals. If you decide to use this approach to detect candidate regions for selection, we recommend carefully checking that the Gamma approximation is reasonable. Because particularly large DRC values may lead to artificially small p-values, you may also want to use the Gamma approximation only to determine an approximate significance threshold, then report raw DRC values larger than this threshold, rather than reporting approximate p-values. Finally, we also suggest that you check that candidate regions detected this way do not overlap problematic regions of the genome, such as regions of very high/low recombination rate or LD, as well as large structural variants, such as inversions, which may affect recombination.

[kapoormanav]

Hi Thanks @jnaitsaada!

The FastSMC manuscript states, "FastSMC does not output the posterior of the TMRCA but provides an IBD quality score, corresponding to the sum of posterior probabilities between generations 0 and T, where T is the user-specified threshold. As the UK Biobank dataset was analyzed for T = 50, the DRC50 statistic at a given site along the genome was estimated by averaging all IBD quality scores obtained from all analyzed pairs of samples (assuming a score of 0 if no segment is present for a pair)."

So should I take the column number 11 (IBD score) from FastSMC output and average it across all pairs of individuals for given segment?

[jnaitsaada]

Hi @kapoormanav

I think what you describe is correct, although I am not entirely sure what you mean by "for a given segment", you need to do what you described for a given genomic site. When you average all the scores, do not forget to account for the pairs of individuals not sharing any IBD segment on the genomic site you are considering (IBD score = 0 for those pairs).
e.g: for site number 1, you need to consider all the IBD segments in the FastSMC output that contain site number 1, sum all the IBD scores of those segments, and divide by the total number of haplotype pairs in your dataset. Then you iterate for site number 2 and so on...

[kapoormanav]

Hi @jnaitsaada,

Thanks a lot for suggestions. I have summed the quality scores for the sites that are 0.3 cM apart. For averaging should I divide it by N pairs (sample size x 2) or N choose 2 (all possible combinations of pairs)?

Also I am seeing extraordinary large values for IBD score at chromosome 15 specifically. Chromosome 15 values are at least 10-20 times larger than other chromosomes. Can you please suggest some probable explanation for this? Can duplications at this chromosome cause larger values?

Manav

[jnaitsaada]

Hi @kapoormanav,

You should divide by N choose 2, where N is the sample size in haploids.

It's hard to tell without looking at the data but as we mentioned in the caveats listed above, there are several things to be careful about. These might include issues with the genetic map, overlap with centromeres (best to analyze chromosome arms separately), presence of structural variant (e.g. inversions, duplications). If in doubt it's probably best to exclude the region.