发那个太丢人

导航

 

1)简介:

DESeq2-package: for differential analysis of count data(对count data 做差异分析)

2)安装

if("DESeq2" %in% rownames(installed.packages()) == FALSE) {source("http://bioconductor.org/biocLite.R");biocLite("DESeq2")}
suppressMessages(library(DESeq2))
ls('package:DESeq2')

 3)对象的使用说明

3.1)coef(Extract a matrix of model coefficients/standard errors,高级用户检验模型系数)

语法:coef(object, SE = FALSE, ...)

参数解释:

object:a DESeqDataSet returned by DESeq, nbinomWaldTest, or nbinomLRT.

例子:

dds <- makeExampleDESeqDataSet(m=4)
dds <- DESeq(dds)
coef(dds)[1,]
coef(dds, SE=TRUE)[1,]

 

3.2) collapseReplicates:Collapse technical replicates in a RangedSummarizedExperiment or DESeqDataSet(用于消除技术重复)
用法:collapseReplicates(object, groupby, run, renameCols = TRUE)
参数:

object:A RangedSummarizedExperiment or DESeqDataSet
groupby:a grouping factor, as long as the columns of object,分组因子
run:optional, the names of each unique column in object. if provided, a new column runsCollapsed will be added to the colData which pastes together the names of run (测序run)
renameCols:whether to rename the columns of the returned object using the levels of the grouping factor

例子:

dds <- makeExampleDESeqDataSet(m=12)
str(dds)
dds$sample <- factor(sample(paste0("sample",rep(1:9, c(2,1,1,2,1,1,2,1,1))))) (#共9个样品:其中 3个样品有2个技术重重)
dds$run <- paste0("run",1:12) #12个run道
ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
# examine the colData and column names of the collapsed data
colData(ddsColl)
colnames(ddsColl)
# check that the sum of the counts for "sample1" is the same
# as the counts in the "sample1" column in ddsColl
matchFirstLevel <- dds$sample == levels(dds$sample)[1]
stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))

 

3.3)counts:Accessors for the ’counts’ slot of a DESeqDataSet object(对表达矩阵进行统计,)

one row for each observational unit (gene or the like), and one column for each sample(行代表观察值(例如基因),列代表样本(例如肝、脾、肾等))

 语法:counts(object, normalized = FALSE,replaced = FALSE)

 参数:

object:a DESeqDataSet object(表达矩阵).
normalized:logical indicating whether or not to divide the counts by the size factors or normalization factors before returning (normalization factors always preempt size factors),(即不同量级的数据要不要归一化)
replaced:返回极端值

dds <- makeExampleDESeqDataSet(m=4)  ##构建一个表达矩阵
head(counts(dds))
dds <- estimateSizeFactors(dds) # run this or DESeq() first  
head(counts(dds, normalized=TRUE))

 

3.4)DESeq:Differential expression analysis based on the Negative Binomial (a.k.a.Gamma-Poisson) distribution(基于负二项分布进行差异分析)

语法:

DESeq(object, test = c("Wald", "LRT"), fitType = c("parametric", "local","mean"), sfType = c("ratio", "poscounts", "iterate"), betaPrior,full = design(object), reduced, quiet = FALSE,minReplicatesForReplace = 7, modelMatrixType, useT = FALSE, minmu = 0.5,
parallel = FALSE, BPPARAM = bpparam())

参数:

object:a DESeqDataSet object(表达矩阵对象)
test:Wald" or "LRT"检验
fitType:either "parametric", "local", or "mean"
sfType:either "ratio", "poscounts", or "iterate" for teh type of size factor estimation.
betaPrior:whether or not to put a zero-mean normal prior on the non-intercept coefficients
reduced:for test="LRT", a reduced formula to compare against
quiet:whether to print messages at each step
minReplicatesForReplace:the minimum number of replicates required
modelMatrixType:either "standard" or "expanded", which describe how the model matrix, X of the GLM formula is formed.
useT:logical, passed to nbinomWaldTest, default is FALSE
minmu:lower bound on the estimated count for fitting gene-wise dispersion
parallel:if FALSE, no parallelization. if TRUE, parallel execution using BiocParallel,
BPPARAM:an optional parameter object passed internally to bplapply when parallel=TRUE.
例子:

# count tables from RNA-Seq data
cnts <- matrix(rnbinom(n=1000, mu=100, size=1/0.5), ncol=10)
cond <- factor(rep(1:2, each=5))

# object construction
dds <- DESeqDataSetFromMatrix(cnts, DataFrame(cond), ~ cond)

# standard analysis
dds <- DESeq(dds)
res <- results(dds)

# moderated log2 fold changes
resultsNames(dds)
resLFC <- lfcShrink(dds, coef=2, type="apeglm")

# an alternate analysis: likelihood ratio test
ddsLRT <- DESeq(dds, test="LRT", reduced= ~ 1)
resLRT <- results(ddsLRT)

 

3.5)DESeqDataSet-class(DESeqDataSet object and constructors)

 语法:

DESeqDataSet(se, design, ignoreRank = FALSE)
DESeqDataSetFromMatrix(countData, colData, design, tidy = FALSE,ignoreRank = FALSE, ...)
DESeqDataSetFromHTSeqCount(sampleTable, directory = ".", design,ignoreRank = FALSE, ...)
DESeqDataSetFromTximport(txi, colData, design, ...)

例子:

countData <- matrix(1:100,ncol=4)
condition <- factor(c("A","A","B","B"))
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), ~ condition)

 3.6)DESeqResults-class:DESeqResults object and constructor

语法:DESeqResults(DataFrame, priorInfo = list())

参数:

DataFrame:a DataFrame of results, standard column names are: baseMean, log2FoldChange,lfcSE, stat, pvalue, padj.
priorInfo:a list giving information on the log fold change prior

3.7)DESeqTransform-class(DESeqTransform object and constructor)

语法:DESeqTransform(SummarizedExperiment)

参数:SummarizedExperiment a RangedSummarizedExperiment

3.8)rlog Apply a ’regularized log’ transformation

用法:
rlog(object, blind = TRUE, intercept, betaPriorVar, fitType = "parametric")
rlogTransformation(object, blind = TRUE, intercept, betaPriorVar,fitType = "parametric")

dds <- makeExampleDESeqDataSet(m=6,betaSD=1)
rld <- rlog(dds)
dists <- dist(t(assay(rld)))
plot(hclust(dists))

 

3.9)plotPCA(Sample PCA plot for transformed data)

用法:plotPCA(object, intgroup = "condition",ntop = 500, returnData = FALSE)

参数:

object:a DESeqTransform object, with data in assay(x), produced for example by either rlog or varianceStabilizingTransformation.
intgroup: interesting groups: a character vector of names in colData(x) to use for grouping
ntop:number of top genes to use for principal components, selected by highest row variance
returnData:should the function only return the data.frame of PC1 and PC2 with intgroup covariates for custom plotting

# using rlog transformed data:
dds <- makeExampleDESeqDataSet(betaSD=1)
rld <- rlog(dds)
plotPCA(rld)

# also possible to perform custom transformation:
dds <- estimateSizeFactors(dds)
# shifted log of normalized counts
se <- SummarizedExperiment(log2(counts(dds, normalized=TRUE) + 1),
colData=colData(dds))
# the call to DESeqTransform() is needed to
# trigger our plotPCA method.
plotPCA( DESeqTransform( se ) )

 

3.10)

 

posted on 2018-09-13 17:12  发那个太丢人  阅读(1798)  评论(0编辑  收藏  举报