MSD plots
Razvan Sultana
1/31/2018
Load the packages used throughout this analysis.
setwd("/share/users/Mike/AutismPaper/RNASeq/R")
suppressMessages({
library('readr')
library('stringr')
library('tibble')
library('dplyr')
library('tidyr')
library('tximport')
library('BiocParallel')
library('DESeq2')
library('ggplot2')
library('RColorBrewer')
library('colorspace')
library('ggrepel')
library('gridExtra')
library('vsn')
library('MASS')
library('plotly')
library('htmlwidgets')
register(MulticoreParam(detectCores()))
options(tibble.print_max=100)
})MSD plots
for Bottomly dataset
2D
metadata.bot <- read.csv('/share/users/Razvan/sleuth_walkthroughs/bottomly/metadata/experiment.csv', stringsAsFactors = FALSE)
extract_metadata <- function(library_name) {
ret <- lapply(strsplit(library_name, '_'),
function(x) {
data.frame(strain = x[1], experiment = x[2], lane = x[3],
stringsAsFactors = FALSE)
})
dplyr::bind_rows(ret)
}
metadata.bot <- dplyr::select(metadata.bot, -strain)
metadata.bot <- dplyr::bind_cols(metadata.bot, extract_metadata(metadata.bot$library_name))
metadata.bot <- dplyr::select(metadata.bot, run_accession, library_name, strain,
experiment, lane)
metadata.bot <- dplyr::mutate(metadata.bot,
path = file.path('/share/users/Razvan/sleuth_walkthroughs/bottomly', 'results', run_accession, 'kallisto','abundance.h5'))
metadata.bot <- dplyr::rename(metadata.bot, sample = run_accession)
files.bot <- metadata.bot %>% dplyr::select('sample', 'path') %>% spread('sample', 'path') %>% as.list() %>% unlist()
txi.kallisto.bot <- tximport(files.bot, type = "kallisto", txOut = TRUE)
metadata.bot.df <- metadata.bot %>% column_to_rownames(var="sample") %>% as.data.frame()
rownames(txi.kallisto.bot$counts) <- as.character(rownames(txi.kallisto.bot$counts))
dds.bot <- DESeqDataSetFromTximport(txi = txi.kallisto.bot, colData = metadata.bot.df, design = ~strain)
dds.bot <- DESeq(dds.bot, test = "Wald", betaPrior = TRUE, parallel = TRUE)
c <- counts(dds.bot) %>% as.tibble()
c[['target_id']] <- rownames(counts(dds.bot))
straind2 <- metadata.bot %>% filter(grepl('^D2', strain)) %>% dplyr::select(sample) %>% as.list() %>% unlist() %>% as.vector()
strainb6 <- metadata.bot %>% filter(grepl('^B6', strain)) %>% dplyr::select(sample) %>% as.list() %>% unlist() %>% as.vector()
c <- c %>% mutate(sum1 = rowSums(dplyr::select(., one_of(straind2))),
sum2 = rowSums(dplyr::select(., one_of(strainb6))))
c <- c %>% filter((sum1 > 1) & (sum2 > 1))
dds.bot <- dds.bot[c$target_id,]
rm(list=c('c', 'straind2', 'strainb6'))
resultsNames(dds.bot)## [1] "Intercept" "strainB6" "strainD2"rld.bot <- rlog(dds.bot, blind = FALSE)
vst.bot <- vst(dds.bot, blind = FALSE)msd0.bot <- meanSdPlot(log2(txi.kallisto.bot$counts+1), plot = FALSE)
msd1.bot <- meanSdPlot(log2(counts(dds.bot, normalized=FALSE)+1), plot = FALSE)
msd2.bot <- meanSdPlot(log2(counts(dds.bot, normalized=TRUE)+1), plot=FALSE)
msd3.bot <- meanSdPlot(assay(rld.bot), plot=FALSE)
msd4.bot <- meanSdPlot(assay(vst.bot), plot=FALSE)
pl.bot <- list(
msd1.bot$gg + ggtitle("Unnormalized raw counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd2.bot$gg + ggtitle("Normalized raw counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd3.bot$gg + ggtitle("Rlog transformed counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd4.bot$gg + ggtitle("Vst transformed counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8))
)
ml.bot <- marrangeGrob(pl.bot, nrow = 2, ncol = 2)
ml.bot
3D
d1.bot <- with(msd1.bot$gg$data, MASS::kde2d(px, py, n = 100))
d2.bot <- with(msd2.bot$gg$data, MASS::kde2d(px, py, n = 100))
d3.bot <- with(msd3.bot$gg$data, MASS::kde2d(px, py, n = 100))
d4.bot <- with(msd4.bot$gg$data, MASS::kde2d(px, py, n = 100))
label1 <- list(
text = "Raw, unnormalized counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.25, y = 0.5, showarrow = FALSE
)
label2 <- list(
text = "Normalized counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.75, y = 0.5, showarrow = FALSE
)
label3 <- list(
text = "rlog transformed counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.25, y = 0, showarrow = FALSE
)
label4 <- list(
text = "vst transformed counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.75, y = 0, showarrow = FALSE
)
p1.bot <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene1') %>%
add_surface(x=d1.bot$y, y=d1.bot$x, z=d1.bot$z, cauto=FALSE, cmax=max(d1.bot$z))
p2.bot <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene2') %>%
add_surface(x=d2.bot$y, y=d2.bot$x, z=d2.bot$z, cauto=FALSE, cmax=max(d1.bot$z))
p3.bot <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene3') %>%
add_surface(x=d3.bot$y, y=d3.bot$x, z=d3.bot$z, cauto=FALSE, cmax=max(d1.bot$z))
p4.bot <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene4') %>%
add_surface(x=d4.bot$y, y=d4.bot$x, z=d4.bot$z, cauto=FALSE, cmax=max(d1.bot$z))
axx <- list(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=TRUE,
backgroundcolor='rgb(230, 230,230)',
range=c(0,4)
)
axy <- list(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=TRUE,
backgroundcolor='rgb(230, 230,230)'
)
axz <- list(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=TRUE,
backgroundcolor='rgb(230, 230,230)',
range=c(0,5e-5)
)
ply.bot <- subplot(p1.bot, p2.bot, p3.bot, p4.bot) %>%
layout(title = "Bottomly dataset",
scene = list(domain=list(x=c(0,0.5), y=c(0.525,1)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
scene2 = list(domain=list(x=c(0.5,1), y=c(0.525,1)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
scene3 = list(domain=list(x=c(0,0.5), y=c(0.025,0.475)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
scene4 = list(domain=list(x=c(0.5,1), y=c(0.025,0.475)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
annotations = list(label1, label2, label3, label4))
ply.bot %>%
onRender(
"function(el, x) {
var gd = document.getElementById(el.id);
var cnt = 0;
var scene = gd._fullLayout['scene']._scene;
var scene2 = gd._fullLayout['scene2']._scene;
var scene3 = gd._fullLayout['scene3']._scene;
var scene4 = gd._fullLayout['scene4']._scene;
gd.on('plotly_relayout', function() {
var camera = scene.getCamera();
scene2.setCamera(camera);
scene3.setCamera(camera);
scene4.setCamera(camera);
})
function run() {
rotate('scene', Math.PI / 180);
rotate('scene2', Math.PI / 180);
rotate('scene3', Math.PI / 180);
rotate('scene4', Math.PI / 180);
requestAnimationFrame(run);
}
//run();
function rotate(id, angle) {
var scene = gd._fullLayout[id]._scene;
var camera = scene.getCamera();
var rtz = xyz2rtz(camera.eye);
rtz.t += angle;
camera.eye = rtz2xyz(rtz);
scene.setCamera(camera);
}
function xyz2rtz(xyz) {
return {
r: Math.sqrt(xyz.x * xyz.x + xyz.y * xyz.y),
t: Math.atan(xyz.y / xyz.x),
z: xyz.z
};
}
function rtz2xyz(rtz) {
return {
x: rtz.r * Math.cos(rtz.t),
y: rtz.r * Math.sin(rtz.t),
z: rtz.z
};
}
}"
)for our dataset
meta_int <- read_tsv('../../metadata/RNASeq.internal.tsv') %>%
mutate(library = 'paired end') %>%
mutate(readlength = 101) %>%
mutate(tissue = 'SVZ') %>%
mutate(path = paste0('/share/users/Mike/AutismPaper/RNASeq/OurData/',
sample, '.refseqnorrna.kallisto.paired/abundance.h5')
)
meta_ext <- read_tsv('../../metadata/RNASeq.external.tsv') %>%
mutate(library = 'single end') %>%
mutate(readlength = 75) %>%
mutate(condition = 'Control') %>%
mutate(path = paste0('/share/users/Mike/AutismPaper/RNASeq/external/',
sample, '.refseqnorrna.kallisto.single/abundance.h5')
)
meta_ext1 <-read_tsv('../../metadata/RNASeq.external2.tsv') %>%
mutate(path = paste0('/share/users/Mike/AutismPaper/RNASeq/external2/',
sample, '.refseqnorrna.kallisto.single/abundance.h5')
)
metadata <- bind_rows(meta_int, meta_ext, meta_ext1) %>%
mutate(age = factor(if_else(!is.na(age_wpc), 'Fetal',
if_else(between(age_years,1,19),'Child','Adult')),
levels = c('Fetal', 'Child', 'Adult')
)
) %>%
mutate(tissue1 = if_else(grepl('SVZ$', tissue), 'SVZ',
if_else(grepl('frontal', tissue), 'FC',
if_else(grepl('temporal', tissue), 'TC',
tissue)
)
)
) %>%
mutate(sample = if_else(grepl('OurData', path),
paste(sep="_",sub("(\\w).*", "\\1", condition),
paste0(age_years, 'y'),
sample),
if_else(grepl('external2',path),
paste(sep='_', sub("(\\w).*", "\\1", condition),
tissue1,
paste0(age_years, 'y'),
sub(".*(\\d{2})$", "\\1", sample)
),
paste(sep='_', tissue,
paste0(age_wpc, 'w'),
sub(".*(\\d{2})$", "\\1", sample)
)
)
)
) %>%
filter(is.na(sex) | sex == 'M')
metadata1 <- metadata %>%
filter(grepl('SVZ',tissue) & !grepl('H16$', sample)) %>%
mutate(condition = factor(condition, levels=c('Control', 'Autism'))) %>%
mutate(grp = factor(paste0(age,condition), levels=c('FetalControl', 'ChildControl', 'ChildAutism', 'AdultControl', 'AdultAutism'))) %>%
arrange(sample) %>%
dplyr::select( c('sample', 'condition', 'library', 'readlength', 'tissue', 'age', 'grp', 'path'))
files <- metadata1 %>% dplyr::select('sample', 'path') %>% spread('sample', 'path') %>% as.list() %>% unlist()
txi.kallisto <- tximport(files, type = "kallisto", txOut = TRUE)
metadata1.df <- metadata1 %>% column_to_rownames(var="sample") %>% as.data.frame()
rownames(txi.kallisto$counts) <- as.character(rownames(txi.kallisto$counts))
ttgfile <- '/share/db/refseq/human/refseq.108.rna.t2g.txt'
ttg <- read_tsv(ttgfile)
dds <- DESeqDataSetFromTximport(txi = txi.kallisto, colData = metadata1.df, design = ~grp)
c <- counts(dds) %>% as.tibble()
c[['target_id']] <- rownames(counts(dds))
c <- c %>% mutate(sum1 = rowSums(dplyr::select(., contains("SVZ"))),
sum2 = rowSums(dplyr::select(., starts_with("A_"))),
sum3 = rowSums(dplyr::select(., starts_with("C_"))))
c <- c %>% filter((sum1 > 1) & (sum2 > 1) & (sum3 > 1))
dds <- dds[c$target_id,]
rm('c')
ddsWald <- DESeq(dds, test = "Wald", betaPrior = TRUE, parallel = TRUE)
resultsChild <- results(ddsWald, alpha = 0.05,
contrast = c('grp', 'ChildControl', 'FetalControl'))
resultsAdult <- results(ddsWald, alpha = 0.05,
contrast = c('grp', 'AdultControl', 'FetalControl'))
resultsChildA <- results(ddsWald, alpha = 0.05,
contrast = c('grp', 'ChildAutism', 'FetalControl'))
resultsAdultA <- results(ddsWald, alpha = 0.05,
contrast = c('grp', 'AdultAutism', 'FetalControl'))
resultsChildCA <- results(ddsWald, alpha = 0.05,
contrast = c('grp','ChildControl','ChildAutism'))
resultsAdultCA <- results(ddsWald, alpha = 0.05,
contrast = c('grp','AdultControl','AdultAutism'))
signifChild <- resultsChild %>%
as.tibble(.) %>% rownames_to_column(var='target_id') %>%
filter(padj<0.05) %>% arrange(padj) %>%
left_join(ttg, by = 'target_id') %>%
dplyr::select(c('target_id', 'baseMean', 'log2FoldChange', 'lfcSE', 'stat', 'pvalue', 'padj', 'gene_id', 'symbol', 'description', 'status', 'molecule_type'))
signifAdult <- resultsAdult %>%
as.tibble(.) %>% rownames_to_column(var='target_id') %>%
filter(padj<0.05) %>% arrange(padj) %>%
left_join(ttg, by = 'target_id') %>%
dplyr::select(c('target_id', 'baseMean', 'log2FoldChange', 'lfcSE', 'stat', 'pvalue', 'padj', 'gene_id', 'symbol', 'description', 'status', 'molecule_type'))
signifChildA <- resultsChildA %>%
as.tibble(.) %>% rownames_to_column(var='target_id') %>%
filter(padj<0.05) %>% arrange(padj) %>%
left_join(ttg, by = 'target_id') %>%
dplyr::select(c('target_id', 'baseMean', 'log2FoldChange', 'lfcSE', 'stat', 'pvalue', 'padj', 'gene_id', 'symbol', 'description', 'status', 'molecule_type'))
signifAdultA <- resultsAdultA %>%
as.tibble(.) %>% rownames_to_column(var='target_id') %>%
filter(padj<0.05) %>% arrange(padj) %>%
left_join(ttg, by = 'target_id') %>%
dplyr::select(c('target_id', 'baseMean', 'log2FoldChange', 'lfcSE', 'stat', 'pvalue', 'padj', 'gene_id', 'symbol', 'description', 'status', 'molecule_type'))
signifChildCA <- resultsChildCA %>%
as.tibble(.) %>% rownames_to_column(var='target_id') %>%
filter(padj<0.05) %>% arrange(padj) %>%
left_join(ttg, by = 'target_id') %>%
dplyr::select(c('target_id', 'baseMean', 'log2FoldChange', 'lfcSE', 'stat', 'pvalue', 'padj', 'gene_id', 'symbol', 'description', 'status', 'molecule_type'))
signifAdultCA <- resultsAdultCA %>%
as.tibble(.) %>% rownames_to_column(var='target_id') %>%
filter(padj<0.05) %>% arrange(padj) %>%
left_join(ttg, by = 'target_id') %>%
dplyr::select(c('target_id', 'baseMean', 'log2FoldChange', 'lfcSE', 'stat', 'pvalue', 'padj', 'gene_id', 'symbol', 'description', 'status', 'molecule_type'))
devChildCA <- intersect(signifChild$target_id, signifChildCA$target_id)
devAdultCA <- intersect(signifAdult$target_id, signifAdultCA$target_id)
devChildCAg <- tibble(target_id=devChildCA) %>% left_join(ttg, by='target_id') %>%
dplyr::select('gene_id') %>% as.list() %>% unlist() %>% unique()
devAdultCAg <- tibble(target_id=devAdultCA) %>% left_join(ttg, by='target_id') %>%
dplyr::select('gene_id') %>% as.list() %>% unlist() %>% unique()
devChildCAs <- tibble(target_id=devChildCA) %>% left_join(ttg, by='target_id') %>%
dplyr::select('symbol') %>% as.list() %>% unlist() %>% unique()
devAdultCAs <- tibble(target_id=devAdultCA) %>% left_join(ttg, by='target_id') %>%
dplyr::select('symbol') %>% as.list() %>% unlist() %>% unique()
intgroup = "grp"
rld <- rlog(ddsWald, blind = FALSE)
vst <- varianceStabilizingTransformation(ddsWald, blind=FALSE)```
2D
msd0 <- meanSdPlot(log2(txi.kallisto$counts+1), plot = FALSE)
msd1 <- meanSdPlot(log2(counts(ddsWald, normalized=FALSE)+1), plot = FALSE)
msd2 <- meanSdPlot(log2(counts(ddsWald, normalized=TRUE)+1), plot=FALSE)
msd3 <- meanSdPlot(assay(rld), plot=FALSE)
msd4 <- meanSdPlot(assay(vst), plot=FALSE)
pl <- list(
msd1$gg + ggtitle("Unnormalized raw counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd2$gg + ggtitle("Normalized raw counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd3$gg + ggtitle("Rlog transformed counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd4$gg + ggtitle("Vst transformed counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8))
)
ml <- marrangeGrob(pl, nrow = 2, ncol = 2)
ml
3D
d1 <- with(msd1$gg$data, MASS::kde2d(px, py, n = 100))
d2 <- with(msd2$gg$data, MASS::kde2d(px, py, n = 100))
d3 <- with(msd3$gg$data, MASS::kde2d(px, py, n = 100))
d4 <- with(msd4$gg$data, MASS::kde2d(px, py, n = 100))
label1 <- list(
text = "Raw, unnormalized counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.25, y = 0.5, showarrow = FALSE
)
label2 <- list(
text = "Normalized counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.75, y = 0.5, showarrow = FALSE
)
label3 <- list(
text = "rlog transformed counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.25, y = 0, showarrow = FALSE
)
label4 <- list(
text = "vst transformed counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.75, y = 0, showarrow = FALSE
)
p1 <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene1') %>%
add_surface(x=d1$y, y=d1$x, z=d1$z, cauto=FALSE, cmax=max(d1$z))
p2 <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene2') %>%
add_surface(x=d2$y, y=d2$x, z=d2$z, cauto=FALSE, cmax=max(d1$z))
p3 <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene3') %>%
add_surface(x=d3$y, y=d3$x, z=d3$z, cauto=FALSE, cmax=max(d1$z))
p4 <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene4') %>%
add_surface(x=d4$y, y=d4$x, z=d4$z, cauto=FALSE, cmax=max(d1$z))
axx <- list(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=TRUE,
backgroundcolor='rgb(230, 230,230)',
range=c(0,4)
)
axy <- list(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=TRUE,
backgroundcolor='rgb(230, 230,230)'
)
axz <- list(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=TRUE,
backgroundcolor='rgb(230, 230,230)',
range=c(0,5e-5)
)
ply <- subplot(p1, p2, p3, p4) %>%
layout(title = "Autism SVZ dataset",
scene = list(domain=list(x=c(0,0.5), y=c(0.525,1)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
scene2 = list(domain=list(x=c(0.5,1), y=c(0.525,1)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
scene3 = list(domain=list(x=c(0,0.5), y=c(0.025,0.475)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
scene4 = list(domain=list(x=c(0.5,1), y=c(0.025,0.475)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
annotations = list(label1, label2, label3, label4))
ply %>%
onRender(
"function(el, x) {
var gd = document.getElementById(el.id);
var cnt = 0;
var scene = gd._fullLayout['scene']._scene;
var scene2 = gd._fullLayout['scene2']._scene;
var scene3 = gd._fullLayout['scene3']._scene;
var scene4 = gd._fullLayout['scene4']._scene;
gd.on('plotly_relayout', function() {
var camera = scene.getCamera();
scene2.setCamera(camera);
scene3.setCamera(camera);
scene4.setCamera(camera);
})
function run() {
rotate('scene', Math.PI / 180);
rotate('scene2', Math.PI / 180);
rotate('scene3', Math.PI / 180);
rotate('scene4', Math.PI / 180);
requestAnimationFrame(run);
}
//run();
function rotate(id, angle) {
var scene = gd._fullLayout[id]._scene;
var camera = scene.getCamera();
var rtz = xyz2rtz(camera.eye);
rtz.t += angle;
camera.eye = rtz2xyz(rtz);
scene.setCamera(camera);
}
function xyz2rtz(xyz) {
return {
r: Math.sqrt(xyz.x * xyz.x + xyz.y * xyz.y),
t: Math.atan(xyz.y / xyz.x),
z: xyz.z
};
}
function rtz2xyz(rtz) {
return {
x: rtz.r * Math.cos(rtz.t),
y: rtz.r * Math.sin(rtz.t),
z: rtz.z
};
}
}"
)Fetal samples
metadata.fetal <- metadata1 %>% filter(age == 'Fetal')
files.fetal <- metadata.fetal %>% dplyr::select('sample', 'path') %>% spread('sample', 'path') %>% as.list() %>% unlist()
txi.kallisto.fetal <- tximport(files.fetal, type = "kallisto", txOut = TRUE)
metadata.fetal.df <- metadata.fetal %>% column_to_rownames(var="sample") %>% as.data.frame()
rownames(txi.kallisto.fetal$counts) <- as.character(rownames(txi.kallisto.fetal$counts))
dds.fetal <- DESeqDataSetFromTximport(txi = txi.kallisto.fetal, colData = metadata.fetal.df, design = ~tissue)
dds.fetal <- DESeq(dds.fetal, test = "Wald", betaPrior = TRUE, parallel = TRUE)
c <- counts(dds.fetal) %>% as.tibble()
c[['target_id']] <- rownames(counts(dds.fetal))
c <- c %>% mutate(sum = rowSums(dplyr::select(., contains('SVZ'))))
c <- c %>% filter((sum > 1))
dds.fetal <- dds.fetal[c$target_id,]
rm('c')
resultsNames(dds.fetal)## [1] "Intercept" "tissueISVZ" "tissueOSVZ"rld.fetal <- tryCatch(rlog(dds.fetal, blind = FALSE), error = function(e) { rlog(dds.fetal, fitType = 'mean', blind = FALSE) })
vst.fetal <- tryCatch(vst(dds.fetal, blind = FALSE), error = function(e) { vst(dds.fetal, fitType = 'mean', blind = FALSE) })2D
msd0.fetal <- meanSdPlot(log2(txi.kallisto.fetal$counts+1), plot = FALSE)
msd1.fetal <- meanSdPlot(log2(counts(dds.fetal, normalized=FALSE)+1), plot = FALSE)
msd2.fetal <- meanSdPlot(log2(counts(dds.fetal, normalized=TRUE)+1), plot=FALSE)
msd3.fetal <- meanSdPlot(assay(rld.fetal), plot=FALSE)
msd4.fetal <- meanSdPlot(assay(vst.fetal), plot=FALSE)
pl.fetal <- list(
msd0.fetal$gg + ggtitle("Raw kallisto estimated counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd1.fetal$gg + ggtitle("Unnormalized raw counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd2.fetal$gg + ggtitle("Normalized raw counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd3.fetal$gg + ggtitle("Rlog transformed counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd4.fetal$gg + ggtitle("Vst transformed counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8))
)
ml.fetal <- marrangeGrob(pl.fetal, nrow = 2, ncol = 3)
ml.fetal
3D
d1.fetal <- with(msd1.fetal$gg$data, MASS::kde2d(px, py, n = 100))
d2.fetal <- with(msd2.fetal$gg$data, MASS::kde2d(px, py, n = 100))
d3.fetal <- with(msd3.fetal$gg$data, MASS::kde2d(px, py, n = 100))
d4.fetal <- with(msd4.fetal$gg$data, MASS::kde2d(px, py, n = 100))
label1 <- list(
text = "Raw, unnormalized counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.25, y = 0.5, showarrow = FALSE
)
label2 <- list(
text = "Normalized counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.75, y = 0.5, showarrow = FALSE
)
label3 <- list(
text = "rlog transformed counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.25, y = 0, showarrow = FALSE
)
label4 <- list(
text = "vst transformed counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.75, y = 0, showarrow = FALSE
)
p1.fetal <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene1') %>%
add_surface(x=d1.fetal$y, y=d1.fetal$x, z=d1.fetal$z, cauto=FALSE, cmax=max(d1.fetal$z))
p2.fetal <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene2') %>%
add_surface(x=d2.fetal$y, y=d2.fetal$x, z=d2.fetal$z, cauto=FALSE, cmax=max(d1.fetal$z))
p3.fetal <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene3') %>%
add_surface(x=d3.fetal$y, y=d3.fetal$x, z=d3.fetal$z, cauto=FALSE, cmax=max(d1.fetal$z))
p4.fetal <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene4') %>%
add_surface(x=d4.fetal$y, y=d4.fetal$x, z=d4.fetal$z, cauto=FALSE, cmax=max(d1.fetal$z))
axx <- list(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=TRUE,
backgroundcolor='rgb(230, 230,230)',
range=c(0,4)
)
axy <- list(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=TRUE,
backgroundcolor='rgb(230, 230,230)'
)
axz <- list(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=TRUE,
backgroundcolor='rgb(230, 230,230)',
range=c(0,5e-5)
)
ply.fetal <- subplot(p1.fetal, p2.fetal, p3.fetal, p4.fetal) %>%
layout(title = "Autism SVZ dataset - fetal",
scene = list(domain=list(x=c(0,0.5), y=c(0.525,1)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
scene2 = list(domain=list(x=c(0.5,1), y=c(0.525,1)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
scene3 = list(domain=list(x=c(0,0.5), y=c(0.025,0.475)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
scene4 = list(domain=list(x=c(0.5,1), y=c(0.025,0.475)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
annotations = list(label1, label2, label3, label4))
ply.fetal %>%
onRender(
"function(el, x) {
var gd = document.getElementById(el.id);
var cnt = 0;
var scene = gd._fullLayout['scene']._scene;
var scene2 = gd._fullLayout['scene2']._scene;
var scene3 = gd._fullLayout['scene3']._scene;
var scene4 = gd._fullLayout['scene4']._scene;
gd.on('plotly_relayout', function() {
var camera = scene.getCamera();
scene2.setCamera(camera);
scene3.setCamera(camera);
scene4.setCamera(camera);
})
function run() {
rotate('scene', Math.PI / 180);
rotate('scene2', Math.PI / 180);
rotate('scene3', Math.PI / 180);
rotate('scene4', Math.PI / 180);
requestAnimationFrame(run);
}
//run();
function rotate(id, angle) {
var scene = gd._fullLayout[id]._scene;
var camera = scene.getCamera();
var rtz = xyz2rtz(camera.eye);
rtz.t += angle;
camera.eye = rtz2xyz(rtz);
scene.setCamera(camera);
}
function xyz2rtz(xyz) {
return {
r: Math.sqrt(xyz.x * xyz.x + xyz.y * xyz.y),
t: Math.atan(xyz.y / xyz.x),
z: xyz.z
};
}
function rtz2xyz(rtz) {
return {
x: rtz.r * Math.cos(rtz.t),
y: rtz.r * Math.sin(rtz.t),
z: rtz.z
};
}
}"
)Adult samples
metadata.adult <- metadata1 %>% filter(age != 'Fetal')
files.adult <- metadata.adult %>% dplyr::select('sample', 'path') %>% spread('sample', 'path') %>% as.list() %>% unlist()
txi.kallisto.adult <- tximport(files.adult, type = "kallisto", txOut = TRUE)
metadata.adult.df <- metadata.adult %>% column_to_rownames(var="sample") %>% as.data.frame()
rownames(txi.kallisto.adult$counts) <- as.character(rownames(txi.kallisto.adult$counts))
dds.adult <- DESeqDataSetFromTximport(txi = txi.kallisto.adult, colData = metadata.adult.df, design = ~condition)
dds.adult <- DESeq(dds.adult, test = "Wald", betaPrior = TRUE, parallel = TRUE)
c <- counts(dds.adult) %>% as.tibble()
c[['target_id']] <- rownames(counts(dds.adult))
c <- c %>% mutate(sum1 = rowSums(dplyr::select(., starts_with('C_'))),
sum2 = rowSums(dplyr::select(., starts_with('A_'))))
c <- c %>% filter((sum1 > 1) & (sum2 > 1))
dds.adult <- dds.adult[c$target_id,]
rm('c')
resultsNames(dds.adult)## [1] "Intercept" "conditionControl" "conditionAutism"rld.adult <- tryCatch(rlog(dds.adult, blind = FALSE), error = function(e) { rlog(dds.adult, fitType = 'mean', blind = FALSE) })
vst.adult <- tryCatch(vst(dds.adult, blind = FALSE), error = function(e) { vst(dds.adult, fitType = 'mean', blind = FALSE) })2D
msd0.adult <- meanSdPlot(log2(txi.kallisto.adult$counts+1), plot = FALSE)
msd1.adult <- meanSdPlot(log2(counts(dds.adult, normalized=FALSE)+1), plot = FALSE)
msd2.adult <- meanSdPlot(log2(counts(dds.adult, normalized=TRUE)+1), plot=FALSE)
msd3.adult <- meanSdPlot(assay(rld.adult), plot=FALSE)
msd4.adult <- meanSdPlot(assay(vst.adult), plot=FALSE)
pl.adult <- list(
msd0.adult$gg + ggtitle("Raw kallisto estimated counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd1.adult$gg + ggtitle("Unnormalized raw counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd2.adult$gg + ggtitle("Normalized raw counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd3.adult$gg + ggtitle("Rlog transformed counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8)),
msd4.adult$gg + ggtitle("Vst transformed counts") + scale_fill_gradient(low = "yellow", high = "darkred") + scale_y_continuous(limits = c(0, 8))
)
ml.adult <- marrangeGrob(pl.adult, nrow = 2, ncol = 3)
ml.adult
3D
d0.adult <- with(msd0.adult$gg$data, MASS::kde2d(px, py, n = 100))
d1.adult <- with(msd1.adult$gg$data, MASS::kde2d(px, py, n = 100))
d2.adult <- with(msd2.adult$gg$data, MASS::kde2d(px, py, n = 100))
d3.adult <- with(msd3.adult$gg$data, MASS::kde2d(px, py, n = 100))
d4.adult <- with(msd4.adult$gg$data, MASS::kde2d(px, py, n = 100))
label1 <- list(
text = "Raw, unnormalized counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.25, y = 0.5, showarrow = FALSE
)
label2 <- list(
text = "Normalized counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.75, y = 0.5, showarrow = FALSE
)
label3 <- list(
text = "rlog transformed counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.25, y = 0, showarrow = FALSE
)
label4 <- list(
text = "vst transformed counts", font = list(size = 12),
xref = "paper", yref = "paper",
xanchor = "center", yanchor = "center", align = "center",
x = 0.75, y = 0, showarrow = FALSE
)
p1.adult <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene1') %>%
add_surface(x=d1.adult$y, y=d1.adult$x, z=d1.adult$z, cauto=FALSE, cmax=max(d1.adult$z))
p2.adult <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene2') %>%
add_surface(x=d2.adult$y, y=d2.adult$x, z=d2.adult$z, cauto=FALSE, cmax=max(d1.adult$z))
p3.adult <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene3') %>%
add_surface(x=d3.adult$y, y=d3.adult$x, z=d3.adult$z, cauto=FALSE, cmax=max(d1.adult$z))
p4.adult <- plot_ly(showscale=FALSE, reversescale=TRUE, scene='scene4') %>%
add_surface(x=d4.adult$y, y=d4.adult$x, z=d4.adult$z, cauto=FALSE, cmax=max(d1.adult$z))
axx <- list(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=TRUE,
backgroundcolor='rgb(230, 230,230)',
range=c(0,4)
)
axy <- list(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=TRUE,
backgroundcolor='rgb(230, 230,230)'
)
axz <- list(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=TRUE,
backgroundcolor='rgb(230, 230,230)',
range=c(0,5e-5)
)
ply.adult <- subplot(p1.adult, p2.adult, p3.adult, p4.adult) %>%
layout(title = "Autism SVZ dataset - adult",
scene = list(domain=list(x=c(0,0.5), y=c(0.525,1)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
scene2 = list(domain=list(x=c(0.5,1), y=c(0.525,1)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
scene3 = list(domain=list(x=c(0,0.5), y=c(0.025,0.475)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
scene4 = list(domain=list(x=c(0.5,1), y=c(0.025,0.475)),
xaxis=axx, yaxis=axy, zaxis=axz,
aspectmode='cube'),
annotations = list(label1, label2, label3, label4))
ply.adult %>%
onRender(
"function(el, x) {
var gd = document.getElementById(el.id);
var cnt = 0;
var scene = gd._fullLayout['scene']._scene;
var scene2 = gd._fullLayout['scene2']._scene;
var scene3 = gd._fullLayout['scene3']._scene;
var scene4 = gd._fullLayout['scene4']._scene;
gd.on('plotly_relayout', function() {
var camera = scene.getCamera();
scene2.setCamera(camera);
scene3.setCamera(camera);
scene4.setCamera(camera);
})
function run() {
rotate('scene', Math.PI / 180);
rotate('scene2', Math.PI / 180);
rotate('scene3', Math.PI / 180);
rotate('scene4', Math.PI / 180);
requestAnimationFrame(run);
}
//run();
function rotate(id, angle) {
var scene = gd._fullLayout[id]._scene;
var camera = scene.getCamera();
var rtz = xyz2rtz(camera.eye);
rtz.t += angle;
camera.eye = rtz2xyz(rtz);
scene.setCamera(camera);
}
function xyz2rtz(xyz) {
return {
r: Math.sqrt(xyz.x * xyz.x + xyz.y * xyz.y),
t: Math.atan(xyz.y / xyz.x),
z: xyz.z
};
}
function rtz2xyz(rtz) {
return {
x: rtz.r * Math.cos(rtz.t),
y: rtz.r * Math.sin(rtz.t),
z: rtz.z
};
}
}"
)