Thanks for watching the channel! t.co/XgH1fe3TE9
In this video, I show you how to obtain mRNA half-life values from SLAM-seq data. We cover how to obtain SLAM-seq data, download the genomic sequences and annotation. We then obtain conversion values using the Slamdunk pipeline. After combining all conversion values from different samples, we use this data to fit an exponential decay model that will allow us to estimate the half-life of any mRNA.
Chapters:
00:00 Intro
00:14 Downloading SLAM-seq data
01:10 Downloading fastq genomic reference
01:38 Downloading the annotation data
02:54 Downloading Slumdunk
03:16 Slumdunk metadata
03:47 Running Slumdunk
04:05 Conversion rate results
04:38 Downloading data from Biomart
05:05 Merging conversion rate files
06:23 Fitting the exponential decay model
07:32 Estimating half-life values
08:40 Outro
Links
Link to SLAM-seq manuscript
www.nature.com/articles/nmeth.4435
Link to UCSC genome browser
genome.ucsc.edu/
Link to Slamdunk repo
github.com/t-neumann/slamdunk
Code
Line of code to run Slumdunk with the data obtained.
Code
library(dplyr)
library(tidyr)
library(purrr)
library(ggplot2)
library(biomaRt)
Code to get data from Biomart
ensembl = useEnsembl(biomart="ensembl")
listDatasets(ensembl)
mart= useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
listAttributes(mart)
symbols_mRNA = getBM(attributes=c("mgi_symbol", "ensembl_transcript_id"), mart=mart)
write.table(symbols_mRNA, "../data/symbols_mRNA.csv", sep=",",
row.names = F, col.names = F, quote = F)
Code to obtain a combined table with conversion values
metadata = read.csv2("../data/sampleInfo_final.tsv",
header = T, sep = "\t", stringsAsFactors = F) %%
rename(sample = name) %% rowwise() %%
mutate(file = paste("../results/slamdunk/count/utrs/", sample, "_trimmed.fq_slamdunk_mapped_filtered_tcount.tsv", sep = "")) %%
select(file, sample, condition, time)
conversion = data.frame(file=character(), mRNA=character(),
conversion=numeric())
for(i in c(1:nrow(metadata))){
file = read.table(as.character(metadata[i,"file"]), stringsAsFactors = F,
skip = 1, header = T, sep="\t") %%
dplyr::select(Name, ConversionRate) %% rowwise() %%
rename(conversion = ConversionRate) %%
mutate(mRNA = unlist(strsplit(Name, "_utr3_"))[1]) %% rowwise() %%
mutate(mRNA = unlist(strsplit(mRNA, "\\."))[1]) %%
dplyr::select(-Name) %% unique() %%
mutate(file = as.character(metadata[i,"file"]))
conversion = bind_rows(conversion, file)
full_table = conversion %% left_join(metadata) %% dplyr::select(-file)
symbols_mRNA = read.csv2("../data/symbols_mRNA.csv",
sep=",", header = F, stringsAsFactors = F) %%
rename(symbol=V1, mRNA=V2)
full_table_symbols = full_table %% left_join(symbols_mRNA) %%
dplyr::select(-mRNA) %% unique() %% filter(!is.na(symbol)) %%
group_by(sample, condition, time, symbol) %%
mutate(conversion = max(conversion)) %% unique() %% ungroup()
normalized = full_table_symbols %% group_by(symbol, condition, time) %%
mutate(mean = mean(conversion)) %% ungroup() %%
group_by(symbol, condition) %% mutate(max = max(mean)) %% ungroup() %%
filter(max != 0) %% rowwise() %% mutate(nomalized = conversion/max) %%
dplyr::select(-max)
write.table(normalized, "../views/normalized.csv")
saveRDS(normalized, "../views/normalized.rds")
Code to fit the exponential decay model and estimate half-life values
normalized = readRDS("../views/normalized.rds")
nested_normalized = normalized %%
dplyr::select(condition, time, symbol, nomalized) %%
group_by(symbol, condition) %% nest()
test = function(df){
nls(nomalized ~ exp(-k * time), start = list(k = 0.01), data=df)}
test_possibly = possibly(test, NULL)
nested_normalized = nested_normalized %%
mutate(z = map(data, test_possibly))
hl_fun = function(z){
tryCatch({hl = log(2)/summary(z)$coefficients[1, 1]
return(hl)},
error=function(error_message){return(NA)})}
R_fun = function(z){
tryCatch({R = z$m$deviance()
return(R)},
error=function(error_message){return(NA)})}
nested_normalized_ = nested_normalized %%
rowwise() %% mutate(hl = hl_fun(z), R = R_fun(z)) %% dplyr::select(-data, -z)