搜索
查看: 9302|回复: 16

[mRNA-seq] 转录组入门(7-8)从差异基因到富集分析

[复制链接]

29

主题

131

帖子

1200

积分

金牌会员

Rank: 6Rank: 6

积分
1200
发表于 2017-8-14 20:19:59 | 显示全部楼层 |阅读模式
本帖最后由 anlan 于 2017-8-14 20:27 编辑

差异基因表达分析
我按照前面的流程转录组入门(1-6)从测序数据到生成count矩阵,将小鼠的4个样本又重新跑了一遍,从而获得一个新的count文件:mouse_all_count.txt,有需要的话,可以下载下来进行后续的差异分析。

一般来说,由于普遍认为高通量的read count符合泊松分布,所以一些差异分析的R包都是基于负二项式分布模型的,比如DESeq、DESeq2和edgeR等,所以这3个R包从整体上来说是类似的(但各自标准化算法是不一样的)。

当然还有一个常用的R包则是Limma包,其中的limma-trend和limma-voom都能用来处理RNA-Seq数据(但对应适用的情况不一样)

下面准备适用DESeq2和edgeR两个R包分别对小鼠的count数据进行差异表达分析,然后取两者的结果的交集的基因作为差异表达基因。

  • DEseq2
    [AppleScript] 纯文本查看 复制代码
    library(DESeq2)
                    ##数据预处理
                    database <- read.table(file = "mouse_all_count.txt", sep = "\t", header = T, row.names = 1)
                    database <- round(as.matrix(database))
                    
                    ##设置分组信息并构建dds对象
                    condition <- factor(c(rep("control",2),rep("Akap95",2)), levels = c("control", "Akap95"))
                    coldata <- data.frame(row.names = colnames(database), condition)
                    dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)
                    dds <- dds[ rowSums(counts(dds)) > 1, ]
                    
                    ##使用DESeq函数估计离散度,然后差异分析获得res对象
                    dds <- DESeq(dds)
                    res <- results(dds)
                    
                    #最后设定阈值,筛选差异基因,导出数据(全部数据。包括标准化后的count数)
                    res <- res[order(res$padj),]
                    diff_gene <- subset(res, padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
                    diff_gene <- row.names(diff_gene)
                    resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by="row.names",sort=FALSE)
                    write.csv(resdata,file = "control_vs_Akap95.csv",row.names = F)

    最终获得572个差异基因(筛选标准为padj < 0.05, |log2FoldChange| > 1)
  • edgeR
    [AppleScript] 纯文本查看 复制代码
    library(edgeR)
                    ##跟DESeq2一样,导入数据,预处理(用了cpm函数)
                    exprSet<- read.table(file = "mouse_all_count.txt", sep = "\t", header = TRUE, row.names = 1, stringsAsFactors = FALSE)
                    group_list <- factor(c(rep("control",2),rep("Akap95",2)))
                    exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2,]
                    
                    ##设置分组信息,并做TMM标准化
                    exprSet <- DGEList(counts = exprSet, group = group_list)
                    exprSet <- calcNormFactors(exprSet)
                    
                    ##使用qCML(quantile-adjusted conditional maximum likelihood)估计离散度(只针对单因素实验设计)
                    exprSet <- estimateCommonDisp(exprSet)
                    exprSet <- estimateTagwiseDisp(exprSet)
                    
                    ##寻找差异gene(这里的exactTest函数还是基于qCML并且只针对单因素实验设计),然后按照阈值进行筛选即可
                    et <- exactTest(exprSet)
                    tTag <- topTags(et, n=nrow(exprSet))
                    diff_gene_edgeR <- subset(tTag$table, FDR < 0.05 & (logFC > 1 | logFC < -1))
                    diff_gene_edgeR <- row.names(diff_gene_edgeR)
                    write.csv(tTag$table,file = "control_vs_Akap95_edgeR.csv")

    最终获得688个差异基因(筛选标准为FDR < 0.05, |log2FC| > 1)   
  • 取DESeq2和edgeR两者结果的交集
    [AppleScript] 纯文本查看 复制代码
    diff_gene <- diff_gene_DESeq2[diff_gene_DESeq2 %in% diff_gene_edgeR]

    最终的差异基因数目为545个
    [AppleScript] 纯文本查看 复制代码
    head(diff_gene)
                    [1] "ENSMUSG00000003309.14" "ENSMUSG00000046323.8"  "ENSMUSG00000001123.15"
                    [4] "ENSMUSG00000023906.2"  "ENSMUSG00000044279.15" "ENSMUSG00000018569.12"
  • 其他两个R包(DESeq和limma)就不在这尝试了,我之前做过对于这4个R包的简单使用笔记,可以参考下:
    简单使用DESeq做差异分析
    简单使用DESeq2/EdgeR做差异分析
    简单使用limma做差异分析

