2. PAT-Seq poly(A) tail length example
2022-04-26
Contents
poly(A) tail length of transcripts can be measured using the PAT-Seq protocol. This protocol produces 3’-end focussed reads that include the poly(A) tail. We examine GSE83162. This is a time-series experiment in which two strains of yeast are released into synchronized cell cycling and observed through two cycles. Yeast are treated with \(\alpha\)-factor, which causes them to stop dividing in antici… pation of a chance to mate. When the \(\alpha\)-factor is washed away, they resume cycling.
1 Read files, extract experimental design from sample names
library(tidyverse) # ggplot2, etc
library(patchwork) # side-by-side plots
library(limma) # differential testing
library(topconfects) # differential testing - top confident effect sizes
library(org.Sc.sgd.db) # Yeast organism info
library(weitrix) # Matrices with precisions weights
# Produce consistent results
set.seed(12345)
# BiocParallel supports multiple backends.
# If the default hangs or errors, try others.
# The most reliable backed is to use serial processing
BiocParallel::register( BiocParallel::SerialParam() )tail <- system.file("GSE83162", "tail.csv.gz", package="weitrix") %>%
read_csv() %>%
column_to_rownames("Feature") %>%
as.matrix()
tail_count <- system.file("GSE83162", "tail_count.csv.gz", package="weitrix") %>%
read_csv() %>%
column_to_rownames("Feature") %>%
as.matrix()
samples <- data.frame(sample=I(colnames(tail))) %>%
extract(sample, c("strain","time"), c("(.+)-(.+)"), remove=FALSE) %>%
mutate(
strain=factor(strain,unique(strain)),
time=factor(time,unique(time)))
rownames(samples) <- colnames(tail)
samples## sample strain time
## WT-tpre WT-tpre WT tpre
## WT-t0m WT-t0m WT t0m
## WT-t15m WT-t15m WT t15m
## WT-t30m WT-t30m WT t30m
## WT-t45m WT-t45m WT t45m
## WT-t60m WT-t60m WT t60m
## WT-t75m WT-t75m WT t75m
## WT-t90m WT-t90m WT t90m
## WT-t105m WT-t105m WT t105m
## WT-t120m WT-t120m WT t120m
## DeltaSet1-tpre DeltaSet1-tpre DeltaSet1 tpre
## DeltaSet1-t0m DeltaSet1-t0m DeltaSet1 t0m
## DeltaSet1-t15m DeltaSet1-t15m DeltaSet1 t15m
## DeltaSet1-t30m DeltaSet1-t30m DeltaSet1 t30m
## DeltaSet1-t45m DeltaSet1-t45m DeltaSet1 t45m
## DeltaSet1-t60m DeltaSet1-t60m DeltaSet1 t60m
## DeltaSet1-t75m DeltaSet1-t75m DeltaSet1 t75m
## DeltaSet1-t90m DeltaSet1-t90m DeltaSet1 t90m
## DeltaSet1-t105m DeltaSet1-t105m DeltaSet1 t105m
## DeltaSet1-t120m DeltaSet1-t120m DeltaSet1 t120m“tpre” is the cells in an unsynchronized state, other times are minutes after release into cycling.
The two strains are a wildtype and a strain with a mutated set1 gene.
2 Create weitrix object
These tail lengths are each the average over many reads. We therefore weight each tail length by the number of reads. This is somewhat overoptimistic as there is biological noise that doesn’t go away with more reads, which we will correct for in the next step.
good <- rowMeans(tail_count) >= 10
table(good)## good
## FALSE TRUE
## 1657 4875wei <- as_weitrix(
tail[good,,drop=FALSE],
weights=tail_count[good,,drop=FALSE])
rowData(wei)$gene <- AnnotationDbi::select(
org.Sc.sgd.db, keys=rownames(wei), columns=c("GENENAME"))$GENENAME
rowData(wei)$total_reads <- rowSums(weitrix_weights(wei))
colData(wei) <- cbind(colData(wei), samples)3 Calibration
Our first step is to calibrate our weights. Our weights are overoptimistic for large numbers of reads, as there is a biological components of noise that does not go away with more reads.
Calibration requires a model explaining non-random effects. We provide a design matrix and a weighted linear model fit is found for each row. The lack of replicates makes life difficult, for simplicity here we will assume time and strain are independent.
design <- model.matrix(~ strain + time, data=colData(wei))
fit <- weitrix_components(wei, design=design)A gamma GLM with log link function can then be fitted to the squared residuals. 1 over the predictions from this models will serve as the new weights.
cal <- weitrix_calibrate_all(
wei,
design = fit,
trend_formula = ~well_knotted_spline(mu,3)+well_knotted_spline(log(weight),3))## mu range 4.402217 77.000000 knots 23.04834 30.07095 37.91549## log(weight) range 0.0000 10.6244 knots 2.841452 4.397124 6.359334For comparison, we’ll also look at completely unweighted residuals.
unwei <- wei
weitrix_weights(unwei) <- weitrix_weights(unwei) > 0
# (missing data still needs weight 0)Calibration should remove any pattern in the weighted residuals, compared to known covariates. The trend line shown in red is based on the squared weighted residuals.
First look for any pattern relative to the linear model prediction (“mu”). A trend has been removed by the calibration.
weitrix_calplot(unwei, fit, covar=mu, guides=FALSE) + labs(title="Unweighted\n") |
weitrix_calplot(wei, fit, covar=mu, guides=FALSE) + labs(title="Weighted by\nread count") |
weitrix_calplot(cal, fit, covar=mu, guides=FALSE) + labs(title="Calibrated\n")
Next look for any pattern relative to the number of reads (recall these were the original weights). Again, a trend has been removed by the calibration.
weitrix_calplot(unwei, fit, covar=log(weitrix_weights(wei)), guides=FALSE) + labs(title="Unweighted\n") |
weitrix_calplot(wei, fit, covar=log(weitrix_weights(wei)), guides=FALSE) + labs(title="Weighted by\nread count") |
weitrix_calplot(cal, fit, covar=log(weitrix_weights(wei)), guides=FALSE) + labs(title="Calibrated\n")
4 Testing
4.1 Top confident effects
We are now ready to test things.
My recommended approach is to find top confident effects, here top confident differential tail length. Core functionality is implemented in my package topconfects. Applying this to a weitrix is done using the weitrix_confects function. Rather than picking “most significant” genes, it will highlight genes with a large effect size.
weitrix_confects(cal, design, coef="strainDeltaSet1", fdr=0.05)## $table
## confect effect se df fdr_zero row_mean typical_obs_err
## 1 -2.569 -6.027 0.6864 30.25 3.915e-06 26.23 1.485
## 2 -2.104 -7.124 1.0472 30.25 3.538e-04 29.08 2.057
## 3 -1.212 -5.070 0.8298 30.25 1.207e-03 14.87 1.773
## 4 -1.064 -3.420 0.5180 30.25 4.108e-04 19.10 1.119
## 5 -0.872 -3.470 0.5833 30.25 1.265e-03 20.29 1.267
## 6 -0.872 -7.505 1.5042 30.25 6.497e-03 38.19 3.062
## 7 -0.872 -3.304 0.5610 30.25 1.279e-03 18.36 1.218
## 8 -0.872 -6.710 1.3544 30.25 6.497e-03 20.48 2.627
## 9 -0.831 -2.867 0.4782 30.25 1.265e-03 26.88 1.055
## 10 0.828 5.472 1.0990 30.25 6.497e-03 29.08 2.290
## name gene total_reads
## 1 YDR170W-A <NA> 4832
## 2 YIL015W BAR1 3898
## 3 YAR009C <NA> 1866
## 4 YJR027W/YJR026W <NA> 6577
## 5 YML045W/YML045W-A <NA> 4553
## 6 YHR051W COX6 1167
## 7 YPL257W-B/YPL257W-A <NA> 5729
## 8 YOR192C-B/YOR192C-A <NA> 864
## 9 YIL053W GPP1 25996
## 10 YGL209W MIG2 1937
## ...
## 112 of 4875 non-zero contrast at FDR 0.05
## Prior df 21.2This lists the largest confident changes in poly(A) tail length.
The confect column is an inner confidence bound on
the difference in tail length,
adjusted for multiple testing.
Note that due to PCR amplification slippage and limited read length, the observed poly(A) tail lengths may be an underestimate. However as all samples have been prepared in the same way, observed differences should indicate the existence of true differences.
If you prefer to rank by signal to noise ratio, it is possible to use Cohen’s f as an effect size. This is similar to ranking by p-value, but Cohen’s f can be interpreted meaningfully as a signal to noise ratio.
weitrix_confects(cal, design, coef="strainDeltaSet1", effect="cohen_f", fdr=0.05)## $table
## confect effect strainDeltaSet1 fdr_zero row_mean typical_obs_err
## 1 0.569 1.963 -6.027 3.915e-06 26.23 1.4851
## 2 0.322 1.521 -7.124 3.538e-04 29.08 2.0570
## 3 0.317 1.476 -3.420 4.108e-04 19.10 1.1193
## 4 0.259 1.366 -5.070 1.207e-03 14.87 1.7733
## 5 0.259 1.341 -2.867 1.265e-03 26.88 1.0552
## 6 0.259 1.330 -3.470 1.265e-03 20.29 1.2672
## 7 0.259 1.317 -3.304 1.279e-03 18.36 1.2180
## 8 0.174 1.190 -2.846 5.584e-03 18.31 1.1618
## 9 0.171 1.170 -2.804 6.426e-03 18.99 1.1637
## 10 0.171 1.144 -2.283 6.497e-03 27.95 0.9934
## name gene total_reads
## 1 YDR170W-A <NA> 4832
## 2 YIL015W BAR1 3898
## 3 YJR027W/YJR026W <NA> 6577
## 4 YAR009C <NA> 1866
## 5 YIL053W GPP1 25996
## 6 YML045W/YML045W-A <NA> 4553
## 7 YPL257W-B/YPL257W-A <NA> 5729
## 8 YLR157C-B/YLR157C-A <NA> 5519
## 9 YPR137C-B/YPR137C-A <NA> 7287
## 10 YOL109W ZEO1 34839
## ...
## 112 of 4875 non-zero Cohen's f at FDR 0.05
## Prior df 21.24.2 Testing with limma
If you prefer a more traditional approach, we can feed our calibrated weitrix to limma.
fit_cal_design <- cal %>%
weitrix_elist() %>%
lmFit(design)
ebayes_fit <- eBayes(fit_cal_design)
topTable(ebayes_fit, "strainDeltaSet1", n=10) %>%
dplyr::select(
gene,diff_tail=logFC,ave_tail=AveExpr,adj.P.Val,total_reads)## gene diff_tail ave_tail adj.P.Val total_reads
## YDR170W-A <NA> -6.027208 25.87662 3.915262e-06 4832
## YJR027W/YJR026W <NA> -3.419778 19.40928 4.107986e-04 6577
## YIL015W BAR1 -7.124118 30.24503 3.538329e-04 3898
## YAR009C <NA> -5.070438 15.56624 1.207189e-03 1866
## YIL053W GPP1 -2.867498 27.19599 1.265338e-03 25996
## YML045W/YML045W-A <NA> -3.470070 20.55073 1.265338e-03 4553
## YPL257W-B/YPL257W-A <NA> -3.304377 18.59545 1.278512e-03 5729
## YLR157C-B/YLR157C-A <NA> -2.845888 18.46280 5.584245e-03 5519
## YPR137C-B/YPR137C-A <NA> -2.803791 19.21016 6.425669e-03 7287
## YOL109W ZEO1 -2.282940 28.11748 6.497358e-03 348394.3 Testing multiple contrasts
weitrix_confects can also be used as an omnibus test of multiple contrasts. Here the default effect size will be standard deviation of observations explained by the part of the model captured by the contrasts. The standardized effect size “Cohen’s f” can also be used.
Here we will look for any step changes between time steps, ignoring the “tpre” timepoint. The exact way these contrasts are specified will not modify the ranking, so long as they specify the same subspace of coefficients.
multiple_contrasts <- limma::makeContrasts(
timet15m-timet0m, timet30m-timet15m, timet45m-timet30m,
timet60m-timet45m, timet75m-timet60m, timet90m-timet75m,
timet105m-timet90m, timet120m-timet105m,
levels=design)## Warning in limma::makeContrasts(timet15m - timet0m, timet30m - timet15m, :
## Renaming (Intercept) to Interceptweitrix_confects(cal, design, contrasts=multiple_contrasts, fdr=0.05)## $table
## confect effect timet15m - timet0m timet30m - timet15m timet45m - timet30m
## 1 3.5 5.858 -9.8373 10.8696 1.5141
## 2 3.5 7.524 -0.3369 -0.3586 8.7815
## 3 3.5 5.785 -0.4785 4.0872 -2.2867
## 4 3.5 5.323 -6.8395 10.2452 0.8396
## 5 3.4 7.410 -11.8301 7.4880 4.4171
## 6 3.2 9.821 -8.3607 7.6550 18.4197
## 7 2.8 6.285 -6.8300 10.1441 6.2579
## 8 2.8 4.078 4.8376 -7.2995 -3.3326
## 9 2.7 7.400 2.3901 6.4600 -0.1309
## 10 2.6 3.722 2.6153 1.9706 -0.3507
## timet60m - timet45m timet75m - timet60m timet90m - timet75m
## 1 -0.1804 1.857 -1.2957
## 2 -0.4211 12.594 -5.5116
## 3 7.2778 8.972 -1.4281
## 4 0.1190 5.537 -2.5875
## 5 7.2070 1.792 -5.5915
## 6 -13.1018 18.661 -13.1212
## 7 -2.3809 4.930 -5.4150
## 8 -0.1799 5.257 4.4367
## 9 3.7622 11.352 -15.9122
## 10 2.4888 3.903 0.1883
## timet105m - timet90m timet120m - timet105m fdr_zero row_mean
## 1 3.97229 -0.40015 2.440e-08 29.37
## 2 3.86528 0.20344 6.670e-06 45.91
## 3 -0.07182 -1.28441 1.210e-07 30.67
## 4 1.64288 -1.47056 1.231e-08 38.24
## 5 0.76415 -0.79431 1.661e-05 29.85
## 6 -2.41187 7.39443 2.320e-04 31.90
## 7 2.28033 -2.57263 2.771e-05 36.33
## 8 1.34827 0.02541 1.774e-09 14.49
## 9 14.50129 -5.28822 2.028e-04 39.64
## 10 1.16199 -1.69823 1.774e-09 22.41
## typical_obs_err name gene total_reads
## 1 1.833 YNL036W NCE103 5712
## 2 3.156 YBL097W BRN1 1680
## 3 1.960 YBL016W FUS3 6287
## 4 1.604 YKR093W PTR2 23012
## 5 3.280 YHL021C AIM17 793
## 6 5.191 YDL204W RTN2 418
## 7 2.885 YML100W TSL1 1780
## 8 1.124 YFL014W HSP12 9509
## 9 3.862 YER185W PUG1 813
## 10 1.013 YOR096W RPS7A 72660
## ...
## 511 of 4875 non-zero standard deviation explained at FDR 0.05
## Prior df 21.24.4 Examine individual genes
Having discovered genes with differential tail length, let’s look at some genes in detail.
view_gene <- function(id) {
gene <- rowData(wei)[id,"gene"]
if (is.na(gene)) gene <- ""
tails <- weitrix_x(cal)[id,]
std_errs <- weitrix_weights(cal)[id,] ^ -0.5
ggplot(samples) +
aes(x=time,color=strain,group=strain,
y=tails, ymin=tails-std_errs, ymax=tails+std_errs) +
geom_errorbar(width=0.2) +
geom_hline(yintercept=0) +
geom_line() +
geom_point(aes(size=tail_count[id,])) +
labs(x="Time", y="Tail length", size="Read count", title=paste(id,gene)) +
theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0.5))
}
caption <- plot_annotation(
caption="Error bars show +/- one standard error of measurement.")
# Top confident differences between WT and deltaSet1
view_gene("YDR170W-A") +
view_gene("YIL015W") +
view_gene("YAR009C") +
view_gene("YJR027W/YJR026W") +
caption