RNA-seq(7): DEseq2筛选差异表达基因并注释(bioMart)

============================================
写在前面:可以参考另外一篇《得到差异基因后怎么做?
============================================

  • 走到这一步,似乎顺畅了许多,最主要时间不用花费那么多了。另外,以前曾经处理过不计其数的芯片,挑选差异表达基因,筛选关键基因,功能富集,还有基于全部数据的WGCNA(当然你也可以用差异基因来做,虽然不推荐,看不少文章也这么发),GSEA,PPI等等,这些后续我会慢慢发出来。
  • 但是,因为以前处理的芯片表达谱数据是符合正态分布,所以可以用t检验来筛选差异表达基因,但RNA-seq的read count普遍认为符合泊松分布。所以筛选DEGs的方法还是不一样
------------要求---------------
  • 载入表达矩阵
  • 设置好分组信息
  • 用DEseq2进行差异分析
  • 输出差异分析结果
    来源于生信技能树
-------------------------参考文章-----------------------------------

-------------------开始------------------------

  • 正式载入数据之前,先看一下数据结构data structure,因为这是我们下面一切分析的起点。
  • 我们需要准备2个table:一个是countData,一个是colData
two kinds of data.png

关于上面两个表的说明

  • countData表示的是count矩阵,行代表gene,列代表样品,中间的数字代表对应count数。colData表示sample的元数据,因为这个表提供了sample的元数据。
    because this table supplies metadata/information about the columns of the countData matrix. Notice that the first column of colData must match the column names of countData (except the first, which is the gene ID column).
  • colData的行名与countData的列名一致(除去代表gene ID的第一列)

1 载入数据(countData和colData)

> library(tidyverse)
> library(DESeq2)
> #import data
> setwd("F:/rna_seq/data/matrix")
> mycounts<-read.csv("readcount.csv")
> head(mycounts)
                   X control1 control2 treat1 treat2
1 ENSMUSG00000000001     1648     2306   2941   2780
2 ENSMUSG00000000003        0        0      0      0
3 ENSMUSG00000000028      835      950   1366   1051
4 ENSMUSG00000000031       65       83     52     53
5 ENSMUSG00000000037       70       53     94     66
6 ENSMUSG00000000049        0        3      4      5
#这里有个x,需要去除,先把第一列当作行名来处理
> rownames(mycounts)<-mycounts[,1]
#把带X的列删除
> mycounts<-mycounts[,-1]
> head(mycounts)
                   control1 control2 treat1 treat2
ENSMUSG00000000001     1648     2306   2941   2780
ENSMUSG00000000003        0        0      0      0
ENSMUSG00000000028      835      950   1366   1051
ENSMUSG00000000031       65       83     52     53
ENSMUSG00000000037       70       53     94     66
ENSMUSG00000000049        0        3      4      5
# 这一步很关键,要明白condition这里是因子,不是样本名称;小鼠数据有对照组和处理组,各两个重复
> condition <- factor(c(rep("control",2),rep("treat",2)), levels = c("control","treat"))
> condition
[1] control control treat   treat  
Levels: control treat
#colData也可以自己在excel做好另存为.csv格式,再导入即可
> colData <- data.frame(row.names=colnames(mycounts), condition)
> colData
         condition
control1   control
control2   control
treat1       treat
treat2       treat

2构建dds对象,开始DESeq流程

注释:dds=DESeqDataSet Object
> dds <- DESeqDataSetFromMatrix(mycounts, colData, design= ~ condition)
> dds <- DESeq(dds)
> # 查看一下dds的内容
> dds

显示为

class: DESeqDataSet 
dim: 6 4 
metadata(1): version
assays(3): counts mu cooks
rownames(6): ENSMUSG00000000001 ENSMUSG00000000003 ... ENSMUSG00000000037 ENSMUSG00000000049
rowData names(21): baseMean baseVar ... deviance maxCooks
colnames(4): control1 control2 treat1 treat2
colData names(2): condition sizeFactor

