Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change in time segments with masking #46

Open
EEPuckett opened this issue Aug 31, 2022 · 3 comments
Open

Change in time segments with masking #46

EEPuckett opened this issue Aug 31, 2022 · 3 comments

Comments

@EEPuckett
Copy link

Dear Dr. Schiffels and other msmc users,

I'm trying to better understand how msmc works. The figure below is a single sample (2 haplotypes; phased) and I've run through msmc2. The black line is the result using only the mask from bamCaller. The purple line has an additional mask on the genes and flank regions; the pink line has the mappability mask from SNPable added. I used the default time segments.

The results with the additional masking push the left time boundary far beyond the speciation age of my study organism. I'm trying to determine if I have a problem with this analysis as I did not expect either masking file to change the results in this way.

I appreciate any insight.
Emily

psmcprime

@stschiff
Copy link
Owner

stschiff commented Sep 5, 2022

OK, this is complicated. MSMC determines the time-boundaries from the average heterozygosity of the input data. Of course, with different masks you may get different thetas, so time will be scaled differently. Not sure if I read the Log-scale correctly, but it looks like you get a factor of 2 or so between the end point of the different curves, which is not really expected. You can check the different thetas in the log file, and perhaps confirm with an external tool from the input files.

Basically you need to investigate why these masks result in different thetas.

@EEPuckett
Copy link
Author

I think I have figured this out. While there was a problem with my mappability mask (related to only running the first segment coming out of splitfa for snpable; thus, only having a positive mask for a fraction of each chr), the real issue was with phasing.

My goal is to run msmc to get input for msmc-IM between two species that have signals of introgression. When I phased all of the samples (beagle v5) from both species together, that moves the time bounds to the right (ie- back in time). However, when I phase the species separately, I see the time bounds as I expect based on published phylogenies.

@stschiff
Copy link
Owner

Ah OK. So incorrect phasing may have moved mutations in biased ways onto haplotypes, so that some haplotypes were much more similar than expected. Interesting, but not too surprising given the complexity of statistical phasing

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants