Corrigendum: Antagonistic actions of two human Pan3 isoforms on global mRNA turnover

Chyi-Ying A. Chen; Yueqiang Zhang; Yu Xiang; Leng Han; Ann-Bin Shyu

doi:10.1261/rna.068387.118

Corrigendum: Antagonistic actions of two human Pan3 isoforms on global mRNA turnover

RNA 23: 1404–1418 (2017)

We are submitting this corrigendum to correct a few minor errors and add more detail about the bioinformatics method for mRNA half-life analysis in this paper. We are also adding Dr. Jeffrey T. Chang, who performed the related bioinformatics analysis for panels C and D of Figure 7 and is providing these corrections with additional information, as a coauthor (Dr. Chang was previously credited in the paper's “Acknowledgments” section). These changes do not affect the interpretations of the data and conclusions of the study. The changes are detailed as follows:

Page 1412, line 4 from the bottom (left): “…using circular plots generated with CIRCOS software (Hu et al. 2014).” The reference should instead be Krzywinski et al. 2009. Circos: an information aesthetic for comparative genomics. Genome Res 19:1639–1645.
Page 1416, under RNA-sequencing and decay rate calculation: TPM is described as transcripts per kilobase million but, instead, it should be noted as transcripts per million.
Page 1417, under RNA-sequencing and decay rate calculation: The statement “Transcripts whose half-life estimates could not be accurately determined due to a high residual value…in linear regression analysis were further filtered out” should be changed as follows: “We estimated the transcript half-life (t_1/2) by fitting the normalized values to the equation Y(t) = X(t)/X(0) = C*exp(−at) using a nonlinear least-squares estimation approach implemented in R.”
We now include the script and instruction for calculation of mRNA half-life (see below: § The R Script for calculating mRNA half-lives from RNA-seq). Briefly, for estimating decay rates from the normalized values from time-course and RNA-seq experiments, we characterized each transcript's decay characteristics by nonlinear least-squares fitting of the normalized values to the equation: Y(t) = X(t)/X(0) = C*exp(−at), where C (initial intensity) and a (decay rate constant) are floating parameters. This equation was modified from a first-order exponential decay model described in Imamachi et al. (2014), where the C term was fixed and set as 1. The C term in the equation we used provided more freedom for variation in the first time point, thus yielding better estimates of the half-lives. We also quantified the goodness of fit by evaluating the root-mean-square deviation (RMSD) between the experimental and fitted values.
We add Dr. Jeffrey T. Chang to the author list in the penultimate position. The updated list of authors and affiliations is:

Chyi-Ying A. Chen,¹ Yueqiang Zhang,^1,4 Yu Xiang,¹ Leng Han,¹ Jeffrey T. Chang,^2,3 and Ann-Bin Shyu¹

¹Department of Biochemistry and Molecular Biology, McGovern Medical School, The University of Texas Health Science Center at Houston, Houston, Texas 77030, USA

²Department of Integrative Biology and Pharmacology, McGovern Medical School, The University of Texas Health Science Center at Houston, Houston, Texas 77030, USA

³School of Biomedical Informatics, The University of Texas Health Science Center at Houston, Houston, Texas 77030, USA

⁴Current Address: Center for Cell and Gene Therapy, Baylor College of Medicine, Houston, Texas 77030, USA

Corresponding author: Ann-Bin.Shyu@uth.tmc.edu

§ The R script for calculating mRNA half-lives from RNA-seq:

This R script will take the pan3.tpm file (see Additional Supplemental Material link for this paper) containing relative abundance of transcripts measured for each gene from time-course experiments under control (control knockdown) and experimental (Pan3S, Pan3L, and Pan3S+L knockdowns) conditions, calculate half-life for each gene, and generate the three scatter plots as shown in Figure 7C of the paper. The specific steps of the R script are as follows: 1. Normalize each of the transcripts to the expression of GAPDH. 2. Calculate the half life using a first-order decay equation: Y(t) = X(t)/X(0) = C * e^{-at} A table containing all the half-lives calculated can be assembled after this step. 3. Filter the genes as follows: - Must be expressed with TPM >=10. - Must have estimated half live >0 hours and <48 hours. - Decay must fit first-order decay curve with residual <0.5. - Residual is defined as \sum (y^{hat}-y)^{2}, where y is the gene expression value, and y^{hat} is the value predicted by the equation in step 2. Three tables listing half-lives from qualified genes can be assembled after this step 4. Make scatter plots: The half lives of the control samples are plotted against the half lives of the experimental samples. An identity line is plotted to make clear changes in the half lives. The R script “calc_decay.R”: library(parallel) data <- read.delim(“pan3.tpm”, header=TRUE, as.is=TRUE, comment.char="", quote="") X <- as.matrix(data[,5:ncol(data)]) # Normalize to GAPDH (NM_002046). I <- which(data[[“RefSeq.ID"]] == “NM_002046”) if(length(I) != 1) stop(“can't find NM_002046”) gapdh <- as.numeric(X[I,]) X.norm <- t(apply(X, 1, function(x) x/gapdh)) # Time points. tm <- c(0, 0.5, 1, 2, 4, 6) SAMPLE.SETS <- list( C=list(I=1:6), L=list(I=7:12), LS=list(I=13:18), S=list(I=19:24) ) calc.decay.m1 <- function(y, tm) { control <- list(maxiter=50, tol=1e-5) a.start <- log(2)/5 # assume 5 hour half life C.start <- 1 start <- list(a=a.start, C=C.start) fit <- nls(y ∼ C*exp(-a*tm), start=start, control=control, trace=FALSE) a.hat <- coef(fit)[1] C.hat <- coef(fit)[2] y.hat <- C.hat * exp(-a.hat * tm) resid <- sum((y-y.hat)**2) T.half.hat <- log(2)/a.hat list(a.hat=a.hat, C.hat=C.hat, y.hat=y.hat, T.half.hat=T.half.hat, fit=fit, resid=resid) } HALF.LIFE <- list() for(j in 1:length(SAMPLE.SETS)) { I <- SAMPLE.SETS[[j]]$I X.analyze <- X.norm[,I] if(length(tm) != ncol(X.analyze)) stop(“unaligned”) start <- proc.time() x <- mclapply(1:nrow(X.analyze), mc.cores=12, FUN=function(i) { if(i %% 1000 == 0) print( sprintf(“[%d:%d] [%d:%d] Calculating half life %.2fs”, j, length(SAMPLE.SETS), i, nrow(X.analyze), (proc.time()-start)[3])) Y.t <- as.numeric(X.analyze[i,]) # Normalize to first time point. Y.t <- Y.t/Y.t[1] m1 <- tryCatch({ calc.decay.m1(Y.t, tm) }, error=function(err) { as.character(err) }) x <- rep("", 5) if(is.character(m1)) { x[5] <- gsub(“[\r\n]”, "", m1) } else { x <- c(m1$a.hat, m1$C.hat, m1$T.half.hat, m1$resid, "") } x }) data.out <- matrix(unlist(x), nrow=length(x), byrow=TRUE) colnames(data.out) <- c("M1 a", "M1 C", "M1 T_half", "M1 residual", "M1 err") rownames(data.out) <- rownames(data) data.out <- cbind(data[,1:4], data.out) HALF.LIFE[[names(SAMPLE.SETS)[j]]] <- data.out } # Filter the data. F.HALF.LIFE <- list() for(i in 2:length(HALF.LIFE)) { # Time 0 must have TPM >= 10 in both samples. I1 <- SAMPLE.SETS[[1]]$I I2 <- SAMPLE.SETS[[i]]$I t0.1 <- X[,I1[1]] t0.2 <- X[,I2[1]] F1 <- (t0.1 >= 10) & (t0.2 >= 10) # 0 < T1/2 < 48 t12.1 <- HALF.LIFE[[1]][[“M1 T_half”]] t12.2 <- HALF.LIFE[[i]][[“M1 T_half”]] t12.1 <- as.numeric(levels(t12.1))[t12.1] t12.2 <- as.numeric(levels(t12.2))[t12.2] I1 <- (t12.1 > 0) & (t12.1 < 48) I2 <- (t12.2 > 0) & (t12.2 < 48) F2 <- I1 & I2 F2[is.na(F2)] <- FALSE # residual < 0.5 r.1 <- HALF.LIFE[[1]][[“M1 residual”]] r.2 <- HALF.LIFE[[i]][[“M1 residual”]] r.1 <- as.numeric(levels(r.1))[r.1] r.2 <- as.numeric(levels(r.2))[r.2] F3 <- (r.1 < 0.5) & (r.2 < 0.5) F3[is.na(F3)] <- FALSE I <- F1 & F2 & F3 F.HALF.LIFE[[names(SAMPLE.SETS)[i]]] <- HALF.LIFE[[i]][I,] } # Make scatter plot. matrix2color <- function(cmatrix, pos) { # pos is [0, 1]. Returns r, g, b where each one is from [0, 1]. if(is.nan(pos)) # this can happen if someone calls matlab.colors(1) pos <- 0.5 breaks <- cmatrix[,1] i1 <- sum(pos >= breaks) x <- cmatrix[i1,2:4] if(i1 < nrow(cmatrix)) { i2 <- i1 + 1 delta <- (pos - cmatrix[i1,1]) / (cmatrix[i2,1]-cmatrix[i1,1]) x <- cmatrix[i1,2:4] + delta*(cmatrix[i2,2:4]-cmatrix[i1,2:4]) } x/255 } matrix2rgb <- function(cmatrix, pos) { # pos is [0, 1]. Returns color, e.g. “#003300”. x <- matrix2color(cmatrix, pos) rgb(x[1], x[2], x[3], maxColorValue=1) } bild.colors <- function(n) { if(n <= 0) stop(“n must be greater than 0”) cmatrix <- t(matrix(c( c(0.000, 49, 50, 114), c(0.050, 61, 69, 137), c(0.100, 62, 84, 154), c(0.150, 67, 89, 160), c(0.200, 85, 108, 176), c(0.250, 115, 145, 201), c(0.300, 160, 205, 240), c(0.350, 180, 220, 243), c(0.400, 169, 216, 211), c(0.450, 160, 208, 164), c(0.500, 179, 213, 112), c(0.550, 203, 220, 61), c(0.600, 232, 231, 61), c(0.650, 255, 234, 47), c(0.700, 250, 180, 50), c(0.750, 243, 136, 54), c(0.800, 231, 80, 61), c(0.850, 218, 54, 55), c(0.900, 204, 55, 59), c(0.950, 160, 52, 52), c(1.000, 114, 39, 44)), 4, 21)) sapply((0:(n-1))/(n-1), function(x) matrix2rgb(cmatrix, x)) } for(i in 1:length(F.HALF.LIFE)) { SAMPLE <- names(F.HALF.LIFE)[i] D1 <- HALF.LIFE[[1]] # control D2 <- F.HALF.LIFE[[i]] # Align the control samples. I <- match(D2[[“RefSeq.ID"]], D1[["RefSeq.ID”]]) if(any(is.na(I))) stop(“unaligned”) D1 <- D1[I,] x <- D1[[“M1 T_half”]] y <- D2[[“M1 T_half”]] x <- as.numeric(levels(x))[x] y <- as.numeric(levels(y))[y] x[x<0] <- 48 x[x>48] <- 48 y[y<0] <- 48 y[y>48] <- 48 xlim <- c(0, 48) ylim <- c(0, 48) nbin <- 256 nrpoints <- length(x)*0.1 # 10% of points pch <- 19 cex <- 0.5 cex.lab <- 2 xlab <- “C T_half” ylab <- sprintf(“%s T_half”, SAMPLE) p <- bild.colors(25) p <- c(rep(“#FFFFFF”, 10), p) colramp <- colorRampPalette(p) OUTFILE <- sprintf(“%s.pdf”, SAMPLE) print(sprintf(“Writing file: %s”, OUTFILE)) bitmap(OUTFILE, type=“pdfwrite”, height=6, width=6, units=“in”, res=600) op <- par( mar=(c(5, 5*1.2, 4, 2)+0.1), mgp=c(3, 1.5*0.8, 0), bg=“#FFFFFF” ) smoothScatter(x, y, nbin=nbin, nrpoints=nrpoints, colramp=colramp, pch=pch, cex=cex, xlab=xlab, ylab=ylab, cex.lab=1.5, cex.axis=2, xlim=xlim, ylim=ylim) lines(c(0, 48), c(0, 48), lty=2, lwd=4, col=“#FF0000”) par(op) dev.off() } print(“Done”)

doi: 10.1261/rna.068387.118