3 总体结果查看

接下来,我们要查看treat versus control的总体结果,并根据p-value进行重新排序。利用summary命令统计显示一共多少个genes上调和下调(FDR0.1)

> res = results(dds, contrast=c("condition", "control", "treat"))
#或下面命令
> res= results(dds)
> res = res[order(res$pvalue),]
> head(res)
> summary(res)
#所有结果先进行输出
> write.csv(res,file="All_results.csv")
> table(res$padj<0.05)
上述代码的结果显示
> res = results(dds2, contrast=c("condition", "control", "treat"))
> res = res[order(res$pvalue),]
> head(res)
log2 fold change (MLE): condition control vs treat 
Wald test p-value: condition control vs treat 
DataFrame with 6 rows and 6 columns
                    baseMean log2FoldChange     lfcSE      stat       pvalue         padj
                   <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>
ENSMUSG00000003309  548.1926       3.231611 0.2658125 12.157485 5.234568e-34 8.193146e-30
ENSMUSG00000046323  404.1894       3.067050 0.2628220 11.669687 1.820923e-31 1.425055e-27
ENSMUSG00000001123  341.8542       2.797485 0.2766499 10.112004 4.887441e-24 2.549941e-20
ENSMUSG00000023906  951.9460       2.382307 0.2510718  9.488551 2.342684e-21 9.116395e-18
ENSMUSG00000018569  485.4839       3.136031 0.3312999  9.465836 2.912214e-21 9.116395e-18
ENSMUSG00000000184  601.0842      -2.827750 0.3154171 -8.965112 3.099648e-19 8.085948e-16
> summary(res)

out of 28335 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)     : 445, 1.6% 
LFC < 0 (down)   : 625, 2.2% 
outliers [1]     : 0, 0% 
low counts [2]   : 12683, 45% 
(mean count < 18)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

> write.csv(res,file="All_results.csv")
> table(res$padj<0.05)

FALSE  TRUE 
14909   743
可见,一共445个genes上调,625个genes下调,没有离群值。padj小于0.05的共有743个。

4 提取差异表达genes(DEGs)并进行gene symbol注释

差异表达基因的界定很不统一,但log2FC是用的最广泛同时也是最不精确的方式,但因为其好理解所以广泛被应用尤其芯片数据处理中,记的是havard universit做过一个统计,FC=2相对比较可靠。但无论怎么,这种界定人为因素太大,过于武断。所以GSEA,WGCNA是拿全部表达数据(可以进行初步过滤)来进行分析,并且WGCNA采取软阈值的方式来挑选genes更加合理。既然挑选差异表达基因,还是采取log2FC和padj来进行。

获取padj(p值经过多重校验校正后的值)小于0.05,表达倍数取以2为对数后大于1或者小于-1的差异表达基因。代码如下
> diff_gene_deseq2 <-subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
#或
#> diff_gene_deseq2 <-subset(res,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
> dim(diff_gene_deseq2)
> head(diff_gene_deseq2)
> write.csv(diff_gene_deseq2,file= "DEG_treat_vs_control.csv")

结果显示如下:

> diff_gene_deseq2 <-subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
> dim(diff_gene_deseq2)
[1] 431   6
> head(diff_gene_deseq2)
log2 fold change (MLE): condition control vs treat 
Wald test p-value: condition control vs treat 
DataFrame with 6 rows and 6 columns
                    baseMean log2FoldChange     lfcSE      stat       pvalue         padj
                   <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>
ENSMUSG00000003309  548.1926       3.231611 0.2658125 12.157485 5.234568e-34 8.193146e-30
ENSMUSG00000046323  404.1894       3.067050 0.2628220 11.669687 1.820923e-31 1.425055e-27
ENSMUSG00000001123  341.8542       2.797485 0.2766499 10.112004 4.887441e-24 2.549941e-20
ENSMUSG00000023906  951.9460       2.382307 0.2510718  9.488551 2.342684e-21 9.116395e-18
ENSMUSG00000018569  485.4839       3.136031 0.3312999  9.465836 2.912214e-21 9.116395e-18
ENSMUSG00000000184  601.0842      -2.827750 0.3154171 -8.965112 3.099648e-19 8.085948e-16

