文章/答案/技术大牛

发布

社区首页 >专栏 >拟南芥的基因ID批量转换？差异基因，GO/KEGG数据库注释（转录组直接送你全套流程）

拟南芥的基因ID批量转换？差异基因，GO/KEGG数据库注释（转录组直接送你全套流程）

生信技能树

发布于 2019-07-04 08:55:35

8.2K00

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

新手遇到的问题都是类似的，比如批量ID转换

虽然我写过大量的教程：ID转换大全 不过都需要R基础，因为是大批量转换啊！

但热心肠的植物生物信息学教学大佬还是友善的给出了解决方案

我也狗尾续貂制作了一个网页工具教程：

简单的3个步骤，不会代码也可以很容易把ID批量转换啦！

不过有趣的是我搜索电脑资料，看到了2年前写的拟南芥教程。

不过我为什么会花时间写拟南芥教程呢？

1 首先加载必要的包

library(ggplot2)
library(stringr)
# source("https://bioconductor.org/biocLite.R")
# biocLite("clusterProfiler")
# biocLite("org.At.tair.db")
library(clusterProfiler)
library(org.At.tair.db)

2 然后导入数据

deg1=read.table('https://raw.githubusercontent.com/jmzeng1314/my-R/master/DEG_scripts/tair/DESeq2_DEG.Day1-Day0.txt')
deg1=na.omit(deg1)
head(deg1)

##              baseMean log2FoldChange    lfcSE      stat       pvalue
## AT3G01430   22.908929      18.989990 2.148261  8.839704 9.597263e-19
## AT1G47405   20.709551      20.958806 2.434574  8.608820 7.381677e-18
## AT2G33830 1938.159722      -2.560609 0.312663 -8.189678 2.619266e-16
## AT5G28627    8.118376     -21.131318 2.875691 -7.348257 2.008078e-13
## AT2G33750    9.789740     -19.989301 2.847513 -7.019915 2.220033e-12
## AT3G54500 2238.314652       2.720430 0.386305  7.042182 1.892517e-12
##                   padj
## AT3G01430 1.858318e-14
## AT1G47405 7.146571e-14
## AT2G33830 1.690562e-12
## AT5G28627 9.720602e-10
## AT2G33750 7.164418e-09
## AT3G54500 7.164418e-09

prefix='Day1-Day0'

3 然后判断显著差异基因

很明显，这个时候用padj来挑选差异基因即可，不需要看foldchange了。

table(deg1$padj<0.05)

## 
## FALSE  TRUE 
## 19166   197

table(deg1$padj<0.01)

## 
## FALSE  TRUE 
## 19303    60

diff_geneList = rownames(deg1[deg1$padj<0.05,])
up_geneList = rownames(deg1[deg1$padj<0.05 & deg1$log2FoldChange >0,])
down_geneList = rownames(deg1[deg1$padj<0.05 & deg1$log2FoldChange <0,])
length(diff_geneList)

## [1] 197

length(up_geneList)

## [1] 89

length(down_geneList)

## [1] 108

3.1 画个火山图看看挑选的差异基因合理与否

colnames(deg1)

## [1] "baseMean"       "log2FoldChange" "lfcSE"          "stat"          
## [5] "pvalue"         "padj"

log2FoldChange_Cutof = 0
## 这里我不准备用log2FoldChange来挑选差异基因，仅仅是用padj即可
deg1$change = as.factor(ifelse(deg1$padj < 0.05 & 
                                       abs(deg1$log2FoldChange) > log2FoldChange_Cutof,
                                     ifelse(deg1$log2FoldChange > log2FoldChange_Cutof ,'UP','DOWN'),'NOT'))