GO&&KEGG富集分析
以前一直没有机会用Y叔写的clusterProfiler包,这次正好看说明用一下。

  • GO富集,加载clusterProfiler包和org.Mm.eg.db包(小鼠嘛),然后将ENSEMBL ID后面的版本号去掉,不然后面不识别这个ID,然后按照clusterProfiler包的教程说明使用函数即可。
    [AppleScript] 纯文本查看 复制代码
    library(clusterProfiler)
                    library(org.Mm.eg.db)
                    
                    ##去除ID的版本号
                    diff_gene_ENSEMBL <- unlist(lapply(diff_gene, function(x){strsplit(x, "\\.")[[1]][1]}))
                    ##GOid mapping + GO富集
                    ego <- enrichGO(gene = diff_gene_ENSEMBL, OrgDb = org.Mm.eg.db, 
                                    keytype = "ENSEMBL", ont = "BP", pAdjustMethod = "BH", 
                                    pvalueCutoff = 0.01, qvalueCutoff = 0.05)
                    ##查看富集结果数据
                    enrich_go <- as.data.frame(ego)
                    ##作图
                    barplot(ego, showCategory=10)
                    dotplot(ego)
                    enrichMap(ego)
                    plotGOgraph(ego)
  • KEGG富集,首先需要将ENSEMBL ID转化为ENTREZ ID,然后使用ENTREZ ID作为kegg id,从而通过enrichKEGG函数从online KEGG上抓取信息,并做富集
    [AppleScript] 纯文本查看 复制代码
    library(clusterProfiler)
                    library(org.Mm.eg.db)
    
                    ##ID转化
                    ids <- bitr(diff_gene_ENSEMBL, fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = "org.Mm.eg.db")
                    kk <- enrichKEGG(gene = ids[,2], organism = "mmu", keyType = "kegg",
                     pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.1)
                    ##查看富集结果数据
                    enrich_kegg <- as.data.frame(kk)
                    ##作图
                    dotplot(kk)
    

到这里为止,转录组的差异表达分析算是做完了,简单的来说,这个过程就是将reads mapping 到reference上,然后计数获得count数,然后做差异分析,最后来个GO KEGG,over了。。。

对于mapping和计数这两部还有其实还有好多软件,具体可见文献:Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis,有时间可以都尝试下。

至于GO && KEGG这两步,对于人、小鼠等模式物种来说,不考虑方便因素来说,完全可以自己写脚本来完成,数据可以从gene ontology官网下载,然后就是GO id与gene id相互转化了。KEGG 也是一样,也可以用脚本去KEGG网站上自行抓取,先知道URL规则,然后爬数据即可。

附上有道笔记链接:http://note.youdao.com/noteshare ... 4acb2e42a0a2885b181




上一篇:转录组入门——5
下一篇:泛基因组测序(学习与探究)
回复

使用道具 举报

7

主题

37

帖子

412

积分

中级会员

Rank: 3Rank: 3

