<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta charset="utf-8" />
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="pandoc" />
<meta name="author" content="Joyce Hsiao" />
<title>Evaluate projected cell time</title>
<script src="site_libs/jquery-1.11.3/jquery.min.js"></script>
<meta name="viewport" content="width=device-width, initial-scale=1" />
<link href="site_libs/bootstrap-3.3.5/css/cosmo.min.css" rel="stylesheet" />
<script src="site_libs/bootstrap-3.3.5/js/bootstrap.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/html5shiv.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/respond.min.js"></script>
<script src="site_libs/navigation-1.1/tabsets.js"></script>
<link href="site_libs/highlightjs-9.12.0/textmate.css" rel="stylesheet" />
<script src="site_libs/highlightjs-9.12.0/highlight.js"></script>
<link href="site_libs/font-awesome-4.5.0/css/font-awesome.min.css" rel="stylesheet" />
<style type="text/css">code{white-space: pre;}</style>
<style type="text/css">
pre:not([class]) {
background-color: white;
}
</style>
<script type="text/javascript">
if (window.hljs) {
hljs.configure({languages: []});
hljs.initHighlightingOnLoad();
if (document.readyState && document.readyState === "complete") {
window.setTimeout(function() { hljs.initHighlighting(); }, 0);
}
}
</script>
<style type="text/css">
h1 {
font-size: 34px;
}
h1.title {
font-size: 38px;
}
h2 {
font-size: 30px;
}
h3 {
font-size: 24px;
}
h4 {
font-size: 18px;
}
h5 {
font-size: 16px;
}
h6 {
font-size: 12px;
}
.table th:not([align]) {
text-align: left;
}
</style>
</head>
<body>
<style type = "text/css">
.main-container {
max-width: 940px;
margin-left: auto;
margin-right: auto;
}
code {
color: inherit;
background-color: rgba(0, 0, 0, 0.04);
}
img {
max-width:100%;
height: auto;
}
.tabbed-pane {
padding-top: 12px;
}
button.code-folding-btn:focus {
outline: none;
}
</style>
<style type="text/css">
/* padding for bootstrap navbar */
body {
padding-top: 51px;
padding-bottom: 40px;
}
/* offset scroll position for anchor links (for fixed navbar) */
.section h1 {
padding-top: 56px;
margin-top: -56px;
}
.section h2 {
padding-top: 56px;
margin-top: -56px;
}
.section h3 {
padding-top: 56px;
margin-top: -56px;
}
.section h4 {
padding-top: 56px;
margin-top: -56px;
}
.section h5 {
padding-top: 56px;
margin-top: -56px;
}
.section h6 {
padding-top: 56px;
margin-top: -56px;
}
</style>
<script>
// manage active state of menu based on current page
$(document).ready(function () {
// active menu anchor
href = window.location.pathname
href = href.substr(href.lastIndexOf('/') + 1)
if (href === "")
href = "index.html";
var menuAnchor = $('a[href="' + href + '"]');
// mark it active
menuAnchor.parent().addClass('active');
// if it's got a parent navbar menu mark it active as well
menuAnchor.closest('li.dropdown').addClass('active');
});
</script>
<div class="container-fluid main-container">
<!-- tabsets -->
<script>
$(document).ready(function () {
window.buildTabsets("TOC");
});
</script>
<!-- code folding -->
<div class="navbar navbar-default navbar-fixed-top" role="navigation">
<div class="container">
<div class="navbar-header">
<button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar">
<span class="icon-bar"></span>
<span class="icon-bar"></span>
<span class="icon-bar"></span>
</button>
<a class="navbar-brand" href="index.html">fucci-seq</a>
</div>
<div id="navbar" class="navbar-collapse collapse">
<ul class="nav navbar-nav">
<li>
<a href="index.html">Home</a>
</li>
<li>
<a href="about.html">About</a>
</li>
<li>
<a href="license.html">License</a>
</li>
</ul>
<ul class="nav navbar-nav navbar-right">
<li>
<a href="https://github.com/jdblischak/fucci-seq">
<span class="fa fa-github"></span>
</a>
</li>
</ul>
</div><!--/.nav-collapse -->
</div><!--/.container -->
</div><!--/.navbar -->
<!-- Add a small amount of space between sections. -->
<style type="text/css">
div.section {
padding-top: 12px;
}
</style>
<div class="fluid-row" id="header">
<h1 class="title toc-ignore">Evaluate projected cell time</h1>
<h4 class="author"><em>Joyce Hsiao</em></h4>
</div>
<div id="TOC">
<ul>
<li><a href="#estimate-cell-time">Estimate cell time</a></li>
<li><a href="#results">Results</a></li>
<li><a href="#save-re-ordered-cell-times">Save re-ordered cell times</a></li>
<li><a href="#session-information">Session information</a></li>
</ul>
</div>
<!-- The file analysis/chunks.R contains chunks that define default settings
shared across the workflowr files. -->
<!-- Update knitr chunk options -->
<!-- Insert the date the file was last updated -->
<p><strong>Last updated:</strong> 2018-04-11</p>
<!-- Insert the code version (Git commit SHA1) if Git repository exists and R
package git2r is installed -->
<p><strong>Code version:</strong> cde0b19</p>
<hr />
<div id="estimate-cell-time" class="section level2">
<h2>Estimate cell time</h2>
<pre class="r"><code>library(Biobase)
# load gene expression
df <- readRDS(file="../data/eset-final.rds")
pdata <- pData(df)
fdata <- fData(df)
log2cpm.all <- t(log2(1+(10^6)*(t(exprs(df))/pdata$molecules)))
macosko <- readRDS("../data/cellcycle-genes-previous-studies/rds/macosko-2015.rds")
pc.fucci <- prcomp(subset(pdata,
select=c("rfp.median.log10sum.adjust",
"gfp.median.log10sum.adjust")),
center = T, scale. = T)
Theta.cart <- pc.fucci$x
library(circular)
Theta.fucci <- coord2rad(Theta.cart)
Theta.fucci <- 2*pi - Theta.fucci</code></pre>
<p>Cluster cell times to move the origin of the cell times</p>
<pre class="r"><code># cluster cell time
library(movMF)
clust.res <- lapply(2:5, function(k) {
movMF(Theta.cart, k=k, nruns = 100, kappa = list(common = TRUE))
})
k.list <- sapply(clust.res, function(x) length(x$theta) + length(x$alpha) + 1)
bic <- sapply(1:length(clust.res), function(i) {
x <- clust.res[[i]]
k <- k.list[i]
n <- nrow(Theta.cart)
-2*x$L + k *(log(n) - log(2*pi)) })
plot(bic)
labs <- predict(clust.res[[2]])
saveRDS(labs, file = "../output/images-time-eval.Rmd/labs.rds")</code></pre>
<pre class="r"><code>labs <- readRDS(file = "../output/images-time-eval.Rmd/labs.rds")
summary(as.numeric(Theta.fucci)[labs==1])</code></pre>
<pre><code> Min. 1st Qu. Median Mean 3rd Qu. Max.
2.885 3.537 3.941 3.896 4.243 5.069 </code></pre>
<pre class="r"><code>summary(as.numeric(Theta.fucci)[labs==2])</code></pre>
<pre><code> Min. 1st Qu. Median Mean 3rd Qu. Max.
1.179 1.646 2.091 2.063 2.455 2.861 </code></pre>
<pre class="r"><code>summary(as.numeric(Theta.fucci)[labs==3])</code></pre>
<pre><code> Min. 1st Qu. Median Mean 3rd Qu. Max.
0.004025 0.304041 1.112189 3.000234 5.923906 6.282360 </code></pre>
<pre class="r"><code># move the origin to 1.17
Theta.fucci.new <- vector("numeric", length(Theta.fucci))
cutoff <- min(Theta.fucci[labs==2])
Theta.fucci.new[Theta.fucci>=cutoff] <- Theta.fucci[Theta.fucci>=cutoff] - cutoff
Theta.fucci.new[Theta.fucci<cutoff] <- Theta.fucci[Theta.fucci<cutoff] - cutoff + 2*pi</code></pre>
<p>Try plotting for one gene</p>
<pre class="r"><code>macosko[macosko$hgnc == "CDK1",]</code></pre>
<pre><code> hgnc phase ensembl
113 CDK1 G2 ENSG00000170312</code></pre>
<pre class="r"><code>cdk1 <- log2cpm.all[rownames(log2cpm.all)=="ENSG00000170312",]
plot(x=Theta.fucci.new, y = cdk1)
points(y=cdk1[labs==1], x=as.numeric(Theta.fucci.new)[labs==1], pch=16, cex=.7, col = "red")
points(y=cdk1[labs==2], x=as.numeric(Theta.fucci.new)[labs==2], pch=16, cex=.7, col = "forestgreen")
points(y=cdk1[labs==3], x=as.numeric(Theta.fucci.new)[labs==3], pch=16, cex=.7, col = "blue")</code></pre>
<p><img src="figure/images-time-eval.Rmd/unnamed-chunk-4-1.png" width="672" style="display: block; margin: auto;" /></p>
<hr />
</div>
<div id="results" class="section level2">
<h2>Results</h2>
<p>Check data points with outlier DAPI values</p>
<pre class="r"><code>ii.min.dapi <- order(pData(df)$dapi.median.log10sum.adjust)[1:2]
pData(df)[ii.min.dapi,]</code></pre>
<pre><code> experiment well cell_number concentration ERCC
20170924-E08 20170924 E08 1 1.125631 50x dilution
20170924-F03 20170924 F03 1 2.157178 50x dilution
individual.1 individual.2 image_individual image_label
20170924-E08 NA18870 NA19160 19160_18870 29
20170924-F03 NA18870 NA19160 19160_18870 33
raw umi mapped unmapped reads_ercc reads_hs
20170924-E08 4494424 3083676 2197773 885903 169730 2027183
20170924-F03 4884157 3348082 2349374 998708 158751 2189654
reads_egfp reads_mcherry molecules mol_ercc mol_hs mol_egfp
20170924-E08 783 77 134694 3220 131433 39
20170924-F03 957 12 156481 3373 153055 46
mol_mcherry detect_ercc detect_hs chip_id chipmix freemix
20170924-E08 2 40 8284 NA19160 0.21325 0.08510
20170924-F03 7 40 8691 NA18870 0.39507 0.15026
snps reads avg_dp min_dp snps_w_min valid_id cut_off_reads
20170924-E08 311848 7531 0.02 1 3354 TRUE TRUE
20170924-F03 311848 8223 0.03 1 3694 TRUE TRUE
unmapped_ratios cut_off_unmapped ercc_percentage cut_off_ercc
20170924-E08 0.2872880 TRUE 0.07722818 TRUE
20170924-F03 0.2982926 TRUE 0.06757162 TRUE
cut_off_genes ercc_conversion conversion conversion_outlier
20170924-E08 TRUE 0.01897131 0.06483529 FALSE
20170924-F03 TRUE 0.02124711 0.06989917 FALSE
molecule_outlier filter_all rfp.median.log10sum
20170924-E08 FALSE TRUE 1.702458
20170924-F03 FALSE TRUE 2.155286
gfp.median.log10sum dapi.median.log10sum
20170924-E08 1.421271 1.176911
20170924-F03 1.421271 1.146540
rfp.median.log10sum.adjust gfp.median.log10sum.adjust
20170924-E08 -0.54086625 -1.333926
20170924-F03 0.07550985 -1.212614
dapi.median.log10sum.adjust size perimeter eccentricity
20170924-E08 -1.561447 703 107 0.9337102
20170924-F03 -1.476240 405 62 0.6811795
theta
20170924-E08 2.070977
20170924-F03 2.382573</code></pre>
<pre class="r"><code>par(mfrow=c(2,2))
ylims <- with(pdata, range(c(dapi.median.log10sum.adjust,
gfp.median.log10sum.adjust,
rfp.median.log10sum.adjust)))
plot(as.numeric(Theta.fucci.new), pdata$dapi.median.log10sum.adjust, col = "blue",
ylab= "adjusted intensity values",
ylim = ylims, main = "DAPI intensity values",
xlab ="Estimated cell time")
points(as.numeric(Theta.fucci.new)[ii.min.dapi],
pdata$dapi.median.log10sum.adjust[ii.min.dapi], pch=4, lwd=2, cex=1)
plot(as.numeric(Theta.fucci.new), pdata$gfp.median.log10sum.adjust, col = "forestgreen",
ylab= "adjusted intensity values",
ylim = ylims, main = "GFP and RFP intensity values",
xlab ="Estimated cell time")
points(as.numeric(Theta.fucci.new)[ii.min.dapi],
pdata$gfp.median.log10sum.adjust[ii.min.dapi], pch=4, lwd=2, cex=1)
points(as.numeric(Theta.fucci.new), pdata$rfp.median.log10sum.adjust, col = "red")
points(as.numeric(Theta.fucci.new)[ii.min.dapi],
pdata$rfp.median.log10sum.adjust[ii.min.dapi], pch=4, lwd=2, cex=1)
plot(as.numeric(Theta.fucci.new), pdata$molecules, main = "Total molecule count",
xlab ="Estimated cell time", ylab = "Total molecule count", col = "gray50")
points(as.numeric(Theta.fucci.new)[ii.min.dapi],
pdata$molecules[ii.min.dapi], pch=4, lwd=2, cex=1)
plot(pdata$dapi.median.log10sum.adjust, pdata$molecules, main = "DAPI vs. molecule count",
xlab = "DAPI intensity adjusted for C1 batch", ylab = "Total molecule count",
col = "gray50")
points(pdata$dapi.median.log10sum.adjust[ii.min.dapi],
pdata$molecules[ii.min.dapi], pch=4, lwd=2, cex=1)</code></pre>
<p><img src="figure/images-time-eval.Rmd/unnamed-chunk-6-1.png" width="672" style="display: block; margin: auto;" /></p>
<p>Test the association between total sample molecule count and DAPI.</p>
<ul>
<li><p>After excluding outliers, pearson correlation is .2</p></li>
<li><p>Consider lm(molecules ~ dapi). The adjusted R-squared is .04</p></li>
<li><p>Consider lm(log10(molecules) ~ dapi). The adjusted R-squared is .04</p></li>
<li><p>Weak linear trend between molecule count and DAPI…</p></li>
</ul>
<pre class="r"><code>xy <- data.frame(dapi=pdata$dapi.median.log10sum.adjust,
molecules=pdata$molecules,
chip_id=pdata$chip_id)
xy <- xy[xy$dapi > -1,]
fit <- lm(molecules~dapi+factor(chip_id), data=xy)
summary(fit)</code></pre>
<pre><code>
Call:
lm(formula = molecules ~ dapi + factor(chip_id), data = xy)
Residuals:
Min 1Q Median 3Q Max
-54492 -15575 -1546 13723 71917
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 128330.1 2191.8 58.550 < 2e-16 ***
dapi 17081.8 3732.2 4.577 5.40e-06 ***
factor(chip_id)NA18855 -7233.1 2838.5 -2.548 0.01100 *
factor(chip_id)NA18870 2640.0 2680.6 0.985 0.32495
factor(chip_id)NA19098 -833.6 2733.1 -0.305 0.76043
factor(chip_id)NA19101 17348.1 2976.0 5.829 7.81e-09 ***
factor(chip_id)NA19160 9631.3 3034.5 3.174 0.00156 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 22160 on 879 degrees of freedom
Multiple R-squared: 0.1335, Adjusted R-squared: 0.1276
F-statistic: 22.57 on 6 and 879 DF, p-value: < 2.2e-16</code></pre>
<pre class="r"><code>fit <- lm(log10(molecules)~dapi+factor(chip_id), data=xy)
summary(fit)</code></pre>
<pre><code>
Call:
lm(formula = log10(molecules) ~ dapi + factor(chip_id), data = xy)
Residuals:
Min 1Q Median 3Q Max
-0.20393 -0.04962 0.00043 0.04868 0.21492
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.103490 0.007289 700.187 < 2e-16 ***
dapi 0.058028 0.012411 4.675 3.39e-06 ***
factor(chip_id)NA18855 -0.027545 0.009439 -2.918 0.00361 **
factor(chip_id)NA18870 0.006707 0.008914 0.752 0.45197
factor(chip_id)NA19098 -0.005623 0.009089 -0.619 0.53626
factor(chip_id)NA19101 0.054951 0.009896 5.553 3.73e-08 ***
factor(chip_id)NA19160 0.030816 0.010091 3.054 0.00233 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0737 on 879 degrees of freedom
Multiple R-squared: 0.1371, Adjusted R-squared: 0.1312
F-statistic: 23.28 on 6 and 879 DF, p-value: < 2.2e-16</code></pre>
<pre class="r"><code>fit <- lm(molecules~dapi, data=xy)
summary(fit)</code></pre>
<pre><code>
Call:
lm(formula = molecules ~ dapi, data = xy)
Residuals:
Min 1Q Median 3Q Max
-49959 -16957 -585 15862 86073
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 131005.0 784.9 166.902 < 2e-16 ***
dapi 20835.5 3794.1 5.492 5.21e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 23350 on 884 degrees of freedom
Multiple R-squared: 0.03299, Adjusted R-squared: 0.0319
F-statistic: 30.16 on 1 and 884 DF, p-value: 5.209e-08</code></pre>
<pre class="r"><code>fit <- lm(log10(molecules)~dapi, data=xy)
summary(fit)</code></pre>
<pre><code>
Call:
lm(formula = log10(molecules) ~ dapi, data = xy)
Residuals:
Min 1Q Median 3Q Max
-0.201411 -0.053816 0.004652 0.056963 0.222330
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.110141 0.002613 1955.683 < 2e-16 ***
dapi 0.071475 0.012630 5.659 2.06e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.07772 on 884 degrees of freedom
Multiple R-squared: 0.03496, Adjusted R-squared: 0.03387
F-statistic: 32.02 on 1 and 884 DF, p-value: 2.057e-08</code></pre>
<pre class="r"><code>cor(xy$dapi, xy$molecules, method = "pearson")</code></pre>
<pre><code>[1] 0.1816292</code></pre>
<pre class="r"><code>cor(xy$dapi, log10(xy$molecules), method = "pearson")</code></pre>
<pre><code>[1] 0.1869765</code></pre>
<pre class="r"><code>par(mfrow=c(1,2))
plot(x=xy$dapi, y = xy$molecules,
xlab = "DAPI intensity", ylab = "Total molecule count",
main = "DAPI vs. molecule count \n adj-Rsq = .03; cor=.2")
abline(lm(molecules~dapi, data=xy), col = "red")
plot(x=xy$dapi, y = log10(xy$molecules),
xlab = "DAPI intensity", ylab = "log10 total molecule count",
main = "DAPI vs. log10 molecule count \n adj-Rsq = .03, cor=.2")
abline(lm(log10(molecules)~dapi, data=xy), col = "red")</code></pre>
<p><img src="figure/images-time-eval.Rmd/unnamed-chunk-7-1.png" width="672" style="display: block; margin: auto;" /></p>
<hr />
</div>
<div id="save-re-ordered-cell-times" class="section level2">
<h2>Save re-ordered cell times</h2>
<pre class="r"><code>theta <- as.numeric(Theta.fucci.new)
names(theta) <- colnames(log2cpm.all)
saveRDS(theta,
file = "../output/images-time-eval.Rmd/theta.rds")</code></pre>
<hr />
</div>
<div id="session-information" class="section level2">
<h2>Session information</h2>
<pre class="r"><code>sessionInfo()</code></pre>
<pre><code>R version 3.4.1 (2017-06-30)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: Scientific Linux 7.2 (Nitrogen)
Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] circular_0.4-93 Biobase_2.38.0 BiocGenerics_0.24.0
loaded via a namespace (and not attached):
[1] Rcpp_0.12.16 mvtnorm_1.0-7 digest_0.6.15 rprojroot_1.3-2
[5] backports_1.1.2 git2r_0.21.0 magrittr_1.5 evaluate_0.10.1
[9] stringi_1.1.7 boot_1.3-19 rmarkdown_1.9 tools_3.4.1
[13] stringr_1.3.0 yaml_2.1.18 compiler_3.4.1 htmltools_0.3.6
[17] knitr_1.20 </code></pre>
</div>
<!-- Adjust MathJax settings so that all math formulae are shown using
TeX fonts only; see
http://docs.mathjax.org/en/latest/configuration.html. This will make
the presentation more consistent at the cost of the webpage sometimes
taking slightly longer to load. Note that this only works because the
footer is added to webpages before the MathJax javascript. -->
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
"HTML-CSS": { availableFonts: ["TeX"] }
});
</script>
<hr>
<p>
This <a href="http://rmarkdown.rstudio.com">R Markdown</a> site was created with <a href="https://github.com/jdblischak/workflowr">workflowr</a>
</p>
<hr>
<!-- To enable disqus, uncomment the section below and provide your disqus_shortname -->
<!-- disqus
<div id="disqus_thread"></div>
<script type="text/javascript">
/* * * CONFIGURATION VARIABLES: EDIT BEFORE PASTING INTO YOUR WEBPAGE * * */
var disqus_shortname = 'rmarkdown'; // required: replace example with your forum shortname
/* * * DON'T EDIT BELOW THIS LINE * * */
(function() {
var dsq = document.createElement('script'); dsq.type = 'text/javascript'; dsq.async = true;
dsq.src = '//' + disqus_shortname + '.disqus.com/embed.js';
(document.getElementsByTagName('head')[0] || document.getElementsByTagName('body')[0]).appendChild(dsq);
})();
</script>
<noscript>Please enable JavaScript to view the <a href="http://disqus.com/?ref_noscript">comments powered by Disqus.</a></noscript>
<a href="http://disqus.com" class="dsq-brlink">comments powered by <span class="logo-disqus">Disqus</span></a>
-->
</div>
<script>
// add bootstrap table styles to pandoc tables
function bootstrapStylePandocTables() {
$('tr.header').parent('thead').parent('table').addClass('table table-condensed');
}
$(document).ready(function () {
bootstrapStylePandocTables();
});
</script>
<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
(function () {
var script = document.createElement("script");
script.type = "text/javascript";
script.src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
document.getElementsByTagName("head")[0].appendChild(script);
})();
</script>
</body>
</html>