Dear community,
I have been using MAGeCK MLE for a couple of analyses previously and also recently.
As I understand the method, the maximum likelihood estimation iteratively takes turns in improving the estimation of efficacies of individual sgRNAs and the beta values, which are the main readout metric provided by the method.
Surprisingly, my sgRNA efficacies, as indicated in the .sgrna_summary.txt file, always end up at a value of "1", and if I recall correctly, this has been the case in any MAGeCK MLE analysis I have performed until now, leaving me three explanations:
- I am always committing the same mistake when setting up the analysis
- The MAGeCK package is buggy
- The MLE is expected to drive the efficacies to 1, but nevertheless works in estimating optimal beta values on the way.
I am pasting below the last command that I was using:
mageck mle \
--threads {threads} \
--count-table {input.counts} \
--design-matrix {input.design_matrix} \
--output-prefix {params.output_prefix} \
--norm-method median \ # gene-level normalization leads to same issue
--permutation-round 5
Also, I should note that I get occasional warnings in the the mageck mle log:
python3.6/site-packages/mageck/mleem.py:73: RuntimeWarning: divide by zero encountered in true_divide nb_r=np.divide(mu_sq,var_vec-mu_estimate)
python3.6/site-packages/scipy/stats/_discrete_distns.py:271: RuntimeWarning: invalid v alue encountered in subtract coeff = gamln(n+x) - gamln(x+1) - gamln(n)
python3.6/site-packages/scipy/stats/_discrete_distns.py:272: RuntimeWarning: invalid v alue encountered in multiply return coeff + n*log(p) + special.xlog1py(x, -p)
python3.6/site-packages/mageck/mleem.py:70: RuntimeWarning: overflow encountered in multiply mu_sq=np.multiply(mu_estimate,mu_estimate)
python3.6/site-packages/mageck/mleem.py:71: RuntimeWarning: overflow encountered in multiply var_vec=mu_estimate+np.multiply(alpha, mu_sq)
python3.6/site-packages/mageck/mleem.py:73: RuntimeWarning: invalid value encountered in true_divide nb_r=np.divide(mu_sq,var_vec-mu_estimate)
If you have any insights for me on the topic, I would be very thankful. Moritz
Thank you for the response. That's relieving to hear.
Actually the --update-efficiency switch is the way to go. Thanks for pointing me in the right direction! I really wonder why it is not enabled by default...