5 用bioMart对差异表达基因进行注释

RNA-seq(6): reads计数,合并矩阵并进行注释代码一样
library('biomaRt')
library("curl")
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
my_ensembl_gene_id<-row.names(diff_gene_deseq2)
#listAttributes(mart)
mms_symbols<- getBM(attributes=c('ensembl_gene_id','external_gene_name',"description"),
                    filters = 'ensembl_gene_id', values = my_ensembl_gene_id, mart = mart)
结果如下:
> library('biomaRt')
> library("curl")
> mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
> my_ensembl_gene_id<-row.names(diff_gene_deseq2)
> mms_symbols<- getBM(attributes=c('ensembl_gene_id','external_gene_name',"description"),
+                     filters = 'ensembl_gene_id', values = my_ensembl_gene_id, mart = mart)
> head(mms_symbols)
     ensembl_gene_id external_gene_name                                                                                  description
1 ENSMUSG00000000120               Ngfr nerve growth factor receptor (TNFR superfamily, member 16) [Source:MGI Symbol;Acc:MGI:97323]
2 ENSMUSG00000000184              Ccnd2                                                  cyclin D2 [Source:MGI Symbol;Acc:MGI:88314]
3 ENSMUSG00000000276               Dgke                           diacylglycerol kinase, epsilon [Source:MGI Symbol;Acc:MGI:1889276]
4 ENSMUSG00000000308              Ckmt1               creatine kinase, mitochondrial 1, ubiquitous [Source:MGI Symbol;Acc:MGI:99441]
5 ENSMUSG00000000320             Alox12                               arachidonate 12-lipoxygenase [Source:MGI Symbol;Acc:MGI:87998]
6 ENSMUSG00000000708              Kat2b                           K(lysine) acetyltransferase 2B [Source:MGI Symbol;Acc:MGI:1343094]

5合并数据:res结果+mms_symbols合并成一个文件

合并的话两个数据必须有共同的列名,我们先看一下

> head(diff_gene_deseq2)
log2 fold change (MLE): condition control vs treat 
Wald test p-value: condition control vs treat 
DataFrame with 6 rows and 6 columns
                    baseMean log2FoldChange     lfcSE      stat       pvalue         padj
                   <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>
ENSMUSG00000003309  548.1926       3.231611 0.2658125 12.157485 5.234568e-34 8.193146e-30
ENSMUSG00000046323  404.1894       3.067050 0.2628220 11.669687 1.820923e-31 1.425055e-27
ENSMUSG00000001123  341.8542       2.797485 0.2766499 10.112004 4.887441e-24 2.549941e-20
ENSMUSG00000023906  951.9460       2.382307 0.2510718  9.488551 2.342684e-21 9.116395e-18
ENSMUSG00000018569  485.4839       3.136031 0.3312999  9.465836 2.912214e-21 9.116395e-18
ENSMUSG00000000184  601.0842      -2.827750 0.3154171 -8.965112 3.099648e-19 8.085948e-16
> head(mms_symbols)
     ensembl_gene_id external_gene_name
1 ENSMUSG00000000120               Ngfr
2 ENSMUSG00000000184              Ccnd2
3 ENSMUSG00000000276               Dgke
4 ENSMUSG00000000308              Ckmt1
5 ENSMUSG00000000320             Alox12
6 ENSMUSG00000000708              Kat2b

可见,两个文件没有共同的列名,所以要先给'diff_gene_deseq2'添加一个‘ensembl_gene_id’的列名,方法如下:(应该有更简便的方法)

