Project: r project
.
This report is meant to help explore the results of the derfinder
r Citep(bib[['derfinder']])
package using the single base-level approach and was generated using the regionReport
r Citep(bib[['regionReport']])
package. While the report is rich, it is meant to just start the exploration of the results and exemplify some of the code used to do so. If you need a more in-depth analysis for your specific data set you might want to use the customCode
argument.
Most plots were made with using ggplot2
r Citep(bib[['ggplot2']])
.
## knitrBoostrap and device chunk options library('knitr') opts_chunk$set(bootstrap.show.code = FALSE, dev = device, crop = NULL) if(!outputIsHTML) opts_chunk$set(bootstrap.show.code = FALSE, dev = device, echo = FALSE)
#### Libraries needed ## Bioconductor library('IRanges') library('GenomicRanges') library('GenomeInfoDb') library('derfinder') library('derfinderPlot') library('ggbio') if(hg19) { library('biovizBase') library('TxDb.Hsapiens.UCSC.hg19.knownGene') } ## CRAN library('ggplot2') if(!is.null(theme)) theme_set(theme) library('grid') library('gridExtra') library('knitr') library('RColorBrewer') library('mgcv') library('DT') library('sessioninfo') #### Code setup ## For ggplot tmp <- fullRegions names(tmp) <- seq_len(length(tmp)) regions.df <- as.data.frame(tmp) regions.df$width <- width(tmp) rm(tmp) nulls.df <- as.data.frame(fullNullSummary) ## Special subsets: need at least 3 points for a density plot keepChr <- table(regions.df$seqnames) > 2 regions.df.plot <- subset(regions.df, seqnames %in% names(keepChr[keepChr])) finite.area <- finite.area[regions.df$seqnames %in% names(keepChr[keepChr])] if(hasSig) { ## Keep only those sig regions.df.sig <- regions.df[idx.sig, ] keepChr <- table(regions.df.sig$seqnames) > 2 regions.df.sig <- subset(regions.df.sig, seqnames %in% names(keepChr[keepChr])) if(nrow(regions.df.sig) > 0) { ## If there's any sig, keep those with finite areas if(hasArea) { regions.df.sig.area <- regions.df.sig[is.finite(regions.df.sig$area), ] keepChr <- table(regions.df.sig.area$seqnames) > 2 regions.df.sig.area <- subset(regions.df.sig.area, seqnames %in% names(keepChr[keepChr])) ## Save the info hasArea <- nrow(regions.df.sig.area) > 0 } } else { print('Not a single chromosome had 2 or more significant regions, so plots for significant regions will be skipped.') hasSig <- hasArea <- FALSE } } ## Get chr lengths if(hg19) { seqlengths(fullRegions) <- seqlengths(getChromInfoFromUCSC('hg19', as.Seqinfo = TRUE))[ mapSeqlevels(names(seqlengths(fullRegions)), 'UCSC')] } ## Find which chrs are present in the data set chrs <- levels(seqnames(fullRegions)) ## Subset the fullCoverage data in case that a subset was used colsubset <- optionsStats$colsubset if(!is.null(fullCov) & !is.null(colsubset)) { fullCov <- lapply(fullCov, function(x) { x[, colsubset] }) } ## Get region coverage for the top regions if(nBestRegions > 0) { regionCoverage <- getRegionCoverage(fullCov = fullCov, regions = fullRegions[seq_len(nBestRegions)], chrsStyle = chrsStyle, species = species, currentStyle = currentStyle, verbose = FALSE) save(regionCoverage, file=file.path(workingDir, 'regionCoverage.Rdata')) } ## Graphical setup: transcription database if(hg19 & is.null(txdb)) { txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene::TxDb.Hsapiens.UCSC.hg19.knownGene } else { stopifnot(!is.null(txdb)) }
Theoretically, the p-values should be uniformly distributed between 0 and 1.
p1 <- ggplot(regions.df.plot, aes(x=pvalues, colour=seqnames)) + geom_line(stat='density') + xlim(0, 1) + labs(title='Density of p-values') + xlab('p-values') + scale_colour_discrete(limits=chrs) + theme(legend.title=element_blank()) p1
## Compare the pvalues summary(fullRegions$pvalues)
This is the numerical summary of the distribution of the p-values. r ifelse(nullExist, '', 'Skipped because there are no null regions.')
summary(fullRegions$qvalues)
This is the numerical summary of the distribution of the q-values. r ifelse(nullExist, '', 'Skipped because there are no null regions.')
qtable <- lapply(c(1e-04, 0.001, 0.01, 0.025, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1), function(x) { data.frame('Cut' = x, 'Count' = sum(fullRegions$qvalues <= x)) }) qtable <- do.call(rbind, qtable) if(outputIsHTML) { kable(qtable, format = 'markdown', align = c('c', 'c')) } else { kable(qtable) }
This table shows the number of candidate Differentially Expressed Regions (DERs) with q-value less or equal than some commonly used cutoff values. r ifelse(nullExist, '', 'Skipped because there are no null regions.')
summary(fullRegions$fwer)
This is the numerical summary of the distribution of the FWER adjusted p-values. r ifelse(fwerExist, '', 'Skipped because there are no FWER-adjusted P-values.')
fwertable <- lapply(c(1e-04, 0.001, 0.01, 0.025, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1), function(x) { data.frame('Cut' = x, 'Count' = sum(fullRegions$fwer <= x)) }) fwertable <- do.call(rbind, fwertable) if(outputIsHTML) { kable(fwertable, format = 'markdown', align = c('c', 'c')) } else { kable(fwertable) }
This table shows the number of candidate Differentially Expressed Regions (DERs) with FWER adjusted p-values less or equal than some commonly used cutoff values. r ifelse(fwerExist, '', 'Skipped because there are no FWER-adjusted P-values.')
xrange <- range(log10(regions.df.plot$width)) * c(0.95, 1.05) p2a <- ggplot(regions.df.plot, aes(x=log10(width), colour=seqnames)) + geom_line(stat='density') + labs(title='Density of region lengths') + xlab('Region width (log10)') + scale_colour_discrete(limits=chrs) + xlim(xrange) + theme(legend.title=element_blank()) p2b <- ggplot(regions.df.sig, aes(x=log10(width), colour=seqnames)) + geom_line(stat='density') + labs(title='Density of region lengths (significant only)') + xlab('Region width (log10)') + scale_colour_discrete(limits=chrs) + xlim(xrange) + theme(legend.title=element_blank()) grid.arrange(p2a, p2b)
p2a <- ggplot(regions.df.plot, aes(x=log10(width), colour=seqnames)) + geom_line(stat='density') + labs(title='Density of region lengths') + xlab('Region width (log10)') + scale_colour_discrete(limits=chrs) + theme(legend.title=element_blank()) p2a
This plot shows the density of the region lengths for all regions. r ifelse(hasSig, paste0('The bottom panel is restricted to significant regions (', pvalText, ' < ', sigCut, ')'), '')
xrange <- range(log10(regions.df.plot$area[finite.area])) * c(0.95, 1.05) if(inf.area > 0) { print(paste('Dropping', inf.area, 'due to Inf values.')) } p3a <- ggplot(regions.df.plot[finite.area, ], aes(x=log10(area), colour=seqnames)) + geom_line(stat='density') + labs(title='Density of region areas') + xlab('Region area (log10)') + scale_colour_discrete(limits=chrs) + xlim(xrange) + theme(legend.title=element_blank()) p3b <- ggplot(regions.df.sig.area, aes(x=log10(area), colour=seqnames)) + geom_line(stat='density') + labs(title='Density of region areas (significant only)') + xlab('Region area (log10)') + scale_colour_discrete(limits=chrs) + xlim(xrange) + theme(legend.title=element_blank()) grid.arrange(p3a, p3b)
if(inf.area > 0) { print(paste('Dropping', inf.area, 'due to Inf values.')) } p3a <- ggplot(regions.df.plot[finite.area, ], aes(x=log10(area), colour=seqnames)) + geom_line(stat='density') + labs(title='Density of region areas') + xlab('Region area (log10)') + scale_colour_discrete(limits=chrs) + theme(legend.title=element_blank()) p3a
This plot shows the density of the region areas for all regions. r ifelse(hasArea, paste0('The bottom panel is restricted to significant regions (', pvalText, ' < ', sigCut, ')'), '')
p4 <- ggplot(nulls.df, aes(x=log10(width), colour=chr)) + geom_line(stat='density') + labs(title='Density of null region lengths') + xlab('Region width (log10)') + scale_colour_discrete(limits=chrs) + theme(legend.title=element_blank()) nulls.inf <- !is.finite(nulls.df$area) if(sum(nulls.inf) > 0) { print(paste('Dropping', sum(nulls.inf), 'due to Inf values.')) } p5 <- ggplot(nulls.df[!nulls.inf, ], aes(x=log10(area), colour=chr)) + geom_line(stat='density') + labs(title='Density of null region areas') + xlab('Region area (log10)') + scale_colour_discrete(limits=chrs) + theme(legend.title=element_blank()) grid.arrange(p4, p5)
This plot shows the density of the null region lengths and areas. r ifelse(nullExist, '', 'Skipped because there are no null regions.')
There were a total of r nrow(nulls.df)
null regions.
xrange <- range(log2(regions.df.plot$meanCoverage)) * c(0.95, 1.05) p6a <- ggplot(regions.df.plot, aes(x=log2(meanCoverage), colour=seqnames)) + geom_line(stat='density') + labs(title='Density of region mean coverage') + xlab('Region mean coverage (log2)') + scale_colour_discrete(limits=chrs) + xlim(xrange) + theme(legend.title=element_blank()) p6b <- ggplot(regions.df.sig, aes(x=log2(meanCoverage), colour=seqnames)) + geom_line(stat='density') + labs(title='Density of region mean coverage (significant only)') + xlab('Region mean coverage (log2)') + scale_colour_discrete(limits=chrs) + xlim(xrange) + theme(legend.title=element_blank()) grid.arrange(p6a, p6b)
p6a <- ggplot(regions.df.plot, aes(x=log2(meanCoverage), colour=seqnames)) + geom_line(stat='density') + labs(title='Density of region mean coverage') + xlab('Region mean coverage (log2)') + scale_colour_discrete(limits=chrs) + theme(legend.title=element_blank()) p6a
This plot shows the density of the region mean coverage for all regions. r ifelse(hasSig, paste0('The bottom panel is restricted to significant regions (', pvalText, ' < ', sigCut, ')'), '')
The following plots are MA-style plots comparing each group vs the first one. The mean coverage is calculated using only two groups at a time and is weighted according to the number of samples on each group. Note that the mean coverage and fold change as calculated here do not taking into account the library sizes.
These plots are only shown when there are two or more groups. A total of r length(grep('log2FoldChange', colnames(values(fullRegions))))
plot(s) were made.
for(j in grep('log2FoldChange', colnames(values(fullRegions)))) { ## Identify the groups groups <- strsplit(gsub('log2FoldChange', '', colnames(values(fullRegions))[j]), 'vs')[[1]] ## Calculate the mean coverage only using the 2 groups in question j.mean <- which(colnames(values(fullRegions)) %in% paste0('mean', groups)) groups.n <- sapply(groups, function(x) { sum(optionsStats$groupInfo == x) }) ma.mean.mat <- as.matrix(values(fullRegions)[, j.mean]) ## Weighted means ma.mean <- drop(ma.mean.mat %*% groups.n) / sum(groups.n) + optionsStats$scalefac ma.fold2 <- drop(log2(ma.mean.mat + optionsStats$scalefac) %*% c(1, -1)) ma <- data.frame(mean=ma.mean, log2FoldChange=ma.fold2) ma2 <- ma[is.finite(ma$log2FoldChange), ] fold.mean <- data.frame(foldMean=mean(ma2$log2FoldChange, na.rm=TRUE)) p.ma <- ggplot(ma, aes(x=log2(mean), y=log2FoldChange)) + geom_point(size=1.5, alpha=1/5) + ylab("Fold Change [log2(x + sf)]\nRed dashed line at mean; blue line is GAM fit: y ~ s(x, bs = 'cs')") + xlab(paste('Mean coverage [log2(x + sf)] using only groups', groups[1], 'and', groups[2])) + labs(title=paste('MA style plot:', groups[1], 'vs ', groups[2])) + geom_hline(aes(yintercept=foldMean), data=fold.mean, colour='#990000', linetype='dashed') + geom_smooth(aes(y=log2FoldChange, x=log2(mean)), data=subset(ma2, mean > 0), method = 'gam', formula = y ~ s(x, bs = 'cs')) print(p.ma) }
The following plots were made using ggbio
r Citep(bib[['ggbio']])
which in turn uses ggplot2
r Citep(bib[['ggplot2']])
. For more details check plotOverview
in derfinder
r Citep(bib[['derfinder']])
.
r pvalText
plotOverview(regions=fullRegions, type=pvalVar, base_size=overviewParams$base_size, areaRel=overviewParams$areaRel, legend.position=c(0.97, 0.12))
This plot shows the genomic locations of the candidate regions found in the analysis. The significant regions (r pvalText
less than r sigCut
) are highlighted and the area of the regions is shown on top of each chromosome. Note that the area is in a relative scale.
pvalueVars <- c('P-values' = 'pvalue', 'FDR adj. P-values' = 'qvalue', 'FWER adj. P-values' = 'FWER') regions <- fullRegions pvalueVars <- pvalueVars[pvalueVars %in% colnames(mcols(regions))] for(i in seq_len(length(pvalueVars))) { densityVarName <- names(pvalueVars[i]) densityVarName <- ifelse(is.null(densityVarName), pvalueVars[i], densityVarName) cat(knit_child(text = whisker.render(templateManhattan, list(varName = pvalueVars[i], densityVarName = densityVarName)), quiet = TRUE), sep = '\n') } rm(regions)
plotOverview(regions=fullRegions, annotation=fullRegions, type='annotation', base_size=overviewParams$base_size, areaRel=overviewParams$areaRel, legend.position=c(0.97, 0.12))
This genomic overview plot shows the annotation region type for the candidate regions. Note that the regions are shown only if the annotation information is available. Below is a table of the actual number of results per annotation region type.
annoReg <- table(fullRegions$region, useNA='always') annoReg.df <- data.frame(Region=names(annoReg), Count=as.vector(annoReg)) if(outputIsHTML) { kable(annoReg.df, format = 'markdown', align=rep('c', 3)) } else { kable(annoReg.df) }
r ifelse(hasSig, '## Annotation (significant)', '')
plotOverview(regions=fullRegions[idx.sig], annotation=fullRegions[idx.sig], type='annotation', base_size=overviewParams$base_size, areaRel=overviewParams$areaRel, legend.position=c(0.97, 0.12))
r ifelse(hasSig, paste0('This genomic overview plot shows the annotation region type for the __candidate__ regions that have a ', pvalText, 'less than ', sigCut, '. Note that the regions are shown only if the annotation information is available.'), '')
Below are the plots for the top r nBestRegions
candidate DERs as ranked by area. For each plot, annotation is shown if the candidate DER has a minimum overlap of r optionsMerge$minoverlap
base pairs with annotation information (strand specific). If present, exons are collapsed and shown in blue. Introns are shown in light blue. The title of each plot is composed of the name of the nearest annotation element, the distance to it, and whether the region of the genome the DER falls into; all three pieces of information are based on bumphunter::matchGenes()
.
The annotation depends on the Genomic State used. For details on which one was used for this report check the call to mergeResults
in the reproducibility details.
if(nBestRegions > 0) { plotRegionCoverage(regions = fullRegions, regionCoverage = regionCoverage, groupInfo = optionsStats$groupInfo, nearestAnnotation = regions.df, annotatedRegions = fullAnnotatedRegions, whichRegions = seq_len(min(nBestRegions, length(fullRegions))), colors = NULL, scalefac = optionsStats$scalefac, ask = FALSE, verbose = TRUE, txdb = txdb) }
Below is a table summarizing the number of genomic states per region.
info <- do.call(rbind, lapply(fullAnnotatedRegions$countTable, function(x) { data.frame(table(x)) })) colnames(info) <- c('Number of Overlapping States', 'Frequency') info$State <- gsub('\\..*', '', rownames(info)) rownames(info) <- NULL if(outputIsHTML) { kable(info, format = 'markdown', align=rep('c', 4)) } else { kable(info) }
The following is a venn diagram showing how many regions overlap known exons, introns, and intergenic segments, none of them, or multiple of these groups.
## Venn diagram for all regions venn <- vennRegions(fullAnnotatedRegions, counts.col = 'blue', main = 'Regions overlapping genomic states')
r ifelse(hasSig, 'The following plot is the genomic states venn diagram only for the significant regions.', '')
## Venn diagram for all regions vennSig <- vennRegions(fullAnnotatedRegions, counts.col = 'blue', main = 'Significant regions overlapping genomic states', subsetIndex = idx.sig)
Below is an interactive table with the top r min(nrow(regions.df), nBestRegions * 5)
regions (out of r nrow(regions.df)
) as ranked by area. Inf and -Inf are shown as 1e100 and -1e100 respectively. r ifelse(outputIsHTML, 'Use the search function to find your region of interest or sort by one of the columns.', 'Since the report is in PDF format, only the top 20 regions are shown.')
topArea <- head(regions.df, nBestRegions * 5) topArea <- data.frame('areaRank'=order(topArea$area, decreasing=TRUE), topArea) ## Clean up -Inf, Inf if present ## More details at https://github.com/ramnathv/rCharts/issues/259 replaceInf <- function(df, colsubset=seq_len(ncol(df))) { for(i in colsubset) { inf.idx <- !is.finite(df[, i]) if(any(inf.idx)) { inf.sign <- sign(df[inf.idx, i]) df[inf.idx, i] <- inf.sign * 1e100 } } return(df) } topArea <- replaceInf(topArea, grep('log2FoldChange|value|area', colnames(topArea))) ## Format p-values in scientific notation for(i in which(grepl('pvalues$|qvalues$|fwer$', colnames(topArea)))) topArea[, i] <- format(topArea[, i], scientific = TRUE) ## Make the table if(outputIsHTML) { datatable(topArea, options = list(pagingType='full_numbers', pageLength=10, scrollX='100%'), rownames = FALSE) %>% formatRound(which(grepl('value$|area$|mean|log2FoldChange', colnames(topArea))), digits) } else { ## Only print the top part if your output is a PDF file df_top <- head(topArea, 20) for(i in which(grepl('value$|area$|mean|log2FoldChange', colnames(topArea)))) df_top[, i] <- round(df_top[, i], digits) kable(df_top) }
The following plots were made using ggbio
r Citep(bib[['ggbio']])
which in turn uses ggplot2
r Citep(bib[['ggplot2']])
. For more details check plotCluster()
in derfinder
r Citep(bib[['derfinder']])
.
## Select clusters by cluster area df <- data.frame(area = fullRegions$area, clusterChr = paste0(as.integer(fullRegions$cluster), chr = as.character(seqnames(fullRegions)))) regionClustAreas <- tapply(df$area, df$clusterChr, sum) bestArea <- sapply(names(head(sort(regionClustAreas, decreasing=TRUE), nBestClusters)), function(y) { which(df$clusterChr == y)[[1]]}) ## Graphical setup: ideograms if(hg19 & is.null(p.ideos)) { ## Load ideogram info data(hg19IdeogramCyto, package = 'biovizBase') ideos.set <- as.character(unique(seqnames(fullRegions[bestArea]))) p.ideos <- lapply(ideos.set, function(xx) { plotIdeogram(hg19IdeogramCyto, mapSeqlevels(xx, 'UCSC')) }) names(p.ideos) <- ideos.set } else { stopifnot(!is.null(p.ideos)) } ## Graphical setup: main plotting function regionClusterPlot <- function(idx, tUse='qval') { ## Chr specific selections chr <- as.character(seqnames(fullRegions[idx])) p.ideo <- p.ideos[[chr]] covInfo <- fullCov[[chr]] ## Make the plot p <- plotCluster(idx, regions = fullRegions, annotation = regions.df, coverageInfo = covInfo, groupInfo = optionsStats$groupInfo, titleUse = tUse, txdb = txdb, p.ideogram = p.ideo) print(p) rm(p.ideo, covInfo) return(invisible(TRUE)) }
Below are the best r nBestClusters
region clusters ordered by cluster area (sum of the area of regions inside a cluster). The region with the highest area in the cluster is shown with a red bar. r ifelse(makeBestClusters, '', 'No plots were generated as requested.')
## Genome plots for(idx in bestArea) { regionClusterPlot(idx, ifelse(nullExist, ifelse(pvalVar == 'fwer', ifelse(fwerExist, 'fwer', 'qval'), pvalVar), 'none')) }
Below is the information on how the samples were permuted.
## Get the permutation information nSamples <- seq_len(length(optionsStats$groupInfo)) permuteInfo <- lapply(seeds, function(x) { set.seed(x) idx <- sample(nSamples) data.frame(optionsStats$groupInfo[idx]) }) permuteInfo <- cbind(data.frame(optionsStats$groupInfo), do.call(cbind, permuteInfo)) colnames(permuteInfo) <- c('original', paste0('perm', seq_len(optionsStats$nPermute))) ## The raw information # permuteInfo n <- names(table(permuteInfo[, 2])) permuteDetail <- data.frame(matrix(NA, nrow=optionsStats$nPermute * length(n), ncol = 2 + length(n))) permuteDetail[, 1] <- rep(seq_len(optionsStats$nPermute), each=length(n)) permuteDetail[, 2] <- rep(n, optionsStats$nPermute) colnames(permuteDetail) <- c('permutation', 'group', as.character(n)) l <- 1 m <- 3:ncol(permuteDetail) for(j in n) { k <- which(permuteInfo[, 1] == j) for(i in 2:(optionsStats$nPermute + 1)) { permuteDetail[l, m] <- table(permuteInfo[k, i]) l <- l + 1 } } ## How many permutations resulted in the original grouping rearrangement obs <- diag(length(m)) * sapply( permuteDetail$group[ permuteDetail$permutation == 1], function(n) { sum(optionsStats$groupInfo == n) }) sameAsObs <- sapply(seq_len(length(seeds)), function(i) { p <- as.matrix(permuteDetail[permuteDetail$permutation == i, m]) all((p - obs) == 0) }) ## Print the summary summary(permuteDetail[, m])
This table shows the summary per group of how many samples were assigned to the group. It can be used for fast detection of anomalies. r ifelse(usedPermutations, paste('Also note that', sum(sameAsObs), 'permutations out of', length(seeds), 'total permutations resulted in the same grouping as in the original observed data.'), 'Skipped because no permutations were used.')
Note that in derfinder
the re-sampling of the samples is done without replacement. This is done to avoid singular model matrices. While the sample balance is the same across the permutations, what changes are the adjusted variables (including the column medians).
The following table shows how the group labels were permuted. r ifelse(usedPermutations, '', 'Skipped because no permutations were used.')
This can be useful to detect whether a permutation in particular had too many samples of a group labeled as another group, meaning that the resulting permuted group label resulted in pretty much a name change.
if(outputIsHTML) { datatable(permuteDetail, options = list(pagingType='full_numbers', pageLength=10, scrollX='100%')) } else { ## Only print the top part if your output is a PDF file kable(head(permuteDetail, 20)) }
The F-statistic cutoff used was r as.character(optionsStats$cutoffFstat)
and type of cutoff used was r optionsStats$cutoffType
. Furthermore, the maximum region (data) gap was set to r optionsStats$maxRegionGap
and the maximum cluster gap was set to r optionsStats$maxClusterGap
.
This analysis was on each chromosome was performed with the following call to analyzeChr()
(shown for one chromosome only):
if('analyzeCall' %in% names(optionsStats)) { optionsStats$analyzeCall } else { 'Skipped since this information was not recorded prior to version 0.0.24' }
The results were merged using the following call to mergeResults()
:
optionsMerge$mergeCall
This report was generated in path r tmpdir
using the following call to derfinderReport()
:
theCall
Date the report was generated.
## Date the report was generated Sys.time()
Wallclock time spent generating the report.
## Processing time in seconds totalTime <- diff(c(startTime, Sys.time())) round(totalTime, digits=3)
R
session information.
## Session info options(width = 120) session_info()
Pandoc version used: r rmarkdown::pandoc_version()
.
This report was created with regionReport
r Citep(bib[['regionReport']])
using rmarkdown
r Citep(bib[['rmarkdown']])
while knitr
r Citep(bib[['knitr']])
and DT
r Citep(bib[['DT']])
were running behind the scenes.
Citations made with r CRANpkg('RefManageR')
r Citep(bib[['RefManageR']])
. The BibTeX file can be found here.
## Print bibliography PrintBibliography(bib, .opts = list(hyperlink = "to.doc", style = "html"))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.