-
Notifications
You must be signed in to change notification settings - Fork 8
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
mcorr-bam - Very large phi_pool values from simulated data #10
Comments
Same problem for me. However, my data is based on the real experiment with recombination of H. pylori. I used whole genome sequencing mapping to reference genome. I have 5 samples, one of them have phi_pool of value = 1.67191e+06. I ma not sure the reason why. |
Same problem for me. I use the mcorr to analyze the metagenomic reads using the custom pipeline (including align the reference genome against the read, filter the bam based on the criteria (NM < 5 & MapQ > 30) and sort the bam). The result from *_bootstrapping_report.txt" is as following: Best, |
Hi, Did you figure out the reason behind this? |
Hi folks, First off, apologies to @sid-krish and @nhuanmc for our slow response to your questions. Looking at the data provided by @sid-krish (thanks for that!), it looks like the issue here is the same as that in issue #16 -- the data is being poorly fit by our coalescent model with recombination, which is resulting in physically unrealistic values for phi_pool. For experimental data like @nhuanmc and @mingleiR, this means that there's either a lack of recombination or the data is too noisy to determine if recombination has taken place. I would recommend re-fitting it using "mcorrFitCompare" (see our ReadMe subsection 2-ii). This will fit the data both with the model with recombination and a model without recombination, and you can determine using AIC values whether or not your data is best fit with the recombination model (my guess is it's best fit by a coalescent model without recombination). For more details, you can see the methods section of our new paper: https://elifesciences.org/articles/78533 For the simulations used by @sid-krish it looks like there just hasn't been sufficient time for enough diversity to build up for you to detect recombination in these simulation data. If you either increase the mutation rates, the simulation time, or the population size, you should see correlation profiles which are well fit by our recombination model. Lastly, regarding @sid-krish about the population recombination rate (rho) vs phi_pool, you can find all the definitions of these terms in the supplement of our lab's nature methods paper describing mcorr (see here). phi_pool is: phi_pool = 2Tgamma, where T is the average coalescence time of a pair of sequences in the pool, and gamma is the recombination rate per basepair per generation. So it corresponds to the number of recombination events which have occurred since coalescence of the pool. Happy to answer additional questions folks might have about this, and I hope that we can respond quicker to future queries! Best, |
Hi Asher, Thanks for your detailed explantion and helpful suggestion. I've checked the results of "mcorrFitCompare", and indeed found the AIC values of zero recombination model is a bit smaller than that of the recombination model. Because I have no prior knowledge about the recombination of the target species in the question, and not sure about my analytic pipeline. I'll scrutinize the pipeline using the published metagenomic datasets (in your work or other work citing your article) to validate it. Thanks again. Best, |
Hi Minglei, Sure thing! That sounds like a good plan. Keep us posted, I'm happy to help out in any way! Best, |
Hi Asher, I re-analyze the gut microbiome metagenome in the paper (Lin et al, 2019), and found the statistic, such as gamma/mu is different from the result shown in the supplmentary (Table S6). The phi_pool value is relative smaller than my environmental samples (and not very large), but its bootstrapping mean and sd value is large, shown as below. " (2) aligning the reads, merge the bam, and filter the alignment. (3) estimating the relative recombination rates Hope you could point the potential problem in the procedure above, which caused the difference in the summary statistics between both analyes. Best, |
Hi Minglei, Apologies for the slow reply. Regarding the bootstrapping, the reason that the mean and standard deviation are so high has to do with the boostrapping procedure. The way it works is it creates bootstrap replicates by randomly sampling with replacement from all the genes analyzed, then creates an average correlation profile over all those genes. Occasionally, one of the replicates is composed of many genes with very noisy correlation profiles. So what ends up happening is the model fits them poorly and yields physically unrealistic values for those parameters. Therefore, we typically take the median and the 95% CI from the bootstraps, but not the mean and standard deviation, which are often overweighted by these poor fits. In fact in the data you shared, both the median and bootstrap CI are a reasonable match for the parameter inferred from the actual data. Regarding your pipeline, I have a few suggestions. It looks like from Mingzhi paper that he did two analyses from that patient (there are two rows listed as patient 10131). One was using this dataset and the other was using this one. The main differences I see in what you've done and what Mingzhi did to obtain his data are that it does not appear that he trimmed the reads for the infant data, I think this was just for the ancient samples. Then for step 2 (aligning the reads, merge the bam, and filter the alignment), he used the following steps:
$db references the PanPhlAn database used (in this case it's for E coli). bamfilter is an in-house script which you can find here. bedtools can be found here. After this step, the reads were mapped with smalt using the following commands:
I am not sure where he specified mapping quality. Step 3 looks to be the same. I hope this helps and happy to keep trying to troubleshoot. Best, |
I am trying to analyse some simulated datasets but am getting very large results for phi_pool (i.e. 1e50).
I have some simulated datasets that I generated using the program fastsimbac that I would like to try with mcorr-bam. With the simulated datasets only the bam and reference fasta is available. For the gff I tried to simply point to the entire reference sequence that was used for alignment.
I have attached an example simulated bam, fasta and gff file where the phi_pool value is really large. The bam has been simulated with a population recombination rate (rho) 0.5, 50 genomes of size 10Kbp and population mutation rate (theta) of 0.1.
I used the following commands for mcorr-bam and mcorr-fit and have attached the bootstrapping_report file along with the other files generated by mcorr.
mcorr-bam test.gff Aligned.bam test_bam
mcorr-fit test_bam.csv test_bam
I was also wondering how phi_pool differs from the population recombination rate (rho)?
data.zip
The text was updated successfully, but these errors were encountered: