lh3 / psmc

Implementation of the Pairwise Sequentially Markovian Coalescent (PSMC) model
Other
146 stars 60 forks source link

How to check the recombination event? #45

Open DuttaAnik opened 2 years ago

DuttaAnik commented 2 years ago

Hello, Thank you for this super useful tool. I was able to run it successfully. However, I have a question regarding choosing the optimum -tvalue and -p model. I was reading in the manual where it is stated that "the -p and -t options are manually chosen such that after 20 rounds of iterations, at least ~10 recombinations are inferred to occur in the intervals each parameter spans."

My question is, how can I check whether ~10 recombinations are inferred or not? I checked the output file after psmc run but could not find it. Maybe I am looking at the wrong place.

It would really help me to choose the correct -p model if you can tell me where to find that information. Thank you.

DuttaAnik commented 1 year ago

Hello, can someone give me an answer to this question above? It really matters to me at the moment. Thank you.

Buffan3369 commented 1 year ago

Hi Dutta,

Not sure you still need it but I'm writing it down here in case some people are also looking for the same thing. To be honest, I was also faced with the same issue, and I found the total absence of practical indications quite frustrating. While working hard on the theoretical basis of the PSMC model (available in the supp. info of the Li and Durbin 2011 paper), I read this between pages 21 and 22: image

On the page 20 of the same document, they come to establish the relation $C{\sigma} \approx \cfrac{1}{\rho C{\pi}}$

If you now have a look at the output .psmc file, you might see the things I indicated below. image

Hence, by taking the parameters (rho, C_pi and pi_k) returned at the last iteration you previously specified (-N argument), you might be able to compute what you are looking for ! :)

Hope it helped

dareen22 commented 1 year ago

This is very helpful, thank you @Buffan3369 for sharing!

DuttaAnik commented 1 year ago

Indeed, this is awesome. Thanks a lot, @Buffan3369. It worked finally!

rapamycin commented 9 months ago

Thank you so much ! @Buffan3369 Does it means that π_K/(ρC_π) should larger than 20, or π_K should larger than 20ρ*C_π?