this_tile <- paste0('Cutoff for log2FoldChange is ',round(log2FoldChange_Cutof,3),
                    '\nThe number of up gene is ',nrow(deg1[deg1$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(deg1[deg1$change =='DOWN',])
)
g_volcano = ggplot(data=deg1, aes(x=log2FoldChange, y=-log10(padj), color=change)) +
  geom_point(alpha=0.4, size=1.75) +
  theme_set(theme_set(theme_bw(base_size=20)))+
  xlab("log2 fold change") + ylab("-log10 p-value") +
  ggtitle( this_tile  ) + 
  theme(plot.title = element_text(size=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red'))  ## corresponding to the levels(res$change)
print(g_volcano)

4 GO/KEGG注释

4.1 首先进行KEGG注释

diff.kk <- enrichKEGG(gene         = diff_geneList,
                      organism     = 'ath',
                      pvalueCutoff = 0.99,
                      qvalueCutoff=0.99
)
kegg_diff_dt <- as.data.frame(setReadable(diff.kk,org.At.tair.db,keytype = 'TAIR'))

up.kk <- enrichKEGG(gene         = up_geneList,
                    organism     = 'ath',
                    pvalueCutoff = 0.99,
                    qvalueCutoff=0.99
)
kegg_up_dt <- as.data.frame(setReadable(up.kk,org.At.tair.db,keytype = 'TAIR'))

down.kk <- enrichKEGG(gene         = down_geneList,
                      organism     = 'ath',
                      pvalueCutoff = 0.99,
                      qvalueCutoff=0.99
)
kegg_down_dt <- as.data.frame(setReadable(down.kk,org.At.tair.db,keytype = 'TAIR'))

4.2 可视化看看KEGG注释结果

## KEGG patheay visulization: 

  down_kegg<-kegg_down_dt[kegg_down_dt$pvalue<0.05,];down_kegg$group=-1
  up_kegg<-kegg_up_dt[kegg_up_dt$pvalue<0.05,];up_kegg$group=1
  dat=rbind(up_kegg,down_kegg)
  colnames(dat)

##  [1] "ID"          "Description" "GeneRatio"   "BgRatio"     "pvalue"     
##  [6] "p.adjust"    "qvalue"      "geneID"      "Count"       "group"

  dat$pvalue = -log10(dat$pvalue)
  dat$pvalue=dat$pvalue*dat$group 

  dat=dat[order(dat$pvalue,decreasing = F),]

  g_kegg<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + 
    geom_bar(stat="identity") + 
    scale_fill_gradient(low="blue",high="red",guide = FALSE) + 
    scale_x_discrete(name ="Pathway names") +
    scale_y_continuous(name ="-log10P-value") +
    coord_flip() +
    ggtitle("Pathway Enrichment")
  print(g_kegg)

4.3 接着进行GO注释

for (onto in c('CC','BP','MF')){

  ego <- enrichGO(gene         = diff_geneList,
                  OrgDb         = org.At.tair.db, 
                  keyType = 'TAIR',
                  ont           =  onto ,
                  pAdjustMethod = "BH",
                  pvalueCutoff  = 0.2,
                  qvalueCutoff  = 0.9)
  ego <- setReadable(ego, org.At.tair.db,keytype = 'TAIR')
  write.csv(as.data.frame(ego),paste0(prefix,"_",onto,".csv"))
  #ego2 <- gofilter(ego,4)
  ego2=ego
  pdf(paste0(prefix,"_",onto,'_barplot.pdf'))
  p=barplot(ego2, showCategory=12)+scale_x_discrete(labels=function(x) str_wrap(x,width=20))
  print(p)
  dev.off() 

}
ggsave(filename = paste0(prefix,"_volcano_plot.pdf"),g_volcano)

## Saving 7 x 5 in image

ggsave(filename = paste0(prefix,"_kegg_plot.pdf"),g_kegg)

## Saving 7 x 5 in image

write.csv(x = kegg_diff_dt,file = paste0(prefix,"_kegg_diff.csv"))
write.csv(x = kegg_up_dt,file = paste0(prefix,"_kegg_up.csv"))
write.csv(x = kegg_down_dt,file = paste0(prefix,"_kegg_down.csv"))

5 基因ID注释

library(org.At.tair.db)
ls('package:org.At.tair.db')

##  [1] "org.At.tair"             "org.At.tair.db"         
##  [3] "org.At.tairARACYC"       "org.At.tairARACYCENZYME"
##  [5] "org.At.tairCHR"          "org.At.tairCHRLENGTHS"  
##  [7] "org.At.tairCHRLOC"       "org.At.tairCHRLOCEND"   
##  [9] "org.At.tairENTREZID"     "org.At.tairENZYME"      
## [11] "org.At.tairENZYME2TAIR"  "org.At.tairGENENAME"    
## [13] "org.At.tairGO"           "org.At.tairGO2ALLTAIRS" 
## [15] "org.At.tairGO2TAIR"      "org.At.tairMAPCOUNTS"   
## [17] "org.At.tairORGANISM"     "org.At.tairPATH"        
## [19] "org.At.tairPATH2TAIR"    "org.At.tairPMID"        
## [21] "org.At.tairPMID2TAIR"    "org.At.tairREFSEQ"      
## [23] "org.At.tairREFSEQ2TAIR"  "org.At.tairSYMBOL"      
## [25] "org.At.tair_dbInfo"      "org.At.tair_dbconn"     
## [27] "org.At.tair_dbfile"      "org.At.tair_dbschema"

## Then draw GO/kegg figures:
deg1$gene_id=rownames(deg1)
id2symbol=toTable(org.At.tairSYMBOL) 
deg1=merge(deg1,id2symbol,by='gene_id')
## 可以看到有一些ID是没有gene symbol的，这样的基因，意义可能不大，也有可能是大部分人漏掉了
head(deg1)

##     gene_id   baseMean log2FoldChange     lfcSE       stat    pvalue
## 1 AT1G01010   58.25657     1.13105335 0.8000274  1.4137683 0.1574300
## 2 AT1G01010   58.25657     1.13105335 0.8000274  1.4137683 0.1574300
## 3 AT1G01020  542.64394    -0.05745554 0.3721650 -0.1543819 0.8773086
## 4 AT1G01030   48.77442    -1.09845263 1.2875862 -0.8531100 0.3935983
## 5 AT1G01040 1708.68949     0.32421734 0.2777530  1.1672865 0.2430947
## 6 AT1G01040 1708.68949     0.32421734 0.2777530  1.1672865 0.2430947
##        padj change  symbol
## 1 0.6008903    NOT ANAC001
## 2 0.6008903    NOT  NAC001
## 3 0.9805661    NOT    ARV1
## 4 0.8144858    NOT    NGA3
## 5 0.6992279    NOT    DCL1
## 6 0.6992279    NOT   EMB60