-
Notifications
You must be signed in to change notification settings - Fork 0
/
R-code
56 lines (46 loc) · 2.52 KB
/
R-code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
############# install the required package for searching the annotation files
BiocManager::install("AnnotationHub")
BiocManager::install("biomaRt")
library(AnnotationHub)
library(biomaRt)
hub <- AnnotationHub::AnnotationHub()
query(hub,"hypogaea") # hypogaea means the name of species you want to annotate
hypogaea.OrgDb <- hub[["AH85392"]] # after query,get the package's name "AH85392"
#############
############# install package for GO and KEGG
BiocManager::install("clusterProfiler")
BiocManager::install("topGO")
BiocManager::install("Rgraphviz")
BiocManager::install("pathview")
library(clusterProfiler)
library(topGO)
library(Rgraphviz)
############ this step is to show what kind of keytype(the Gene name type) the package have.
columns(hypogaea.OrgDb)
############
############ This example clearly show what every keytype looks like.cols=XXX, XXX is the keytype you want to match
select(hypogaea.OrgDb,keys = "112768463",columns = cols,keytype = "ENTREZID")
ENTREZID ACCNUM ALIAS CHR
112768463 XR_003185888.1 LOC112768463 Arahy.17
GENENAME GID PMID REFSEQ SYMBOL
28S ribosomal RNA 112768463 <NA> XR_003185888.1 LOC112768463
########### The keytype of different speices may be different, you can use below command to know the detail of every keytype of your downloaded database
head(keys(hypogaea.OrgDb,keytype="XXX"))
############ GO enrichment
enrich <- enrichGO(gene = DEG.peanut,OrgDb = hypogaea.OrgDb,keyType = "ACCNUM",ont = "MF",readable = T)
dotplot(enrich)
# sometimes it shows GOALL error, it is because the package don't have GO information. So we can't do GO analysis
############
############ ID transform.
KEGGID <- mapIds(x=hypogaea.OrgDb,keys = KEGGIDA,keytype = "ACCNUM",column = "ENTREZID")
############ KEGG analysis
enKEGG <- enrichKEGG(gene = KEGGIDA,organism = "ahf",pvalueCutoff = 0.05,qvalueCutoff = 0.05)
# organism is the species' name in KEGG. Ahf means human. You can search from "https://www.genome.jp/kegg/catalog/org_list.html"
dotplot(enKEGG,title='KEGG pathway of DEGs', showCategory=20)
# organism: every specie has different organism number , you can get it by searching at KEGG website. This command don't need local database. It works online.
############ subest can filter data for specify condition
MD_P <- subset(MD,MD$GMM_logit_pvalue<0.05)
################# Every time you want to do GO or KEGG, you should rerun below document to load the data even you have done some step below.
hub <- AnnotationHub::AnnotationHub()
query(hub,"hypogaea")
hypogaea.OrgDb <- hub[["AH85392"]]