DESeq2转录组差异表达分析实例

用户7010445

发布于 2020-03-03 15:13:53

2.2K0

发布于 2020-03-03 15:13:53

文章被收录于专栏：小明的数据分析笔记本小明的数据分析笔记本

参考文章

生信技能树B站转录组数据分析视频
https://github.com/jmzeng1314/my-R/blob/master/8-DEG/example_input_output/DESeq2/github-R-code.txt
https://github.com/jmzeng1314/GEO/blob/master/task1-check-specific-genes/step3-DEG.R
DESeq2帮助文档

我的R语言版本是3.6.1

安装分析过程需要用的的R包

DESeq2 差异表达分析

BiocManager::install("DESeq2")

使用library(DESeq2)加载的时候遇到报错：载入了名字空间‘rlang’ 0.4.0，但需要的是>= 0.4.2 解决办法：将rlang包手动删除，rlang所在的路径是\R-3.6.1\library\rlang。然后使用命令install.packages("rlang")重新安装就可以了

pasilla 使用这个R包中的数据

BiocManager::install("pasilla")

读入数据

library(pasilla)
pasCts<-system.file("extdata",
                    "pasilla_gene_counts.tsv",
                    package="pasilla",
                    mustWork = T)
pasCts
pasAnno<-system.file("extdata",
                     "pasilla_sample_annotation.csv",
                     package="pasilla",
                     mustWork=T)
pasAnno
df<-read.csv(pasCts,sep="\t",row.names = "gene_id")
cts<-as.matrix(df)
head(cts)
coldata<-read.csv(pasAnno,row.names = 1)
head(coldata)
coldata<-coldata[,c("condition","type")]
rownames(coldata)<-sub("fb","",rownames(coldata))
all(rownames(coldata) %in% colnames(cts))
cts<-cts[,rownames(coldata)]
all(rownames(coldata) == colnames(cts))

获得了两个数据集

> coldata
           condition        type
treated1     treated single-read
treated2     treated  paired-end
treated3     treated  paired-end
untreated1 untreated single-read
untreated2 untreated single-read
untreated3 untreated  paired-end
untreated4 untreated  paired-end
> head(cts)
            treated1 treated2 treated3 untreated1 untreated2 untreated3
FBgn0000003        0        0        1          0          0          0
FBgn0000008      140       88       70         92        161         76
FBgn0000014        4        0        0          5          1          0
FBgn0000015        1        0        0          0          2          1
FBgn0000017     6205     3072     3334       4664       8714       3564
FBgn0000018      722      299      308        583        761        245
            untreated4
FBgn0000003          0
FBgn0000008         70
FBgn0000014          0
FBgn0000015          2
FBgn0000017       3150
FBgn0000018        310

cts是表达矩阵
coldata是用来指定样本分组的数据集

DESeq2差异表达分析

library(DESeq2)
dds<-DESeqDataSetFromMatrix(countData = cts,
                            colData = coldata,
                            design= ~ condition)
dds
keep<-rowSums(counts(dds))>=10
dds<-dds[keep,]
dds<-DESeq(dds)
res<-results(dds,contrast=c("condition","treated","untreated"))

火山图

resdata<-merge(as.data.frame(res),
               as.data.frame(counts(dds,normalized=T)),
               by="row.names",sort=F)
DEG<-resdata
logFC_cutoff<-with(DEG,mean(abs(log2FoldChange))+2*sd(abs(log2FoldChange)))
logFC_cutoff
DEG$change<-as.factor(ifelse(DEG$pvalue<0.05&abs(DEG$log2FoldChange)>logFC_cutoff,
                             ifelse(DEG$log2FoldChange>logFC_cutoff,"UP","DOWN"),
                             "NOT"))
this_title <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',]))
DEG<-na.omit(DEG)
library(ggplot2)
ggplot(data=DEG,aes(x=log2FoldChange,y=-log10(pvalue),color=change))+
  geom_point(alpha=0.4,size=1.75)+
  labs(x="log2 fold change")+ ylab("-log10 pvalue")+
  ggtitle(this_title)+theme_bw(base_size = 20)+
  theme(plot.title = element_text(size=15,hjust=0.5))+
  scale_color_manual(values=c('blue','black','red'))