> ensembl_gene_id<-rownames(diff_gene_deseq2)
> diff_gene_deseq2<-cbind(ensembl_gene_id,diff_gene_deseq2)
> colnames(diff_gene_deseq2)[1]<-c("ensembl_gene_id")
> diff_name<-merge(diff_gene_deseq2,mms_symbols,by="ensembl_gene_id")
>head(diff_name)
#查看Akap8的情况
Akap8 <- diff_name[diff_name$external_gene_name=="Akap8",]
中间显示过程如下:
> ensembl_gene_id<-rownames(diff_gene_deseq2)
> diff_gene_deseq2<-cbind(ensembl_gene_id,diff_gene_deseq2)
> colnames(diff_gene_deseq2)[1]<-c("ensembl_gene_id")
> diff_name<-merge(diff_gene_deseq2,mms_symbols,by="ensembl_gene_id")
> head(diff_name)
DataFrame with 6 rows and 9 columns
     ensembl_gene_id   baseMean log2FoldChange     lfcSE      stat       pvalue         padj external_gene_name
         <character>  <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>        <character>
1 ENSMUSG00000000120  434.04177      -1.293648 0.2713146 -4.768072 1.859970e-06 2.064698e-04               Ngfr
2 ENSMUSG00000000184  601.08417      -2.827750 0.3154171 -8.965112 3.099648e-19 8.085948e-16              Ccnd2
3 ENSMUSG00000000276  668.12500      -1.071362 0.2445381 -4.381168 1.180446e-05 9.603578e-04               Dgke
4 ENSMUSG00000000308  207.46719       1.944949 0.3427531  5.674489 1.391035e-08 3.819733e-06              Ckmt1
5 ENSMUSG00000000320   61.96266       1.451927 0.4637101  3.131109 1.741473e-03 4.105051e-02             Alox12
6 ENSMUSG00000000708 1070.03203      -1.046546 0.2049062 -5.107440 3.265530e-07 5.056107e-05              Kat2b
                                                                                   description
                                                                                   <character>
1 nerve growth factor receptor (TNFR superfamily, member 16) [Source:MGI Symbol;Acc:MGI:97323]
2                                                  cyclin D2 [Source:MGI Symbol;Acc:MGI:88314]
3                           diacylglycerol kinase, epsilon [Source:MGI Symbol;Acc:MGI:1889276]
4               creatine kinase, mitochondrial 1, ubiquitous [Source:MGI Symbol;Acc:MGI:99441]
5                               arachidonate 12-lipoxygenase [Source:MGI Symbol;Acc:MGI:87998]
6                           K(lysine) acetyltransferase 2B [Source:MGI Symbol;Acc:MGI:1343094]
> Akap8 <- diff_name[diff_name$external_gene_name=="Akap8",]
> Akap8
DataFrame with 1 row and 9 columns
     ensembl_gene_id  baseMean log2FoldChange     lfcSE      stat       pvalue         padj external_gene_name
         <character> <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>        <character>
1 ENSMUSG00000024045  2281.053      -1.276329 0.2428315 -5.256028 1.471996e-07 2.775865e-05              Akap8
                                                           description
                                                           <character>
1 A kinase (PRKA) anchor protein 8 [Source:MGI Symbol;Acc:MGI:1928488]
至此,差异表达基因提取并注释完毕,下一步
  • 先进行数据可视化(Data visulization)
  • 然后进行富集分分析及可视化

后记:也可以用Y叔的ClusterProfiler或Annotation包进行基因名转换,很方便。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 214,313评论 6 496
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,369评论 3 389
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 159,916评论 0 349
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,333评论 1 288
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,425评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,481评论 1 292
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,491评论 3 412
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,268评论 0 269
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,719评论 1 307
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,004评论 2 328
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,179评论 1 342
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,832评论 4 337
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,510评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,153评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,402评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,045评论 2 365
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,071评论 2 352

推荐阅读更多精彩内容