积分
412
发表于 2017-12-1 10:07:00 | 显示全部楼层
在GO富集上, 用你的代码怎么会出现  
Error in enrichGO(gene = diff_gene_ENSEMBL, OrgDb = org.Mm.eg.db, keytype = "ENSEMBL",  :
  unused argument (keytype = "ENSEMBL")

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1200

积分

金牌会员

Rank: 6Rank: 6

积分
1200
 楼主| 发表于 2017-12-1 15:05:16 | 显示全部楼层
xiao7462 发表于 2017-12-1 10:07
在GO富集上, 用你的代码怎么会出现  
Error in enrichGO(gene = diff_gene_ENSEMBL, OrgDb = org.Mm.eg.db ...

我拿我之前的count数据的前200行gene作为输入,做了下enrichGO,并没有发现报错,出的结果也是正常的,要不可能是你的代码文字的格式问题?


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

7

主题

37

帖子

412

积分

中级会员

Rank: 3Rank: 3

积分
412
发表于 2017-12-1 19:07:23 | 显示全部楼层
照搬您的代码 ,可是到diff_gene这里就是 null了

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1200

积分

金牌会员

Rank: 6Rank: 6

积分
1200
 楼主| 发表于 2017-12-1 21:29:02 | 显示全部楼层
xiao7462 发表于 2017-12-1 19:07
照搬您的代码 ,可是到diff_gene这里就是 null了

0-0 我粗略看了下,应该是data$X这里出问题,因为我输入的mouse_all_count.txt的第一列是没有设置列名的,所以R默认给了X作为列名,所以这里的X是大写的,我看你的代码那里是小写的,所以diff_gene是空值,你看看是不是这个原因
回复 支持 反对

使用道具 举报

7

主题

37

帖子

412

积分

中级会员

Rank: 3Rank: 3

积分
412
发表于 2017-12-1 22:07:03 | 显示全部楼层
得出结果还是error

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1200

积分

金牌会员

Rank: 6Rank: 6

积分
1200
 楼主| 发表于 2017-12-2 20:31:25 | 显示全部楼层
xiao7462 发表于 2017-12-1 22:07
得出结果还是error

我也不知道了,可能是包没装好?你可以试试重新安装下包


或者你去https://bioconductor.org/package ... lusterProfiler.html网站
按照作者给的测试的基因,试一下网站上的enrichGO函数,看看会不会报错,如果正常出结果的话,那么这个包是没问题的,不然的话就重新安装包直到不会报错为止吧。。
回复 支持 反对

使用道具 举报

7

主题

37

帖子

412

积分

中级会员

Rank: 3Rank: 3

积分
412
发表于 2017-12-3 13:55:57 | 显示全部楼层
anlan 发表于 2017-12-2 20:31
我也不知道了,可能是包没装好?你可以试试重新安装下包

是不是你的转GENE ID的 代码 省略了 没有写出来 ?
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1200

积分

金牌会员

Rank: 6Rank: 6

积分
1200
 楼主| 发表于 2017-12-3 15:27:11 | 显示全部楼层
xiao7462 发表于 2017-12-3 13:55
是不是你的转GENE ID的 代码 省略了 没有写出来 ?

用这个包 做GO分析,如果是Ensembl id的话,我都不转化ID的,因为本来这个包就支持ensembl,你可以看看说明文档,里面写的了。。。。你用entrez id就不报错了?
回复 支持 反对

使用道具 举报

7

主题

37

帖子

412

积分

中级会员

Rank: 3Rank: 3

积分
412
发表于 2017-12-3 19:20:20 | 显示全部楼层
加了这个  
[AppleScript] 纯文本查看 复制代码
gene=bitr(diff_gene_ENSEMBL,fromType = 'ENSEMBL',toType = 'ENTREZID',OrgDb = 'org.Mm.eg.db')


没有报错了 ,但是得不到结果
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2018-11-16 22:49 , Processed in 0.037159 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.