前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >表达差异基因分析

表达差异基因分析

原创
作者头像
爱学习的小明明
修改2020-09-23 10:15:21
1.4K0
修改2020-09-23 10:15:21
举报
文章被收录于专栏:R语言学习

1安装BiocManage,再安装DESeq2包

代码语言:javascript
复制
> # <差异基因分析>
> # 1.判断是否有BiocManager包,若不存在则安装
> options(repos=structure(c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))) #设置清华镜像,加速下载
> 
> if (!requireNamespace("BiocManager", quietly = TRUE))
+   install.packages("BiocManager")
> if (!requireNamespace("DESeq2", quietly = TRUE))
+   BiocManager::install('DESeq2')  #通过BiocManager安装DESeq2 
> library(DESeq2) #加载library

2输入数据

代码语言:javascript
复制
> #输入数据要求
> # DEseq2要求输入数据是由整数组成的矩阵
> # DESeq2要求矩阵是没有标准化的,一定记住用readcount
> 
> ##2.读入所有基因原始readscount表达矩阵,行为基因,列为样品
> A <- read.table(p, header = T, row.names = 1)
> B <- as.matrix(A) #转换成矩阵格式,保证都是数值

这里用的read_count
这里用的read_count

3实验分组信息

代码语言:javascript
复制
> coldata <- read.table("/home/shijm/Rlearning/R-Online-learning/data/sample_info.txt",header = T,row.names = 1)
> coldata <- coldata[, c("condition", "type")]

coldata:这里的行名必须和read_count的列明对应
coldata:这里的行名必须和read_count的列明对应

4.制作dds对象,构建差异基因分析所需的数据格式

代码语言:javascript
复制
> dds <- DESeqDataSetFromMatrix(countData = B, colData = coldata, design = ~ condition);
> 
> # countData = B,readscount矩阵
> # colData = coldata,分组信息,根据这个才能在2组之间比较
> # design = ~ condition,公式,表示按照condition进行分析

5.差异分析结果

代码语言:javascript
复制
> dds <- DESeq(dds)

6.提取结果,在treated和untreated组进行比较

代码语言:javascript
复制
res <- results(dds, contrast = c("condition", "treated", "untreated")) 
> resordered <- res[order(res$padj),]	 #按照矫正后的P-value padj从小到大排序
> sum(res$padj < 0.05, na.rm = TRUE)	#总结padj小于0.05显著差异的基因
[1] 356

7.输出图片

代码语言:javascript
复制
plotMA(res)	#画火山图,横轴是标准化后的平均readscount,纵轴是差异倍数,大于0是上调,小于0是下调,蓝色点表示显著差异的基因

8使用ggplot自己绘制火山图

代码语言:javascript
复制
> pvalue<-res$padj
> pvalue<-as.matrix(pvalue)
> fc<-res$log2FoldChange
> fc<-as.matrix(fc)
> df<-data.frame(pvalue,fc)
> df<-na.omit(df,cols=c("pvalue","fc"))   这里去掉pvalue或fc列有NA值得行
> df$threshold[df$pvalue < 0.05 & df$fc>0 ] = "up"
> df$threshold[df$pvalue < 0.05 & df$fc<0 ] = "down"
> df$threshold[df$pvalue > 0.05 & df$fc>=0|df$fc<=0 ] = "non"
>

构建的df用于火山图绘制
构建的df用于火山图绘制
代码语言:javascript
复制
> ggplot(df,aes(x=fc,y=-log10(pvalue)))+
+   xlab("log2 Fold Change")+ylab("-log10P-Value")+
+   
+   geom_point(aes(color=threshold),size=1,alpha=0.6,)+
+   
+   scale_color_manual(values=c("red","black","blue")) +
+   
+   geom_hline(aes(yintercept=-log10(0.05)),colour="grey",size=1.2 ,linetype=2) + #增加水平间隔线
+   
+   geom_vline(aes(xintercept=0), colour="grey",size=1.2 ,linetype=2)+ #增加垂直间隔线
+   
+   theme_few()+theme(legend.title = element_blank()) #去掉网格背景和图注标签

原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。

如有侵权,请联系 cloudcommunity@tencent.com 删除。

原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。

如有侵权,请联系 cloudcommunity@tencent.com 删除。

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1安装BiocManage,再安装DESeq2包
  • 2输入数据
  • 3实验分组信息
  • 4.制作dds对象,构建差异基因分析所需的数据格式
  • 5.差异分析结果
  • 6.提取结果,在treated和untreated组进行比较
  • 7.输出图片
  • 8使用ggplot自己绘制火山图
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档