搜索
查看: 6458|回复: 4

DESeq2包的系统学习

[复制链接]

13

主题

33

帖子

249

积分

版主

Rank: 7Rank: 7Rank: 7

积分
249
发表于 2017-1-17 18:41:42 | 显示全部楼层 |阅读模式
DESeq2包的系统学习
DESeq2包主要是针对RNA-seq的数据进行分析的,但还是想在这个版块发下,以便让大家与limma包对比下。。。
由于需要还是对DESeq2包进行了系统的学习,应用里面的测试数据,得到的提升是,知道了怎么做多因子的比较,代码的学习也了解的更加清楚,
过程具体如下:
[AppleScript] 纯文本查看 复制代码
# ###Typically, we recommend users to run samples from all groups together, and then
# use thecontrast argument of the results function to extract comparisons of
# interest after fitting the model using DESeq.
setwd("C:/Users/lixia/Desktop/learn/DESeq2/")
rm(list=ls())
library(DESeq2)
########
source("http://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("airway")
biocLite("IHW")
library(airway)
data("airway")
se <- airway
se
ddsSE <- DESeqDataSet(se, design = ~ cell + dex)
ddsSE
library("pasilla")
pasCts <- system.file("extdata", "pasilla_gene_counts.tsv",
                      package="pasilla", mustWork=TRUE)
pasAnno <- system.file("extdata", "pasilla_sample_annotation.csv",
                       package="pasilla", mustWork=TRUE)
countData <- as.matrix(read.csv(pasCts,sep="\t",row.names="gene_id"))  
###表达矩阵,如countData中相关reads数没有,记作0,不是NA,切记!!!
colData <- read.csv(pasAnno, row.names=1)
colData <- colData[,c("condition","type")]  ###select two column   ###分组信息
head(countData)
head(colData)
rownames(colData) <- sub("fb","",rownames(colData)) ###去除rownames(colData)中的“fb”,以匹配header的信息。。。
all(rownames(colData) %in% colnames(countData))     ###判断两者是否相等
######表达矩阵的列(样品),样本分组信息的行,必须一致
countData <- countData[, rownames(colData)]
all(rownames(colData) == colnames(countData))
dds <- DESeqDataSetFromMatrix(countData = countData,
                              colData = colData,
                              design = ~ condition)                         
dds
featureData <- data.frame(gene=rownames(countData))  ####添加一些额外的信息。。。gene
(mcols(dds) <- DataFrame(mcols(dds), featureData))  ###mcols,...???

dds <- DESeq(dds)  ##第二步,直接用DESeq函数即可  此步的目的!!!
dds
res <- results(dds)
res
write.table(res, file="ddseq2.xls", sep="\t",quote=F)
library("BiocParallel")
register(MulticoreParam(4)) 
# ##???
#These steps should take less than 30 seconds for most analyses. For experiments
# with many samples (e.g. 100 samples), one can take advantage of parallelized
# computation. Both of the above functions have an argument parallel which
# if set to TRUE can be used to distribute computation across cores specified by
# the register function of BiocParallel. For example, the following chunk (not
# evaluated here), would register 4 cores, and then the two functions above, with
# parallel=TRUE, would split computation over these cores.
res <- results(dds)
resOrdered <- res[order(res$padj),]
summary(res)
sum(res$padj < 0.1, na.rm=TRUE)
res05 <- results(dds, alpha=0.05)  ###alpha默认值(结合了P值+FDR)为0.1,如果需要更改则需要将其另外赋值,与上述比较
summary(res05)
sum(res05$padj < 0.05, na.rm=TRUE)
library("IHW")
resIHW <- results(dds, filterFun=ihw)
summary(resIHW)
###A new Bioconductor package,IHW, is now available that implements the
#method ofIndependent Hypothesis Weighting
sum(resIHW$padj < 0.1, na.rm=TRUE)
metadata(resIHW)$ihwResult
######
plotMA(res, main="DESeq2", ylim=c(-2,2))  ###log2FC,normalize的图 具体作用是啥???
###
idx <- identify(res$baseMean, res$log2FoldChange)
rownames(res)[idx]
###????

plotCounts(dds, gene=which.min(res$padj), intgroup="condition")

d <- plotCounts(dds, gene=which.min(res$padj), intgroup="condition",
                returnData=TRUE)
###d为啥?? condition啥时候引入的
library("ggplot2")
ggplot(d, aes(x=condition, y=count)) +
  geom_point(position=position_jitter(w=0.1,h=0)) +
  scale_y_log10(breaks=c(25,100,400))
mcols(res)$description
mcols(res)
res
write.csv(as.data.frame(resOrdered),file="condition_treated_results.csv")
resSig <- subset(resOrdered, padj < 0.1)
resSig
write.csv(as.data.frame(resSig),file="condition_treated_results resSig.csv")
######上述为一个因子的比较
### design = ~ condition
# dds <- DESeq(dds)  
# dds
# res <- results(dds)
# res
###下面讲解多个因子的比较
ddsMF <- dds
design(ddsMF) <- formula(~ type + condition)  
###不知道如何组合的,2与2的组合应该有四种。。。
##解答,默认情况下是后面的(本例中为conditon)进行比较,
##除非额外指定具体的contrast进行比较,
##如resMFType <- results(ddsMF,contrast=c("type", "single-read", "paired-end"))
ddsMF <- DESeq(ddsMF)
colData(dds)
##res的值是
head(res)
##condition untreated vs treated
resMF <- results(ddsMF)
head(resMF)
### 注意同head(res)输出结果的比较,是有差别的;condition untreated vs treated
resMFType <- results(ddsMF,contrast=c("type", "single-read", "paired-end"))
head(resMFType)
###transformations of the form  
##y= log2(n+ 1) or more generally, y= log2(n+n0);
#不知道两者的区别!!!

还需具体实战!!!






上一篇:缺失值和重复值的处理办法
下一篇:线性回归中“回归”的含义
回复

使用道具 举报

633

主题

1189

帖子

4054

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4054
发表于 2017-1-18 11:04:45 | 显示全部楼层
请把代码高亮一下,这样对读者友好,谢谢
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

4

主题

19

帖子

179

积分

注册会员

Rank: 2

积分
179
发表于 2017-6-15 17:42:24 | 显示全部楼层
我的数据不是多因子,而是五个组织,每个组织三个重复。怎么得到两两比较后的差异基因呢?一直没搞懂
回复 支持 反对

使用道具 举报

633

主题

1189

帖子

4054

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4054
发表于 2017-6-16 14:43:04 | 显示全部楼层
zhongguoren2016 发表于 2017-6-15 17:42
我的数据不是多因子,而是五个组织,每个组织三个重复。怎么得到两两比较后的差异基因呢?一直没搞懂 ...

makeContrasts
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

4

主题

19

帖子

179

积分

注册会员

Rank: 2

积分
179
发表于 2017-6-16 16:07:28 | 显示全部楼层
本帖最后由 zhongguoren2016 于 2017-6-16 16:08 编辑

谢谢Jimmy,这个问题我上午刚刚搞明白了。现在想问一下关于result函数结果文件中第二列数据为baseMean,应该指的是均一化后该基因在两个条件下的所有样本中的平均count值吧?如果我想知道该基因在各个处理中的均一化count值,是不是还要单独再计算出来呢?我觉得这点DEseq2就没有DEseq好,DEseq直接给出了每个条件下的均一化后的count值。
result结果:
3[MDUP1RW5P285QYUTXB$_2.png

DEseq结果:(列名向左多移了一列)
5IM`Q8@E2MIY@Z8KL_{D2.png
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2020-6-1 14:18 , Processed in 0.035032 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.