拟南芥的差异分析结果注释

11/9/2017

  • 1 首先加载必要的包

  • 2 然后导入数据

  • 3 然后判断显著差异基因

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

  • 4 GO/KEGG注释

    • 4.1 首先进行KEGG注释

    • 4.2 可视化看看KEGG注释结果

    • 4.3 接着进行GO注释

  • 5 基因ID注释

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-09prefix='Day1-Day0'

3 然后判断显著差异基因

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

table(deg1$padj<0.05)## ## FALSE  TRUE ## 19166   197table(deg1$padj<0.01)## ## FALSE  TRUE ## 19303    60diff_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] 197length(up_geneList)## [1] 89length(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 imageggsave(filename = paste0(prefix,"_kegg_plot.pdf"),g_kegg)## Saving 7 x 5 in imagewrite.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

可以看到基因ID被注释到了基因名，勉强可以认识了。

赞 (0)

相关推荐