t-neumann / slamdunk

Streamlining SLAM-seq analysis with ultra-high sensitivity
GNU Affero General Public License v3.0
38 stars 23 forks source link

How to obtain half-life table? #82

Closed osamnsmr closed 4 years ago

osamnsmr commented 4 years ago

Hi, I want to get a table of half-life per gene from the slamdunk results. And I found an R script under the git plot folder that looks appropriate. slamdunk/plot/compute_halflifes.R

However, the number and content of the columns defined in this script are different from the results of the current version of slamdunk.

Therefore, I was hoping that the script would work with the following processing,

  1. Remove 3 header lines tail -n +4 filtered_tcount.cut.tsv > filtered_tcount.no-head.tsv

  2. Extract the columns according to the script definition cut -f 1-4,6-10,12-14 filtered_tcount.no-head.tsv > filtered_tcount.no-head.cut.tsv

but it stops with an error.

Error in summary(fit) : object 'fit' not found
Calls: computeHalfLife -> summary
Execution halted

This is an error when I remove the tryCatch().

Error in nlsModel(formula, mf, start, wts) :
  singular gradient matrix at initial parameter estimates
Calls: computeHalfLife -> nls -> nlsModel
Execution halted

Is there a better way?

t-neumann commented 4 years ago

Hi,

we never fully implemented the half-life estimation in slamdunk, that script is part of legacy code and should actually be removed.

What is your setup, do you have time-course data or a single pulse timepoint? Depending on this I could give you directions.

osamnsmr commented 4 years ago

Hi Tobias, Thank you for your quick reply.

I have time-course data of 0,1,3,6,12,24 hrs. The analysis I'm imagining is the suppl table in the SLAMseq paper, and that's exactly what compute_halflifes.R was perfect for.

https://static-content.springer.com/esm/art%3A10.1038%2Fnmeth.4435/MediaObjects/41592_2017_BFnmeth4435_MOESM4_ESM.xls

Nishi

osamnsmr commented 4 years ago

Sorry, the information above was insufficient.

[S4U pulse] > 0,1,3,6,12,24, [Uridine chase] > 1,3,6,12,24 hrs.

Nishi

t-neumann commented 4 years ago

Hi Nishi,

hm so we never used the scripts actually for calculating the half-lifes from the paper. What I remember is that they used the chase timepoints and the ConversionRate columns for fitting a linear curve on the logarithmic values because they have then also a proper endpoint (0).

Can you work with that info?

osamnsmr commented 4 years ago

Hi Tobias,

The script worked with slight modifications based on the content of this paper. https://www.nature.com/articles/s41592-020-0935-4

fit = nls(rates ~ I(a*exp(-b*timepoints)), start=list(a=1, b=0), control=nls.control(maxiter=1000), na.action=na.exclude)

Thank you so much for your help!

Nishi

t-neumann commented 4 years ago

Great to hear!

Cheers,

Tobi

ZoeyYang912 commented 11 months ago

Hi,

Would you mind sharing the code for getting the half-life table? I am new to this analysis and I don't how to use the slamdunk/plot/compute_halflifes.R under the git plot folder and set up the analysis. Thanks a lot!

t-neumann commented 11 months ago

Hi @ZoeyYang912 - well the R script is basically the code, since it's a simple single-exponential fit.

It should in principle work if you run it as Rscript slamdunk/plot/compute_halflifes.R -f slamdunktable1,slamdunktable2,slamdunktable3 -t timepoint1,timepoint2,timepoint3 (all in mins) -o output.txt