Create weitrix object
From experience, noise scales with the length of the tail. Therefore to stabilize the variance we will be examing log2 of the tail length.
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
## FALSE TRUE
## 1657 4875
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.
The dispersion is then calculated for each row. A spline curve is fitted to the log dispersion. Weights in each row are then scaled so as to flatten this curve.
Examining the effect of calibration
Let’s unpack what happens in that weitrix_calibrate_trend step.
Consider first the dispersion if weights were uniform. Missing data is weighted 0, all other weights are 1.
There is information from the number of reads that we can remove.
Here are the dispersions using read-counts as weights.
In general it is improved, but now we have the opposite problem. There is a component of noise that does not go away with more and more reads. The red line is a trend fitted to this, which weitrix_calibrate_trend divides out of the weights.
Finally, here are the dispersion from the calibrated weitrix.
This is reasonably close to uniform, with no trend from the number of reads. limma can now estimate the variability of dispersions between genes, and apply its Emprical-Bayes squeezed dispersion based testing.
Testing
We are now ready to test things. We feed our calibrated weitrix to limma.
## gene diff_log2_tail ave_log2_tail adj.P.Val total_reads
## YDR170W-A <NA> -0.3317487 4.679897 3.440790e-06 4832
## YJR027W/YJR026W <NA> -0.2564543 4.271790 2.287980e-05 6577
## YAR009C <NA> -0.4673765 3.931000 2.287980e-05 1866
## YIL015W BAR1 -0.3653600 4.893617 6.233632e-05 3898
## YPL257W-B/YPL257W-A <NA> -0.2558543 4.209088 1.151177e-04 5729
## YML045W/YML045W-A <NA> -0.2406881 4.354082 1.953731e-04 4553
## YLR157C-B/YLR157C-A <NA> -0.2224134 4.200638 3.860716e-04 5519
## YGR027W-B/YGR027W-A <NA> -0.2093877 4.206404 4.750243e-04 5685
## YLR035C-A <NA> -0.3804618 3.918629 5.945463e-04 1374
## YPR137C-B/YPR137C-A <NA> -0.2105818 4.257820 6.441298e-04 7287
## YER138C/YER137C-A <NA> -0.2053578 4.187405 7.537499e-04 5652
## YML039W/YML040W <NA> -0.2343804 4.213771 7.928725e-04 5851
## YMR045C/YMR046C <NA> -0.2087868 4.266254 8.387704e-04 4579
## YJR029W/YJR028W <NA> -0.2062738 4.254686 1.128685e-03 5979
## YNL284C-B/YNL284C-A <NA> -0.1899551 4.370570 1.157811e-03 4830
## YPL217C BMS1 -0.3051718 4.086526 1.254310e-03 2071
## YPR158C-D/YPR158C-C <NA> -0.1975358 4.227524 1.254310e-03 5381
## YIL053W GPP1 -0.1463690 4.760486 1.480738e-03 25996
## YOR192C-B/YOR192C-A <NA> -0.4356656 4.397721 1.480738e-03 864
## YER133W GLC7 -0.7326781 4.084170 7.537499e-04 264
My package topconfects can be used to find top confident differential tail length. Rather than picking “most significant” genes, it will highlight genes with a large effect size.
## gene diff_log2_tail confect total_reads
## YAR009C <NA> -0.4673765 -0.157 1866
## YDR170W-A <NA> -0.3317487 -0.151 4832
## YER133W GLC7 -0.7326781 -0.151 264
## YIL015W BAR1 -0.3653600 -0.131 3898
## YJR027W/YJR026W <NA> -0.2564543 -0.107 6577
## YLR035C-A <NA> -0.3804618 -0.102 1374
## YCR014C POL4 -0.4968279 -0.098 606
## YPL257W-B/YPL257W-A <NA> -0.2558543 -0.094 5729
## YOR192C-B/YOR192C-A <NA> -0.4356656 -0.094 864
## YML045W/YML045W-A <NA> -0.2406881 -0.085 4553
## YPL217C BMS1 -0.3051718 -0.074 2071
## YLR157C-B/YLR157C-A <NA> -0.2224134 -0.074 5519
## YGR027W-B/YGR027W-A <NA> -0.2093877 -0.068 5685
## YML039W/YML040W <NA> -0.2343804 -0.068 5851
## YDR476C <NA> -0.3794150 -0.067 4196
## YPR137C-B/YPR137C-A <NA> -0.2105818 -0.067 7287
## YER138C/YER137C-A <NA> -0.2053578 -0.063 5652
## YMR045C/YMR046C <NA> -0.2087868 -0.063 4579
## YJR029W/YJR028W <NA> -0.2062738 -0.060 5979
## YMR238W DFG5 -0.2478286 -0.060 2502
## 80 genes significantly non-zero at FDR 0.05
This lists the largest confident log2 fold changes in poly(A) tail length. The confect column is an inner confidence bound on the difference in log2 tail length, adjusted for multiple testing.
We discover some genes with less total reads, but large change in tail length.
Note that due to PCR amplification slippage and limited read length, the observed log2 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.
Examine individual genes
Having discovered genes with differential tail length, let’s look at some genes in detail.
view_gene <- function(id, title="") {
ggplot(samples, aes(x=time,color=strain,group=strain, y=tail[id,])) +
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,title))
}
# Top "significant" genes
view_gene("YDR170W-A")
view_gene("YJR027W/YJR026W")
view_gene("YAR009C")
view_gene("YIL015W","BAR1")
# topconfects has highlighted some genes with lower total reads
view_gene("YER133W","GLC7")
view_gene("YCR014C","POL4")
Exploratory analysis
The test we’ve performed was somewhat unsatisfactory. Due to the design of the experiment it’s difficul to specify differential tests that fully interrogate this dataset: the lack of replicates, and the difficult specifying apriori how tail length will change over time.
Perhaps we should let the data speak for itself.
Perhaps this is what we should have done first!
The weitrix package allows us to look for components of variation. We’ll try to explain the data with different numbers of components (from 1 to 10 components).
weitrix_seq_screeplot shows how much additional variation in the data is explained as each further component is allowed. However the ultimate decision of how many components to examine is a matter of judgement.
Looking at three components shows some of the major trends in this data-set.
comp <- comp_seq[[3]]
matrix_long(comp$col[,-1], row_info=samples, varnames=c("sample","component")) %>%
ggplot(aes(x=time, y=value, color=strain, group=strain)) +
geom_hline(yintercept=0) +
geom_line() +
geom_point(alpha=0.75, size=3) +
facet_grid(component ~ .) +
labs(title="Sample scores for each component", y="Sample score", x="Time", color="Strain")
We observe:
- C1 - A gradual lengthening of tails after release into cell cycling. (The reason for the divergence between strains at the end is unclear.)
- C2 - A lengthening of poly(A) tails in the set1 mutant.
- C3 - Variation in poly(A) tail length with the cell cycle.
The log2 tail lengths are approximated by comp$row %*% t(comp$col) where comp$col is an \(n_\text{sample} \times (p+1)\) matrix of scores (shown above), and comp$row is an \(n_\text{gene} \times (p+1)\) matrix of gene loadings, which we will now examine. (The \(+1\) is the intercept “component”, allowing each gene to have a different baseline tail length.)
Treat these results with caution. Confindence bounds take into account uncertainty in the loadings but not in the scores! What follows is best regarded as exploratory rather than a final result.
