GO/KEGG Pathway analysis
Kelly Yan Ren
03 September, 2020
true
1 Major analysis
1.1 Loading packages
library(readr)
library(ggplot2)
library(dplyr)
library(pander)
library(ggrepel)
library(limma)
library(tibble)
library(AnnotationHub)
library(magrittr)
library(biomaRt)
library(GOfuncR)
library(tidyr)1.2 Gene Data import
setwd("~/Documents/Kelly_2020/Research/Six_M_Work_2020/Jo_RNAseq/Kelly_JoRNAseq_DEG_Parthway_R_codes")ah <- read_rds("All_annotation.rds")
#mcols(ah)
#subset(ah, rdataclass == "EnsDb" & species == "Bos taurus")
ensDb <- ah[["AH74957"]]
ensDb## EnsDb for Ensembl:
## |Backend: SQLite
## |Db type: EnsDb
## |Type of Gene ID: Ensembl Gene ID
## |Supporting package: ensembldb
## |Db created by: ensembldb package from Bioconductor
## |script_version: 0.3.5
## |Creation time: Sun Nov 17 21:16:03 2019
## |ensembl_version: 98
## |ensembl_host: localhost
## |Organism: Bos taurus
## |taxonomy_id: 9913
## |genome_build: ARS-UCD1.2
## |DBSCHEMAVERSION: 2.1
## | No. of genes: 27607.
## | No. of transcripts: 43984.
## |Protein data available.genesGR <- genes(ensDb)
genesGR## GRanges object with 27607 ranges and 8 metadata columns:
## seqnames ranges strand | gene_id
## <Rle> <IRanges> <Rle> | <character>
## ENSBTAG00000006648 1 339070-350389 - | ENSBTAG00000006648
## ENSBTAG00000049697 1 475398-475516 + | ENSBTAG00000049697
## ENSBTAG00000047028 1 477378-477504 + | ENSBTAG00000047028
## ENSBTAG00000053686 1 479417-479535 + | ENSBTAG00000053686
## ENSBTAG00000054829 1 626079-668340 + | ENSBTAG00000054829
## ... ... ... ... . ...
## ENSBTAG00000013341 X 138696888-138714389 - | ENSBTAG00000013341
## ENSBTAG00000046328 X 138722235-138751873 - | ENSBTAG00000046328
## ENSBTAG00000008114 X 138754709-138785187 - | ENSBTAG00000008114
## ENSBTAG00000049411 X 138855866-139001761 + | ENSBTAG00000049411
## ENSBTAG00000044391 X 138862220-138862348 - | ENSBTAG00000044391
## gene_name gene_biotype seq_coord_system
## <character> <character> <character>
## ENSBTAG00000006648 protein_coding primary_assembly
## ENSBTAG00000049697 5S_rRNA rRNA primary_assembly
## ENSBTAG00000047028 5S_rRNA rRNA primary_assembly
## ENSBTAG00000053686 5S_rRNA rRNA primary_assembly
## ENSBTAG00000054829 protein_coding primary_assembly
## ... ... ... ...
## ENSBTAG00000013341 GYG2 protein_coding primary_assembly
## ENSBTAG00000046328 XG protein_coding primary_assembly
## ENSBTAG00000008114 CD99 protein_coding primary_assembly
## ENSBTAG00000049411 protein_coding primary_assembly
## ENSBTAG00000044391 snoRNA primary_assembly
## description
## <character>
## ENSBTAG00000006648 NULL
## ENSBTAG00000049697 5S ribosomal RNA [Source:RFAM;Acc:RF00001]
## ENSBTAG00000047028 5S ribosomal RNA [Source:RFAM;Acc:RF00001]
## ENSBTAG00000053686 5S ribosomal RNA [Source:RFAM;Acc:RF00001]
## ENSBTAG00000054829 NULL
## ... ...
## ENSBTAG00000013341 glycogenin 2 [Source:VGNC Symbol;Acc:VGNC:29729]
## ENSBTAG00000046328 Xg glycoprotein [Source:VGNC Symbol;Acc:VGNC:56960]
## ENSBTAG00000008114 CD99 molecule [Source:VGNC Symbol;Acc:VGNC:50554]
## ENSBTAG00000049411 dehydrogenase/reductase X-linked [Source:NCBI gene;Acc:513482]
## ENSBTAG00000044391 NULL
## gene_id_version symbol entrezid
## <character> <character> <list>
## ENSBTAG00000006648 ENSBTAG00000006648.6 NA
## ENSBTAG00000049697 ENSBTAG00000049697.1 5S_rRNA NA
## ENSBTAG00000047028 ENSBTAG00000047028.2 5S_rRNA NA
## ENSBTAG00000053686 ENSBTAG00000053686.1 5S_rRNA NA
## ENSBTAG00000054829 ENSBTAG00000054829.1 NA
## ... ... ... ...
## ENSBTAG00000013341 ENSBTAG00000013341.6 GYG2 505258
## ENSBTAG00000046328 ENSBTAG00000046328.2 XG 101902671
## ENSBTAG00000008114 ENSBTAG00000008114.6 CD99 509230
## ENSBTAG00000049411 ENSBTAG00000049411.1 513482
## ENSBTAG00000044391 ENSBTAG00000044391.2 NA
## -------
## seqinfo: 178 sequences from ARS-UCD1.2 genomecols2Keep <- c("gene_id","gene_name", "gene_biotype", "description", "entrezid")
mcols(genesGR) <- mcols(genesGR)[, cols2Keep]
mcols(genesGR)## DataFrame with 27607 rows and 5 columns
## gene_id gene_name gene_biotype
## <character> <character> <character>
## ENSBTAG00000006648 ENSBTAG00000006648 protein_coding
## ENSBTAG00000049697 ENSBTAG00000049697 5S_rRNA rRNA
## ENSBTAG00000047028 ENSBTAG00000047028 5S_rRNA rRNA
## ENSBTAG00000053686 ENSBTAG00000053686 5S_rRNA rRNA
## ENSBTAG00000054829 ENSBTAG00000054829 protein_coding
## ... ... ... ...
## ENSBTAG00000013341 ENSBTAG00000013341 GYG2 protein_coding
## ENSBTAG00000046328 ENSBTAG00000046328 XG protein_coding
## ENSBTAG00000008114 ENSBTAG00000008114 CD99 protein_coding
## ENSBTAG00000049411 ENSBTAG00000049411 protein_coding
## ENSBTAG00000044391 ENSBTAG00000044391 snoRNA
## description
## <character>
## ENSBTAG00000006648 NULL
## ENSBTAG00000049697 5S ribosomal RNA [Source:RFAM;Acc:RF00001]
## ENSBTAG00000047028 5S ribosomal RNA [Source:RFAM;Acc:RF00001]
## ENSBTAG00000053686 5S ribosomal RNA [Source:RFAM;Acc:RF00001]
## ENSBTAG00000054829 NULL
## ... ...
## ENSBTAG00000013341 glycogenin 2 [Source:VGNC Symbol;Acc:VGNC:29729]
## ENSBTAG00000046328 Xg glycoprotein [Source:VGNC Symbol;Acc:VGNC:56960]
## ENSBTAG00000008114 CD99 molecule [Source:VGNC Symbol;Acc:VGNC:50554]
## ENSBTAG00000049411 dehydrogenase/reductase X-linked [Source:NCBI gene;Acc:513482]
## ENSBTAG00000044391 NULL
## entrezid
## <list>
## ENSBTAG00000006648 NA
## ENSBTAG00000049697 NA
## ENSBTAG00000047028 NA
## ENSBTAG00000053686 NA
## ENSBTAG00000054829 NA
## ... ...
## ENSBTAG00000013341 505258
## ENSBTAG00000046328 101902671
## ENSBTAG00000008114 509230
## ENSBTAG00000049411 513482
## ENSBTAG00000044391 NAALL_entrezID <- genesGR %>%
subset(!is.na(entrezid)) %>%
mcols() %>%
.[["entrezid"]] %>%
unlist() %>%
unique()
length(ALL_entrezID) #18699## [1] 18699#Bos_taurus = useMart(biomart="ensembl", dataset="btaurus_gene_ensembl") #ARS-UCD1.2
### GO terms match their entrezgene_id and ensembl_gene_id
#ALL_GO2gene <- getBM(
# attributes= c("chromosome_name","start_position","end_position","strand","ensembl_gene_id","external_gene_name","gene_biotype", "description","entrezgene_id", "go_id"),mart= Bos_taurus)
# saveRDS(ALL_GO2gene, "ALL_GO2gene.rds")
ALL_GO2gene <- read_rds("ALL_GO2gene.rds")
diff_GO2gene <- subset(ALL_GO2gene, entrezgene_id %in% setdiff(ALL_GO2gene$entrezgene_id, ALL_entrezID))%>%drop_na()
diff_GO2gene$strand <- diff_GO2gene$strand%>%gsub("-1", "-", .)%>%gsub("1", "+", .)
diff_GO2gene$width <- diff_GO2gene[,"end_position"]-diff_GO2gene[,"start_position"]
diff_GO2gene <- diff_GO2gene[,c(1:3)]%>%
cbind(diff_GO2gene$width)%>%
cbind(diff_GO2gene[,c(4:9)])%>%
set_colnames(colnames(data.frame(genesGR)))data_genesGR <- as.data.frame(genesGR)%>%
rbind(diff_GO2gene)
ALL_entrezID <- data_genesGR %>%
subset(!is.na(entrezid)) %>%
.[["entrezid"]] %>%
unlist() %>%
unique()
length(ALL_entrezID) #18730 combine two datasets## [1] 187301.3 phenotype
1.3.0.1 frontal_skin
polledvshornedFS_DEG <- read_csv("polledvshornedForfrontal_skin-ARS_jo-DEG_all.csv")polledvshornedFS_DEG[!((polledvshornedFS_DEG$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
polledvshornedFS_entrezID_DEG <- polledvshornedFS_DEG%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(polledvshornedFS_entrezID_DEG) ## [1] 52641.3.1 Go entrichment analysis
goRes <- goana(polledvshornedFS_entrezID_DEG, ALL_entrezID, species = "Bt")
goRes <- goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
goRes%>%print(n= dim(goRes)[1])## # A tibble: 16 x 7
## goid Term Ont N DE P.DE adjP
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 GO:00055… extracellular region CC 546 225 1.85e-11 2.62e-7
## 2 GO:00076… sensory perception BP 114 58 2.27e- 7 1.01e-3
## 3 GO:00076… visual perception BP 63 37 3.56e- 7 1.01e-3
## 4 GO:00509… sensory perception of light s… BP 63 37 3.56e- 7 1.01e-3
## 5 GO:00056… extracellular space CC 243 105 2.97e- 7 1.01e-3
## 6 GO:00071… G protein-coupled receptor si… BP 149 70 7.12e- 7 1.68e-3
## 7 GO:00508… nervous system process BP 169 76 2.03e- 6 4.10e-3
## 8 GO:00030… system process BP 311 123 7.81e- 6 1.23e-2
## 9 GO:00058… intermediate filament CC 45 27 7.74e- 6 1.23e-2
## 10 GO:00108… regulation of hormone levels BP 93 46 1.03e- 5 1.46e-2
## 11 GO:00049… G protein-coupled receptor ac… MF 68 36 1.38e- 5 1.77e-2
## 12 GO:00095… fertilization BP 50 28 3.19e- 5 3.39e-2
## 13 GO:00380… signaling receptor activity MF 146 64 3.35e- 5 3.39e-2
## 14 GO:00600… molecular transducer activity MF 146 64 3.35e- 5 3.39e-2
## 15 GO:00072… neuropeptide signaling pathway BP 25 17 3.94e- 5 3.72e-2
## 16 GO:00069… acute-phase response BP 15 12 4.50e- 5 3.98e-2write.csv(goRes%>%as.data.frame(),file="polledvshornedForfrontal_skin-ARS_jo-GOresults_p0.05_lfc1.csv")1.3.2 KEGG entrichment analysis
keggRes <- kegga(polledvshornedFS_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
keggRes <- keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
keggRes%>%print(n= dim(keggRes)[1])## # A tibble: 72 x 6
## keggid Pathway N DE P.DE adjP
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 path:bta0… Neuroactive ligand-receptor interac… 345 188 2.76e-25 9.19e-23
## 2 path:bta0… Staphylococcus aureus infection 101 71 1.44e-18 2.40e-16
## 3 path:bta0… Steroid hormone biosynthesis 70 52 1.33e-15 1.48e-13
## 4 path:bta0… Cytokine-cytokine receptor interact… 296 145 1.51e-14 1.26e-12
## 5 path:bta0… Linoleic acid metabolism 33 28 1.72e-11 1.14e- 9
## 6 path:bta0… Retinol metabolism 65 44 3.71e-11 2.06e- 9
## 7 path:bta0… Fat digestion and absorption 49 36 5.59e-11 2.66e- 9
## 8 path:bta0… Protein digestion and absorption 113 63 5.97e-10 2.48e- 8
## 9 path:bta0… Arachidonic acid metabolism 78 48 7.19e-10 2.66e- 8
## 10 path:bta0… Pancreatic secretion 100 57 1.25e- 9 4.17e- 8
## 11 path:bta0… Natural killer cell mediated cytoto… 133 70 2.01e- 9 6.09e- 8
## 12 path:bta0… Complement and coagulation cascades 83 49 3.58e- 9 9.68e- 8
## 13 path:bta0… Graft-versus-host disease 56 37 3.78e- 9 9.68e- 8
## 14 path:bta0… Hematopoietic cell lineage 103 57 5.61e- 9 1.33e- 7
## 15 path:bta0… Type I diabetes mellitus 54 34 9.70e- 8 2.15e- 6
## 16 path:bta0… Allograft rejection 50 32 1.32e- 7 2.75e- 6
## 17 path:bta0… Drug metabolism - cytochrome P450 60 36 2.38e- 7 4.67e- 6
## 18 path:bta0… Bile secretion 94 50 2.54e- 7 4.70e- 6
## 19 path:bta0… Chemical carcinogenesis 77 43 2.90e- 7 5.08e- 6
## 20 path:bta0… Salivary secretion 91 48 6.09e- 7 1.01e- 5
## 21 path:bta0… Autoimmune thyroid disease 60 35 9.02e- 7 1.43e- 5
## 22 path:bta0… ABC transporters 56 33 1.34e- 6 2.03e- 5
## 23 path:bta0… Amoebiasis 107 52 5.31e- 6 7.69e- 5
## 24 path:bta0… cAMP signaling pathway 219 92 6.29e- 6 8.06e- 5
## 25 path:bta0… JAK-STAT signaling pathway 179 78 6.28e- 6 8.06e- 5
## 26 path:bta0… Intestinal immune network for IgA p… 54 31 5.98e- 6 8.06e- 5
## 27 path:bta0… Tyrosine metabolism 38 24 6.98e- 6 8.31e- 5
## 28 path:bta0… Inflammatory bowel disease 69 37 6.98e- 6 8.31e- 5
## 29 path:bta0… Cell adhesion molecules 152 68 8.27e- 6 9.18e- 5
## 30 path:bta0… Nicotine addiction 36 23 8.20e- 6 9.18e- 5
## 31 path:bta0… Asthma 34 22 9.55e- 6 1.03e- 4
## 32 path:bta0… alpha-Linolenic acid metabolism 26 18 1.57e- 5 1.58e- 4
## 33 path:bta0… Maturity onset diabetes of the young 26 18 1.57e- 5 1.58e- 4
## 34 path:bta0… Metabolism of xenobiotics by cytoch… 66 35 1.69e- 5 1.66e- 4
## 35 path:bta0… Carbohydrate digestion and absorpti… 42 25 2.01e- 5 1.91e- 4
## 36 path:bta0… Vascular smooth muscle contraction 136 61 2.14e- 5 1.98e- 4
## 37 path:bta0… Antigen processing and presentation 80 40 2.72e- 5 2.45e- 4
## 38 path:bta0… Viral protein interaction with cyto… 91 43 7.76e- 5 6.80e- 4
## 39 path:bta0… Ovarian steroidogenesis 57 30 8.02e- 5 6.85e- 4
## 40 path:bta0… Chagas disease 114 51 1.06e- 4 8.81e- 4
## 41 path:bta0… Rheumatoid arthritis 98 45 1.25e- 4 1.02e- 3
## 42 path:bta0… Viral myocarditis 69 34 1.54e- 4 1.22e- 3
## 43 path:bta0… Th17 cell differentiation 110 49 1.63e- 4 1.26e- 3
## 44 path:bta0… Prolactin signaling pathway 83 39 1.90e- 4 1.41e- 3
## 45 path:bta0… Primary immunodeficiency 39 22 1.91e- 4 1.41e- 3
## 46 path:bta0… Leishmaniasis 75 36 1.95e- 4 1.41e- 3
## 47 path:bta0… Ether lipid metabolism 50 26 3.04e- 4 2.13e- 3
## 48 path:bta0… Tuberculosis 188 75 3.07e- 4 2.13e- 3
## 49 path:bta0… Inflammatory mediator regulation of… 97 43 4.52e- 4 3.07e- 3
## 50 path:bta0… Phagosome 161 65 5.10e- 4 3.40e- 3
## 51 path:bta0… African trypanosomiasis 39 21 6.22e- 4 4.06e- 3
## 52 path:bta0… Mineral absorption 55 27 7.66e- 4 4.91e- 3
## 53 path:bta0… Estrogen signaling pathway 129 53 9.97e- 4 6.26e- 3
## 54 path:bta0… Phototransduction 26 15 1.48e- 3 9.10e- 3
## 55 path:bta0… Phenylalanine metabolism 24 14 1.83e- 3 1.11e- 2
## 56 path:bta0… Insulin secretion 83 36 2.02e- 3 1.20e- 2
## 57 path:bta0… Drug metabolism - other enzymes 75 33 2.30e- 3 1.32e- 2
## 58 path:bta0… Th1 and Th2 cell differentiation 95 40 2.31e- 3 1.32e- 2
## 59 path:bta0… T cell receptor signaling pathway 105 43 3.05e- 3 1.72e- 2
## 60 path:bta0… Nitrogen metabolism 18 11 3.47e- 3 1.92e- 2
## 61 path:bta0… Ascorbate and aldarate metabolism 23 13 3.90e- 3 2.13e- 2
## 62 path:bta0… Folate biosynthesis 36 18 4.38e- 3 2.32e- 2
## 63 path:bta0… Type II diabetes mellitus 44 21 4.38e- 3 2.32e- 2
## 64 path:bta0… Histidine metabolism 21 12 4.90e- 3 2.55e- 2
## 65 path:bta0… Yersinia infection 141 54 5.47e- 3 2.80e- 2
## 66 path:bta0… Serotonergic synapse 111 44 5.62e- 3 2.84e- 2
## 67 path:bta0… Calcium signaling pathway 195 71 6.91e- 3 3.43e- 2
## 68 path:bta0… Regulation of lipolysis in adipocyt… 54 24 7.42e- 3 3.58e- 2
## 69 path:bta0… Pertussis 74 31 7.40e- 3 3.58e- 2
## 70 path:bta0… Pentose and glucuronate interconver… 30 15 9.01e- 3 4.27e- 2
## 71 path:bta0… Thyroid hormone synthesis 72 30 9.11e- 3 4.27e- 2
## 72 path:bta0… Renin-angiotensin system 25 13 9.84e- 3 4.55e- 2write.csv(keggRes%>%as.data.frame(),file="polledvshornedForfrontal_skin-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.3.2.1 hornbud
polledvshornedHB_DEG <- read_csv("polledvshornedForhornbud-ARS_jo-DEG_all.csv")polledvshornedHB_DEG[!((polledvshornedHB_DEG$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
polledvshornedHB_entrezID_DEG <- polledvshornedHB_DEG%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(polledvshornedHB_entrezID_DEG) ## [1] 56001.3.3 Go entrichment analysis
goRes <- goana(polledvshornedHB_entrezID_DEG, ALL_entrezID, species = "Bt")
goRes <- goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
goRes%>%print(n= dim(goRes)[1])## # A tibble: 20 x 7
## goid Term Ont N DE P.DE adjP
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 GO:0005… extracellular region CC 546 240 1.38e-12 1.96e-8
## 2 GO:0007… visual perception BP 63 44 7.02e-11 3.31e-7
## 3 GO:0050… sensory perception of light st… BP 63 44 7.02e-11 3.31e-7
## 4 GO:0007… sensory perception BP 114 65 1.53e- 9 5.43e-6
## 5 GO:0005… intermediate filament CC 45 33 2.30e- 9 6.50e-6
## 6 GO:0005… extracellular space CC 243 114 1.52e- 8 3.57e-5
## 7 GO:0050… nervous system process BP 169 84 4.90e- 8 9.90e-5
## 8 GO:0045… intermediate filament cytoskel… CC 49 33 6.88e- 8 1.22e-4
## 9 GO:0003… system process BP 311 134 4.85e- 7 7.62e-4
## 10 GO:0019… cytolysis BP 17 15 9.60e- 7 1.36e-3
## 11 GO:0003… lysozyme activity MF 12 11 1.49e- 5 1.75e-2
## 12 GO:0061… peptidoglycan muralytic activi… MF 12 11 1.49e- 5 1.75e-2
## 13 GO:0005… collagen trimer CC 25 18 1.71e- 5 1.87e-2
## 14 GO:0044… multi-organism reproductive pr… BP 207 90 2.26e- 5 2.28e-2
## 15 GO:0019… sexual reproduction BP 183 81 2.63e- 5 2.48e-2
## 16 GO:0045… intermediate filament-based pr… BP 16 13 3.23e- 5 2.69e-2
## 17 GO:0045… intermediate filament cytoskel… BP 16 13 3.23e- 5 2.69e-2
## 18 GO:0051… multi-organism process BP 210 90 4.37e- 5 3.25e-2
## 19 GO:0015… sodium ion transmembrane trans… MF 22 16 4.22e- 5 3.25e-2
## 20 GO:0007… G protein-coupled receptor sig… BP 149 67 6.98e- 5 4.94e-2write.csv(goRes%>%as.data.frame(),file="polledvshornedForhornbud-ARS_jo-GOresults_p0.05_lfc1.csv")1.3.4 KEGG entrichment analysis
keggRes <- kegga(polledvshornedHB_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
keggRes <- keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
keggRes%>%print(n= dim(keggRes)[1])## # A tibble: 51 x 6
## keggid Pathway N DE P.DE adjP
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 path:bta0… Neuroactive ligand-receptor interac… 345 192 9.67e-24 3.22e-21
## 2 path:bta0… Staphylococcus aureus infection 101 76 5.62e-21 9.35e-19
## 3 path:bta0… Steroid hormone biosynthesis 70 52 2.17e-14 1.93e-12
## 4 path:bta0… Protein digestion and absorption 113 73 2.32e-14 1.93e-12
## 5 path:bta0… Cytokine-cytokine receptor interact… 296 150 3.96e-14 2.64e-12
## 6 path:bta0… Linoleic acid metabolism 33 28 8.65e-11 4.80e- 9
## 7 path:bta0… Arachidonic acid metabolism 78 49 1.74e- 9 7.25e- 8
## 8 path:bta0… Retinol metabolism 65 43 1.64e- 9 7.25e- 8
## 9 path:bta0… Fat digestion and absorption 49 35 2.36e- 9 8.72e- 8
## 10 path:bta0… Graft-versus-host disease 56 38 4.94e- 9 1.64e- 7
## 11 path:bta0… Complement and coagulation cascades 83 49 3.29e- 8 9.95e- 7
## 12 path:bta0… Chemical carcinogenesis 77 46 5.19e- 8 1.44e- 6
## 13 path:bta0… Natural killer cell mediated cytoto… 133 69 8.80e- 8 2.25e- 6
## 14 path:bta0… Drug metabolism - cytochrome P450 60 37 3.36e- 7 7.99e- 6
## 15 path:bta0… Type I diabetes mellitus 54 33 1.93e- 6 4.29e- 5
## 16 path:bta0… Pancreatic secretion 100 52 3.00e- 6 6.24e- 5
## 17 path:bta0… Amoebiasis 107 54 6.30e- 6 1.23e- 4
## 18 path:bta0… Allograft rejection 50 30 9.56e- 6 1.77e- 4
## 19 path:bta0… Autoimmune thyroid disease 60 34 1.42e- 5 2.48e- 4
## 20 path:bta0… ABC transporters 56 32 2.00e- 5 3.17e- 4
## 21 path:bta0… Antigen processing and presentation 80 42 1.92e- 5 3.17e- 4
## 22 path:bta0… Metabolism of xenobiotics by cytoch… 66 36 2.51e- 5 3.49e- 4
## 23 path:bta0… Phagosome 161 73 2.35e- 5 3.49e- 4
## 24 path:bta0… Nicotine addiction 36 23 2.51e- 5 3.49e- 4
## 25 path:bta0… Viral protein interaction with cyto… 91 46 2.82e- 5 3.62e- 4
## 26 path:bta0… Salivary secretion 91 46 2.82e- 5 3.62e- 4
## 27 path:bta0… Bile secretion 94 47 3.33e- 5 4.11e- 4
## 28 path:bta0… alpha-Linolenic acid metabolism 26 18 3.97e- 5 4.56e- 4
## 29 path:bta0… Maturity onset diabetes of the young 26 18 3.97e- 5 4.56e- 4
## 30 path:bta0… Hematopoietic cell lineage 103 50 5.16e- 5 5.73e- 4
## 31 path:bta0… JAK-STAT signaling pathway 179 78 6.83e- 5 7.34e- 4
## 32 path:bta0… Estrogen signaling pathway 129 59 1.01e- 4 1.05e- 3
## 33 path:bta0… Primary immunodeficiency 39 23 1.51e- 4 1.53e- 3
## 34 path:bta0… Phototransduction 26 17 1.94e- 4 1.90e- 3
## 35 path:bta0… Viral myocarditis 69 35 2.26e- 4 2.15e- 3
## 36 path:bta0… Tyrosine metabolism 38 22 3.04e- 4 2.81e- 3
## 37 path:bta0… Ether lipid metabolism 50 27 3.15e- 4 2.84e- 3
## 38 path:bta0… Calcium signaling pathway 195 81 3.39e- 4 2.97e- 3
## 39 path:bta0… Asthma 34 20 4.27e- 4 3.64e- 3
## 40 path:bta0… Intestinal immune network for IgA p… 54 28 5.92e- 4 4.93e- 3
## 41 path:bta0… Cell adhesion molecules 152 64 8.88e- 4 7.21e- 3
## 42 path:bta0… Vascular smooth muscle contraction 136 58 1.04e- 3 8.27e- 3
## 43 path:bta0… cAMP signaling pathway 219 87 1.14e- 3 8.83e- 3
## 44 path:bta0… Inflammatory bowel disease 69 33 1.27e- 3 9.65e- 3
## 45 path:bta0… Drug metabolism - other enzymes 75 35 1.59e- 3 1.17e- 2
## 46 path:bta0… Carbohydrate digestion and absorpti… 42 22 1.88e- 3 1.36e- 2
## 47 path:bta0… Ascorbate and aldarate metabolism 23 14 2.00e- 3 1.42e- 2
## 48 path:bta0… Ovarian steroidogenesis 57 27 4.02e- 3 2.79e- 2
## 49 path:bta0… Malaria 55 26 4.88e- 3 3.32e- 2
## 50 path:bta0… Renin-angiotensin system 25 14 5.77e- 3 3.84e- 2
## 51 path:bta0… Pentose and glucuronate interconver… 30 16 6.11e- 3 3.99e- 2write.csv(keggRes%>%as.data.frame(),file="polledvshornedForhornbud-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.3.4.1 forebrain
polledvshornedFB_DEG <- read_csv("polledvshornedForforebrain-ARS_jo-DEG_all.csv")polledvshornedFB_DEG[!((polledvshornedFB_DEG$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
polledvshornedFB_entrezID_DEG <- polledvshornedFB_DEG%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(polledvshornedFB_entrezID_DEG) ## [1] 34721.3.5 Go entrichment analysis
goRes <- goana(polledvshornedFB_entrezID_DEG, ALL_entrezID, species = "Bt")
goRes <- goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
goRes%>%print(n= dim(goRes)[1])## # A tibble: 2 x 7
## goid Term Ont N DE P.DE adjP
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 GO:0005576 extracellular region CC 546 158 0.00000000117 0.0000165
## 2 GO:0005615 extracellular space CC 243 74 0.00000408 0.0288write.csv(goRes%>%as.data.frame(),file="polledvshornedForforebrain-ARS_jo-GOresults_p0.05_lfc1.csv")1.3.6 KEGG entrichment analysis
keggRes <- kegga(polledvshornedFB_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
keggRes <- keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
keggRes%>%print(n= dim(keggRes)[1])## # A tibble: 27 x 6
## keggid Pathway N DE P.DE adjP
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 path:bta0… Systemic lupus erythematosus 150 66 5.55e-13 1.85e-10
## 2 path:bta0… Protein digestion and absorption 113 50 2.74e-10 4.56e- 8
## 3 path:bta0… Alcoholism 195 71 2.60e- 9 2.88e- 7
## 4 path:bta0… Linoleic acid metabolism 33 20 1.03e- 7 8.56e- 6
## 5 path:bta0… Arachidonic acid metabolism 78 34 2.95e- 7 1.96e- 5
## 6 path:bta0… Neuroactive ligand-receptor interac… 345 102 3.60e- 7 2.00e- 5
## 7 path:bta0… Cytokine-cytokine receptor interact… 296 87 3.29e- 6 1.57e- 4
## 8 path:bta0… Viral protein interaction with cyto… 91 35 6.60e- 6 2.44e- 4
## 9 path:bta0… Vascular smooth muscle contraction 136 47 6.33e- 6 2.44e- 4
## 10 path:bta0… Protein export 23 14 8.30e- 6 2.76e- 4
## 11 path:bta0… Viral carcinogenesis 213 64 2.95e- 5 8.94e- 4
## 12 path:bta0… Staphylococcus aureus infection 101 36 3.48e- 5 9.65e- 4
## 13 path:bta0… alpha-Linolenic acid metabolism 26 14 5.59e- 5 1.43e- 3
## 14 path:bta0… Spliceosome 128 42 7.65e- 5 1.82e- 3
## 15 path:bta0… Fat digestion and absorption 49 20 2.39e- 4 5.31e- 3
## 16 path:bta0… Amoebiasis 107 35 3.12e- 4 6.48e- 3
## 17 path:bta0… Pancreatic secretion 100 33 3.81e- 4 7.46e- 3
## 18 path:bta0… Graft-versus-host disease 56 21 6.61e- 4 1.22e- 2
## 19 path:bta0… Antigen processing and presentation 80 27 8.46e- 4 1.41e- 2
## 20 path:bta0… Thyroid hormone synthesis 72 25 8.12e- 4 1.41e- 2
## 21 path:bta0… Inflammatory mediator regulation of… 97 31 1.03e- 3 1.60e- 2
## 22 path:bta0… Type I diabetes mellitus 54 20 1.06e- 3 1.60e- 2
## 23 path:bta0… Phototransduction 26 12 1.17e- 3 1.69e- 2
## 24 path:bta0… Estrogen signaling pathway 129 38 1.69e- 3 2.25e- 2
## 25 path:bta0… Oxytocin signaling pathway 146 42 1.66e- 3 2.25e- 2
## 26 path:bta0… Ether lipid metabolism 50 18 2.66e- 3 3.41e- 2
## 27 path:bta0… Phagosome 161 44 3.84e- 3 4.74e- 2write.csv(keggRes%>%as.data.frame(),file="polledvshornedForforebrain-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.3.6.1 midbrain
polledvshornedMB_DEG <- read_csv("polledvshornedFormidbrain-ARS_jo-DEG_all.csv")polledvshornedMB_DEG[!((polledvshornedMB_DEG$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
polledvshornedMB_entrezID_DEG <- polledvshornedMB_DEG%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(polledvshornedMB_entrezID_DEG) ## [1] 26521.3.7 Go entrichment analysis
goRes <- goana(polledvshornedMB_entrezID_DEG, ALL_entrezID, species = "Bt")
goRes <- goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
goRes%>%print(n= dim(goRes)[1])## # A tibble: 17 x 7
## goid Term Ont N DE P.DE adjP
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 GO:0005… extracellular region CC 546 122 1.17e-7 0.00166
## 2 GO:0007… visual perception BP 63 25 5.54e-7 0.00261
## 3 GO:0050… sensory perception of light sti… BP 63 25 5.54e-7 0.00261
## 4 GO:0099… supramolecular complex CC 287 70 2.46e-6 0.00868
## 5 GO:0032… multicellular organismal process BP 1326 245 3.82e-6 0.0108
## 6 GO:0003… RNA binding MF 257 62 1.30e-5 0.0306
## 7 GO:0005… protein binding MF 1383 250 1.56e-5 0.0316
## 8 GO:0010… regulation of gene expression BP 755 147 2.64e-5 0.0381
## 9 GO:0060… regulation of macromolecule met… BP 1105 204 2.69e-5 0.0381
## 10 GO:0008… RNA splicing BP 113 33 2.63e-5 0.0381
## 11 GO:0031… extracellular matrix CC 59 21 3.15e-5 0.0405
## 12 GO:0019… regulation of metabolic process BP 1203 218 4.79e-5 0.0463
## 13 GO:0048… system development BP 788 151 4.76e-5 0.0463
## 14 GO:0051… regulation of nitrogen compound… BP 995 185 4.25e-5 0.0463
## 15 GO:0007… cytoskeleton organization BP 246 58 4.91e-5 0.0463
## 16 GO:2000… regulation of multicellular org… BP 319 71 5.65e-5 0.0470
## 17 GO:0005… collagen trimer CC 25 12 5.49e-5 0.0470write.csv(goRes%>%as.data.frame(),file="polledvshornedFormidbrain-ARS_jo-GOresults_p0.05_lfc1.csv")1.3.8 KEGG entrichment analysis
keggRes <- kegga(polledvshornedMB_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
keggRes <- keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
keggRes%>%print(n= dim(keggRes)[1])## # A tibble: 7 x 6
## keggid Pathway N DE P.DE adjP
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 path:bta030… Spliceosome 128 41 1.95e-7 3.25e-5
## 2 path:bta049… Protein digestion and absorption 113 38 1.31e-7 3.25e-5
## 3 path:bta040… Cytokine-cytokine receptor inter… 296 69 1.46e-5 1.62e-3
## 4 path:bta046… Antigen processing and presentat… 80 25 7.12e-5 4.74e-3
## 5 path:bta049… Fat digestion and absorption 49 18 7.02e-5 4.74e-3
## 6 path:bta041… Phagosome 161 40 2.13e-4 1.18e-2
## 7 path:bta052… Proteoglycans in cancer 197 45 6.64e-4 3.16e-2write.csv(keggRes%>%as.data.frame(),file="polledvshornedFormidbrain-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.4 tissues
1.4.1 forebrain vs frontal_skin
1.4.1.1 For horned
FBvsFS_DGE_result <- read_csv("forebrainvsfrontal_skinForhorned-ARS_jo-DEG_all.csv")FBvsFS_DGE_result[!((FBvsFS_DGE_result$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
FBvsFS_entrezID_DEG <- FBvsFS_DGE_result%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(FBvsFS_entrezID_DEG) ## [1] 20381.4.1.1.1 Go entrichment analysis
FBvsFS_goRes <- goana(FBvsFS_entrezID_DEG, ALL_entrezID, species = "Bt")
FBvsFS_goRes <- FBvsFS_goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
FBvsFS_goRes%>%print(n= dim(FBvsFS_goRes)[1])## # A tibble: 10 x 7
## goid Term Ont N DE P.DE adjP
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 GO:00055… extracellular region CC 546 116 7.56e-13 1.07e-8
## 2 GO:00076… sensory perception BP 114 36 1.86e- 9 1.31e-5
## 3 GO:00508… nervous system process BP 169 43 7.36e- 8 3.47e-4
## 4 GO:00076… visual perception BP 63 22 3.64e- 7 1.03e-3
## 5 GO:00509… sensory perception of light … BP 63 22 3.64e- 7 1.03e-3
## 6 GO:00056… extracellular space CC 243 52 1.30e- 6 3.07e-3
## 7 GO:00198… cytolysis BP 17 10 2.17e- 6 4.38e-3
## 8 GO:00513… response to mineralocorticoid BP 5 5 1.52e- 5 2.16e-2
## 9 GO:19034… response to dehydroepiandros… BP 5 5 1.52e- 5 2.16e-2
## 10 GO:19034… response to 11-deoxycorticos… BP 5 5 1.52e- 5 2.16e-2write.csv(FBvsFS_goRes%>%as.data.frame(),file="forebrainvsfrontal_skinForhorned-ARS_jo-GOresults_p0.05_lfc1.csv")1.4.1.1.2 KEGG entrichment analysis
FBvsFS_keggRes <- kegga(FBvsFS_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
FBvsFS_keggRes <- FBvsFS_keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
FBvsFS_keggRes%>%print(n= dim(FBvsFS_keggRes)[1])## # A tibble: 30 x 6
## keggid Pathway N DE P.DE adjP
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 path:bta040… Neuroactive ligand-receptor inter… 345 92 1.27e-16 4.23e-14
## 2 path:bta001… Steroid hormone biosynthesis 70 33 2.22e-14 3.69e-12
## 3 path:bta008… Retinol metabolism 65 30 6.76e-13 7.50e-11
## 4 path:bta046… Complement and coagulation cascad… 83 29 5.10e- 9 4.25e- 7
## 5 path:bta049… Protein digestion and absorption 113 33 7.19e- 8 4.79e- 6
## 6 path:bta051… Staphylococcus aureus infection 101 30 1.85e- 7 1.02e- 5
## 7 path:bta049… Bile secretion 94 28 4.39e- 7 2.09e- 5
## 8 path:bta000… Ascorbate and aldarate metabolism 23 12 1.14e- 6 4.75e- 5
## 9 path:bta005… Arachidonic acid metabolism 78 24 1.51e- 6 5.58e- 5
## 10 path:bta009… Drug metabolism - cytochrome P450 60 20 2.79e- 6 9.30e- 5
## 11 path:bta052… Chemical carcinogenesis 77 23 4.36e- 6 1.32e- 4
## 12 path:bta000… Pentose and glucuronate interconv… 30 12 3.52e- 5 9.01e- 4
## 13 path:bta009… Drug metabolism - other enzymes 75 21 3.33e- 5 9.01e- 4
## 14 path:bta049… Type I diabetes mellitus 54 16 1.36e- 4 3.22e- 3
## 15 path:bta053… Autoimmune thyroid disease 60 17 1.56e- 4 3.46e- 3
## 16 path:bta009… Metabolism of xenobiotics by cyto… 66 18 1.72e- 4 3.58e- 3
## 17 path:bta020… ABC transporters 56 16 2.17e- 4 4.25e- 3
## 18 path:bta046… Hematopoietic cell lineage 103 24 2.31e- 4 4.28e- 3
## 19 path:bta049… Ovarian steroidogenesis 57 16 2.71e- 4 4.55e- 3
## 20 path:bta050… Nicotine addiction 36 12 2.73e- 4 4.55e- 3
## 21 path:bta040… Cytokine-cytokine receptor intera… 296 52 3.37e- 4 5.34e- 3
## 22 path:bta003… Tyrosine metabolism 38 12 4.81e- 4 6.97e- 3
## 23 path:bta008… Porphyrin and chlorophyll metabol… 38 12 4.81e- 4 6.97e- 3
## 24 path:bta053… Rheumatoid arthritis 98 22 7.12e- 4 9.59e- 3
## 25 path:bta053… Graft-versus-host disease 56 15 7.20e- 4 9.59e- 3
## 26 path:bta053… Inflammatory bowel disease 69 17 9.29e- 4 1.19e- 2
## 27 path:bta007… Folate biosynthesis 36 11 1.12e- 3 1.38e- 2
## 28 path:bta049… Maturity onset diabetes of the yo… 26 9 1.17e- 3 1.39e- 2
## 29 path:bta005… Linoleic acid metabolism 33 10 2.00e- 3 2.29e- 2
## 30 path:bta046… Intestinal immune network for IgA… 54 13 4.43e- 3 4.92e- 2write.csv(FBvsFS_keggRes%>%as.data.frame(),file="forebrainvsfrontal_skinForhorned-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.4.1.2 For polled
FBvsFS_DGE_result <- read_csv("forebrainvsfrontal_skinForpolled-ARS_jo-DEG_all.csv")FBvsFS_DGE_result[!((FBvsFS_DGE_result$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
FBvsFS_entrezID_DEG <- FBvsFS_DGE_result%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(FBvsFS_entrezID_DEG) ## [1] 2141.4.1.2.1 Go entrichment analysis
FBvsFS_goRes <- goana(FBvsFS_entrezID_DEG, ALL_entrezID, species = "Bt")
FBvsFS_goRes <- FBvsFS_goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
FBvsFS_goRes%>%print(n= dim(FBvsFS_goRes)[1])## # A tibble: 0 x 7
## # … with 7 variables: goid <chr>, Term <chr>, Ont <chr>, N <dbl>, DE <dbl>,
## # P.DE <dbl>, adjP <dbl>write.csv(FBvsFS_goRes%>%as.data.frame(),file="forebrainvsfrontal_skinForpolled-ARS_jo-GOresults_p0.05_lfc1.csv")1.4.1.2.2 KEGG entrichment analysis
FBvsFS_keggRes <- kegga(FBvsFS_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
FBvsFS_keggRes <- FBvsFS_keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
FBvsFS_keggRes%>%print(n= dim(FBvsFS_keggRes)[1])## # A tibble: 7 x 6
## keggid Pathway N DE P.DE adjP
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 path:bta04740 Olfactory transduction 983 45 7.97e-16 2.66e-13
## 2 path:bta05034 Alcoholism 195 21 8.46e-15 9.39e-13
## 3 path:bta05322 Systemic lupus erythematosus 150 19 8.37e-15 9.39e-13
## 4 path:bta05203 Viral carcinogenesis 213 14 1.68e- 7 1.40e- 5
## 5 path:bta04512 ECM-receptor interaction 80 6 3.08e- 4 1.71e- 2
## 6 path:bta05014 Amyotrophic lateral sclerosis 365 13 3.06e- 4 1.71e- 2
## 7 path:bta04217 Necroptosis 158 8 4.74e- 4 2.25e- 2write.csv(FBvsFS_keggRes%>%as.data.frame(),file="forebrainvsfrontal_skinForpolled-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.4.2 forebrain vs hornbud
1.4.2.1 For horned
FBvsHB_DGE_result <- read_csv("forebrainvshornbudForhorned-ARS_jo-DEG_all.csv")FBvsHB_DGE_result[!((FBvsHB_DGE_result$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
FBvsHB_entrezID_DEG <- FBvsHB_DGE_result%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(FBvsHB_entrezID_DEG) ## [1] 37881.4.2.1.1 Go entrichment analysis
FBvsHB_goRes <- goana(FBvsHB_entrezID_DEG, ALL_entrezID, species = "Bt")
FBvsHB_goRes <- FBvsHB_goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
FBvsHB_goRes%>%print(n= dim(FBvsHB_goRes)[1])## # A tibble: 12 x 7
## goid Term Ont N DE P.DE adjP
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 GO:0005576 extracellular region CC 546 192 1.30e-16 1.84e-12
## 2 GO:0005615 extracellular space CC 243 92 1.45e-10 1.03e- 6
## 3 GO:0005882 intermediate filament CC 45 28 9.77e-10 4.61e- 6
## 4 GO:0045111 intermediate filament cytoske… CC 49 28 1.47e- 8 5.18e- 5
## 5 GO:0007600 sensory perception BP 114 47 2.43e- 7 4.91e- 4
## 6 GO:0003796 lysozyme activity MF 12 11 2.26e- 7 4.91e- 4
## 7 GO:0061783 peptidoglycan muralytic activ… MF 12 11 2.26e- 7 4.91e- 4
## 8 GO:0019835 cytolysis BP 17 13 9.77e- 7 1.73e- 3
## 9 GO:0007601 visual perception BP 63 29 3.51e- 6 4.97e- 3
## 10 GO:0050953 sensory perception of light s… BP 63 29 3.51e- 6 4.97e- 3
## 11 GO:0050877 nervous system process BP 169 59 5.87e- 6 7.55e- 3
## 12 GO:0005581 collagen trimer CC 25 15 1.55e- 5 1.83e- 2write.csv(FBvsHB_goRes%>%as.data.frame(),file="forebrainvshornbudForhorned-ARS_jo-GOresults_p0.05_lfc1.csv")1.4.2.1.2 KEGG entrichment analysis
FBvsHB_keggRes <- kegga(FBvsHB_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
FBvsHB_keggRes <- FBvsHB_keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
FBvsHB_keggRes%>%print(n= dim(FBvsHB_keggRes)[1])## # A tibble: 44 x 6
## keggid Pathway N DE P.DE adjP
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 path:bta0… Neuroactive ligand-receptor interac… 345 144 3.50e-20 1.17e-17
## 2 path:bta0… Staphylococcus aureus infection 101 59 3.87e-17 6.44e-15
## 3 path:bta0… Cytokine-cytokine receptor interact… 296 120 6.79e-16 7.54e-14
## 4 path:bta0… Steroid hormone biosynthesis 70 42 3.70e-13 3.08e-11
## 5 path:bta0… Olfactory transduction 983 279 2.04e-10 1.36e- 8
## 6 path:bta0… Graft-versus-host disease 56 32 1.35e- 9 7.47e- 8
## 7 path:bta0… Protein digestion and absorption 113 51 1.87e- 9 8.91e- 8
## 8 path:bta0… Arachidonic acid metabolism 78 39 4.17e- 9 1.73e- 7
## 9 path:bta0… Retinol metabolism 65 34 9.26e- 9 3.42e- 7
## 10 path:bta0… Complement and coagulation cascades 83 40 1.03e- 8 3.42e- 7
## 11 path:bta0… Chemical carcinogenesis 77 37 4.01e- 8 1.21e- 6
## 12 path:bta0… Linoleic acid metabolism 33 21 7.00e- 8 1.94e- 6
## 13 path:bta0… JAK-STAT signaling pathway 179 66 1.72e- 7 4.40e- 6
## 14 path:bta0… Drug metabolism - cytochrome P450 60 30 2.53e- 7 6.01e- 6
## 15 path:bta0… Nicotine addiction 36 21 5.82e- 7 1.29e- 5
## 16 path:bta0… Tyrosine metabolism 38 21 1.97e- 6 4.10e- 5
## 17 path:bta0… Metabolism of xenobiotics by cytoch… 66 30 3.26e- 6 6.38e- 5
## 18 path:bta0… Type I diabetes mellitus 54 26 3.91e- 6 7.22e- 5
## 19 path:bta0… Natural killer cell mediated cytoto… 133 49 6.48e- 6 1.14e- 4
## 20 path:bta0… Inflammatory bowel disease 69 30 9.84e- 6 1.64e- 4
## 21 path:bta0… Hematopoietic cell lineage 103 40 1.05e- 5 1.67e- 4
## 22 path:bta0… Fat digestion and absorption 49 23 2.36e- 5 3.57e- 4
## 23 path:bta0… Maturity onset diabetes of the young 26 15 2.96e- 5 4.28e- 4
## 24 path:bta0… Autoimmune thyroid disease 60 26 4.06e- 5 5.64e- 4
## 25 path:bta0… Drug metabolism - other enzymes 75 29 1.80e- 4 2.40e- 3
## 26 path:bta0… Allograft rejection 50 21 3.69e- 4 4.72e- 3
## 27 path:bta0… Intestinal immune network for IgA p… 54 22 4.51e- 4 5.56e- 3
## 28 path:bta0… Influenza A 175 54 5.45e- 4 6.48e- 3
## 29 path:bta0… Ascorbate and aldarate metabolism 23 12 6.59e- 4 6.96e- 3
## 30 path:bta0… alpha-Linolenic acid metabolism 26 13 6.69e- 4 6.96e- 3
## 31 path:bta0… Viral protein interaction with cyto… 91 32 6.37e- 4 6.96e- 3
## 32 path:bta0… Antigen processing and presentation 80 29 6.39e- 4 6.96e- 3
## 33 path:bta0… Amoebiasis 107 36 7.84e- 4 7.91e- 3
## 34 path:bta0… Bile secretion 94 32 1.20e- 3 1.17e- 2
## 35 path:bta0… Salivary secretion 91 31 1.40e- 3 1.29e- 2
## 36 path:bta0… Asthma 34 15 1.37e- 3 1.29e- 2
## 37 path:bta0… ABC transporters 56 21 2.11e- 3 1.90e- 2
## 38 path:bta0… Rheumatoid arthritis 98 32 2.57e- 3 2.25e- 2
## 39 path:bta0… Folate biosynthesis 36 15 2.72e- 3 2.33e- 2
## 40 path:bta0… Tryptophan metabolism 47 18 3.28e- 3 2.73e- 2
## 41 path:bta0… Pentose and glucuronate interconver… 30 13 3.42e- 3 2.78e- 2
## 42 path:bta0… Phenylalanine metabolism 24 11 4.12e- 3 3.27e- 2
## 43 path:bta0… Prolactin signaling pathway 83 27 5.64e- 3 4.37e- 2
## 44 path:bta0… Cell adhesion molecules 152 44 6.32e- 3 4.78e- 2write.csv(FBvsHB_keggRes%>%as.data.frame(),file="forebrainvshornbudForhorned-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.4.2.2 For polled
FBvsHB_DGE_result <- read_csv("forebrainvshornbudForpolled-ARS_jo-DEG_all.csv")FBvsHB_DGE_result[!((FBvsHB_DGE_result$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
FBvsHB_entrezID_DEG <- FBvsHB_DGE_result%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(FBvsHB_entrezID_DEG) ## [1] 1091.4.2.2.1 Go entrichment analysis
FBvsHB_goRes <- goana(FBvsHB_entrezID_DEG, ALL_entrezID, species = "Bt")
FBvsHB_goRes <- FBvsHB_goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
FBvsHB_goRes%>%print(n= dim(FBvsHB_goRes)[1])## # A tibble: 13 x 7
## goid Term Ont N DE P.DE adjP
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 GO:00055… extracellular region CC 546 16 1.08e-7 0.00152
## 2 GO:00484… platelet-derived growth factor… MF 4 3 7.65e-7 0.00541
## 3 GO:00015… ossification BP 80 6 7.19e-6 0.0260
## 4 GO:00198… growth factor binding MF 22 4 7.34e-6 0.0260
## 5 GO:00487… system development BP 788 16 1.28e-5 0.0263
## 6 GO:00422… response to chemical BP 787 16 1.26e-5 0.0263
## 7 GO:00055… protein binding MF 1383 22 1.30e-5 0.0263
## 8 GO:00485… animal organ development BP 544 13 1.67e-5 0.0296
## 9 GO:00072… multicellular organism develop… BP 920 17 2.26e-5 0.0355
## 10 GO:00708… cellular response to chemical … BP 578 13 3.15e-5 0.0366
## 11 GO:00702… protein heterotrimerization BP 2 2 3.36e-5 0.0366
## 12 GO:00331… V(D)J recombination BP 2 2 3.36e-5 0.0366
## 13 GO:00443… non-sequence-specific DNA bind… MF 2 2 3.36e-5 0.0366write.csv(FBvsHB_goRes%>%as.data.frame(),file="forebrainvshornbudForpolled-ARS_jo-GOresults_p0.05_lfc1.csv")1.4.2.2.2 KEGG entrichment analysis
FBvsHB_keggRes <- kegga(FBvsHB_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
FBvsHB_keggRes <- FBvsHB_keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
FBvsHB_keggRes%>%print(n= dim(FBvsHB_keggRes)[1])## # A tibble: 5 x 6
## keggid Pathway N DE P.DE adjP
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 path:bta05034 Alcoholism 195 17 1.39e-15 4.61e-13
## 2 path:bta05322 Systemic lupus erythematosus 150 14 2.01e-13 3.35e-11
## 3 path:bta05203 Viral carcinogenesis 213 15 1.66e-12 1.84e-10
## 4 path:bta04512 ECM-receptor interaction 80 6 7.16e- 6 5.96e- 4
## 5 path:bta04974 Protein digestion and absorption 113 5 5.19e- 4 3.46e- 2write.csv(FBvsHB_keggRes%>%as.data.frame(),file="forebrainvshornbudForpolled-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.4.3 midbrain vs forebrain
1.4.3.1 For horned
MBvsFB_DGE_result <- read_csv("midbrainvsforebrainForhorned-ARS_jo-DEG_all.csv")MBvsFB_DGE_result[!((MBvsFB_DGE_result$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
MBvsFB_entrezID_DEG <- MBvsFB_DGE_result%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(MBvsFB_entrezID_DEG) ## [1] 831.4.3.1.1 Go entrichment analysis
MBvsFB_goRes <- goana(MBvsFB_entrezID_DEG, ALL_entrezID, species = "Bt")
MBvsFB_goRes <- MBvsFB_goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
MBvsFB_goRes%>%print(n= dim(MBvsFB_goRes)[1])## # A tibble: 0 x 7
## # … with 7 variables: goid <chr>, Term <chr>, Ont <chr>, N <dbl>, DE <dbl>,
## # P.DE <dbl>, adjP <dbl>write.csv(MBvsFB_goRes%>%as.data.frame(),file="midbrainvsforebrainForhorned-ARS_jo-GOresults_p0.05_lfc1.csv")1.4.3.1.2 KEGG entrichment analysis
MBvsFB_keggRes <- kegga(MBvsFB_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
MBvsFB_keggRes <- MBvsFB_keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
MBvsFB_keggRes%>%print(n= dim(MBvsFB_keggRes)[1])## # A tibble: 0 x 6
## # … with 6 variables: keggid <chr>, Pathway <chr>, N <dbl>, DE <dbl>,
## # P.DE <dbl>, adjP <dbl>write.csv(MBvsFB_keggRes%>%as.data.frame(),file="midbrainvsforebrainForhorned-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.4.3.2 For polled
MBvsFB_DGE_result <- read_csv("midbrainvsforebrainForpolled-ARS_jo-DEG_all.csv")MBvsFB_DGE_result[!((MBvsFB_DGE_result$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
MBvsFB_entrezID_DEG <- MBvsFB_DGE_result%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(MBvsFB_entrezID_DEG) ## [1] 61.4.3.2.1 Go entrichment analysis
MBvsFB_goRes <- goana(MBvsFB_entrezID_DEG, ALL_entrezID, species = "Bt")
MBvsFB_goRes <- MBvsFB_goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
MBvsFB_goRes%>%print(n= dim(MBvsFB_goRes)[1])## # A tibble: 0 x 7
## # … with 7 variables: goid <chr>, Term <chr>, Ont <chr>, N <dbl>, DE <dbl>,
## # P.DE <dbl>, adjP <dbl>write.csv(MBvsFB_goRes%>%as.data.frame(),file="midbrainvsforebrainForpolled-ARS_jo-GOresults_p0.05_lfc1.csv")1.4.3.2.2 KEGG entrichment analysis
MBvsFB_keggRes <- kegga(FBvsHB_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
MBvsFB_keggRes <- MBvsFB_keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
MBvsFB_keggRes%>%print(n= dim(MBvsFB_keggRes)[1])## # A tibble: 5 x 6
## keggid Pathway N DE P.DE adjP
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 path:bta05034 Alcoholism 195 17 1.39e-15 4.61e-13
## 2 path:bta05322 Systemic lupus erythematosus 150 14 2.01e-13 3.35e-11
## 3 path:bta05203 Viral carcinogenesis 213 15 1.66e-12 1.84e-10
## 4 path:bta04512 ECM-receptor interaction 80 6 7.16e- 6 5.96e- 4
## 5 path:bta04974 Protein digestion and absorption 113 5 5.19e- 4 3.46e- 2write.csv(MBvsFB_keggRes%>%as.data.frame(),file="midbrainvsforebrainForpolled-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.4.4 hornbud vs frontal_skin
1.4.4.1 For horned
HBvsFS_DGE_result <- read_csv("hornbudvsfrontal_skinForhorned-ARS_jo-DEG_all.csv")HBvsFS_DGE_result[!((HBvsFS_DGE_result$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
HBvsFS_entrezID_DEG <- HBvsFS_DGE_result%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(HBvsFS_entrezID_DEG) ## [1] 901.4.4.1.1 Go entrichment analysis
HBvsFS_goRes <- goana(HBvsFS_entrezID_DEG, ALL_entrezID, species = "Bt")
HBvsFS_goRes <- HBvsFS_goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
HBvsFS_goRes%>%print(n= dim(HBvsFS_goRes)[1])## # A tibble: 0 x 7
## # … with 7 variables: goid <chr>, Term <chr>, Ont <chr>, N <dbl>, DE <dbl>,
## # P.DE <dbl>, adjP <dbl>write.csv(HBvsFS_goRes%>%as.data.frame(),file="hornbudvsfrontal_skinForhorned-ARS_jo-GOresults_p0.05_lfc1.csv")1.4.4.1.2 KEGG entrichment analysis
HBvsFS_keggRes <- kegga(HBvsFS_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
HBvsFS_keggRes <- HBvsFS_keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
HBvsFS_keggRes%>%print(n= dim(HBvsFS_keggRes)[1])## # A tibble: 0 x 6
## # … with 6 variables: keggid <chr>, Pathway <chr>, N <dbl>, DE <dbl>,
## # P.DE <dbl>, adjP <dbl>write.csv(HBvsFS_keggRes%>%as.data.frame(),file="hornbudvsfrontal_skinForhorned-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.4.4.2 For polled
HBvsFS_DGE_result <- read_csv("hornbudvsfrontal_skinForpolled-ARS_jo-DEG_all.csv")HBvsFS_DGE_result[!((HBvsFS_DGE_result$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
HBvsFS_entrezID_DEG <- HBvsFS_DGE_result%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(HBvsFS_entrezID_DEG) ## [1] 11.4.4.2.1 Go entrichment analysis
HBvsFS_goRes <- goana(HBvsFS_entrezID_DEG, ALL_entrezID, species = "Bt")
HBvsFS_goRes <- HBvsFS_goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
HBvsFS_goRes%>%print(n= dim(HBvsFS_goRes)[1])## # A tibble: 0 x 7
## # … with 7 variables: goid <chr>, Term <chr>, Ont <chr>, N <dbl>, DE <dbl>,
## # P.DE <dbl>, adjP <dbl>write.csv(HBvsFS_goRes%>%as.data.frame(),file="hornbudvsfrontal_skinForpolled-ARS_jo-GOresults_p0.05_lfc1.csv")1.4.4.2.2 KEGG entrichment analysis
HBvsFS_keggRes <- kegga(HBvsFS_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
HBvsFS_keggRes <- HBvsFS_keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
HBvsFS_keggRes%>%print(n= dim(HBvsFS_keggRes)[1])## # A tibble: 0 x 6
## # … with 6 variables: keggid <chr>, Pathway <chr>, N <dbl>, DE <dbl>,
## # P.DE <dbl>, adjP <dbl>write.csv(HBvsFS_keggRes%>%as.data.frame(),file="hornbudvsfrontal_skinForpolled-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.4.5 midbrain vs frontal_skin
1.4.5.1 For horned
MBvsFS_DGE_result <- read_csv("midbrainvsfrontal_skinForhorned-ARS_jo-DEG_all.csv")MBvsFS_DGE_result[!((MBvsFS_DGE_result$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
MBvsFS_entrezID_DEG <- MBvsFS_DGE_result%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(MBvsFS_entrezID_DEG) ## [1] 30761.4.5.1.1 Go entrichment analysis
MBvsFS_goRes <- goana(MBvsFS_entrezID_DEG, ALL_entrezID, species = "Bt")
MBvsFS_goRes <- MBvsFS_goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
MBvsFS_goRes%>%print(n= dim(MBvsFS_goRes)[1])## # A tibble: 5 x 7
## goid Term Ont N DE P.DE adjP
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 GO:0005576 extracellular region CC 546 160 1.84e-14 2.60e-10
## 2 GO:0007600 sensory perception BP 114 43 3.54e- 8 1.67e- 4
## 3 GO:0005615 extracellular space CC 243 74 3.54e- 8 1.67e- 4
## 4 GO:0050877 nervous system process BP 169 51 5.96e- 6 2.10e- 2
## 5 GO:0007218 neuropeptide signaling pathway BP 25 14 7.42e- 6 2.10e- 2write.csv(MBvsFS_goRes%>%as.data.frame(),file="midbrainvsfrontal_skinForhorned-ARS_jo-GOresults_p0.05_lfc1.csv")1.4.5.1.2 KEGG entrichment analysis
MBvsFS_keggRes <- kegga(MBvsFS_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
MBvsFS_keggRes <- MBvsFS_keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
MBvsFS_keggRes%>%print(n= dim(MBvsFS_keggRes)[1])## # A tibble: 43 x 6
## keggid Pathway N DE P.DE adjP
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 path:bta04… Neuroactive ligand-receptor intera… 345 139 1.30e-26 4.32e-24
## 2 path:bta00… Steroid hormone biosynthesis 70 40 1.05e-14 1.76e-12
## 3 path:bta00… Retinol metabolism 65 38 1.76e-14 1.95e-12
## 4 path:bta04… Complement and coagulation cascades 83 43 1.16e-13 9.64e-12
## 5 path:bta04… Protein digestion and absorption 113 47 1.78e-10 1.19e- 8
## 6 path:bta05… Staphylococcus aureus infection 101 43 4.34e-10 2.41e- 8
## 7 path:bta05… Chemical carcinogenesis 77 35 2.18e- 9 1.04e- 7
## 8 path:bta00… Drug metabolism - cytochrome P450 60 29 9.08e- 9 3.78e- 7
## 9 path:bta05… Nicotine addiction 36 21 1.39e- 8 5.16e- 7
## 10 path:bta00… Arachidonic acid metabolism 78 33 5.44e- 8 1.81e- 6
## 11 path:bta04… Bile secretion 94 37 8.52e- 8 2.58e- 6
## 12 path:bta05… Graft-versus-host disease 56 26 1.47e- 7 4.09e- 6
## 13 path:bta04… Cytokine-cytokine receptor interac… 296 83 2.91e- 7 7.46e- 6
## 14 path:bta00… Metabolism of xenobiotics by cytoc… 66 28 5.10e- 7 1.21e- 5
## 15 path:bta00… Linoleic acid metabolism 33 18 6.08e- 7 1.35e- 5
## 16 path:bta00… Drug metabolism - other enzymes 75 30 9.27e- 7 1.93e- 5
## 17 path:bta04… Hematopoietic cell lineage 103 36 3.59e- 6 7.03e- 5
## 18 path:bta00… Tyrosine metabolism 38 18 8.57e- 6 1.58e- 4
## 19 path:bta04… Fat digestion and absorption 49 21 1.10e- 5 1.85e- 4
## 20 path:bta05… Autoimmune thyroid disease 60 24 1.11e- 5 1.85e- 4
## 21 path:bta00… Ascorbate and aldarate metabolism 23 13 1.37e- 5 2.18e- 4
## 22 path:bta05… Inflammatory bowel disease 69 26 1.73e- 5 2.62e- 4
## 23 path:bta04… Type I diabetes mellitus 54 22 1.82e- 5 2.63e- 4
## 24 path:bta04… Natural killer cell mediated cytot… 133 41 2.66e- 5 3.69e- 4
## 25 path:bta02… ABC transporters 56 22 3.57e- 5 4.76e- 4
## 26 path:bta04… Ovarian steroidogenesis 57 22 4.92e- 5 6.30e- 4
## 27 path:bta00… Pentose and glucuronate interconve… 30 14 1.06e- 4 1.31e- 3
## 28 path:bta05… Rheumatoid arthritis 98 30 3.47e- 4 4.13e- 3
## 29 path:bta04… Maturity onset diabetes of the you… 26 12 3.75e- 4 4.31e- 3
## 30 path:bta04… Prolactin signaling pathway 83 26 5.67e- 4 6.29e- 3
## 31 path:bta04… Cell adhesion molecules 152 41 6.52e- 4 6.79e- 3
## 32 path:bta05… Allograft rejection 50 18 6.32e- 4 6.79e- 3
## 33 path:bta04… Inflammatory mediator regulation o… 97 28 1.48e- 3 1.49e- 2
## 34 path:bta00… alpha-Linolenic acid metabolism 26 11 1.59e- 3 1.56e- 2
## 35 path:bta04… Intestinal immune network for IgA … 54 18 1.76e- 3 1.67e- 2
## 36 path:bta04… Salivary secretion 91 26 2.51e- 3 2.32e- 2
## 37 path:bta04… Th17 cell differentiation 110 30 2.70e- 3 2.43e- 2
## 38 path:bta00… Glycine, serine and threonine meta… 44 15 3.24e- 3 2.76e- 2
## 39 path:bta04… JAK-STAT signaling pathway 179 44 3.18e- 3 2.76e- 2
## 40 path:bta04… Calcium signaling pathway 195 47 3.56e- 3 2.96e- 2
## 41 path:bta04… Pancreatic secretion 100 27 4.97e- 3 4.04e- 2
## 42 path:bta00… Porphyrin and chlorophyll metaboli… 38 13 5.75e- 3 4.56e- 2
## 43 path:bta05… Asthma 34 12 5.92e- 3 4.59e- 2write.csv(MBvsFS_keggRes%>%as.data.frame(),file="midbrainvsfrontal_skinForhorned-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.4.5.2 For polled
MBvsFS_DGE_result <- read_csv("midbrainvsfrontal_skinForpolled-ARS_jo-DEG_all.csv")MBvsFS_DGE_result[!((MBvsFS_DGE_result$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
MBvsFS_entrezID_DEG <- MBvsFS_DGE_result%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(MBvsFS_entrezID_DEG) ## [1] 2401.4.5.2.1 Go entrichment analysis
MBvsFS_goRes <- goana(MBvsFS_entrezID_DEG, ALL_entrezID, species = "Bt")
MBvsFS_goRes <- MBvsFS_goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
MBvsFS_goRes%>%print(n= dim(MBvsFS_goRes)[1])## # A tibble: 20 x 7
## goid Term Ont N DE P.DE adjP
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 GO:00015… ossification BP 80 11 5.89e-9 8.34e-5
## 2 GO:00312… biomineral tissue development BP 44 8 7.74e-8 3.65e-4
## 3 GO:01101… biomineralization BP 44 8 7.74e-8 3.65e-4
## 4 GO:00015… skeletal system development BP 78 9 6.82e-7 2.41e-3
## 5 GO:00426… muscle cell differentiation BP 62 8 1.21e-6 3.43e-3
## 6 GO:00512… cartilage development BP 36 6 5.89e-6 1.39e-2
## 7 GO:00055… extracellular region CC 546 21 8.04e-6 1.46e-2
## 8 GO:00484… platelet-derived growth fact… MF 4 3 8.25e-6 1.46e-2
## 9 GO:00610… muscle structure development BP 113 9 1.51e-5 1.95e-2
## 10 GO:00302… regulation of ossification BP 41 6 1.29e-5 1.95e-2
## 11 GO:00076… sensory perception of sound BP 25 5 1.43e-5 1.95e-2
## 12 GO:00098… animal organ morphogenesis BP 145 10 1.81e-5 2.01e-2
## 13 GO:00509… sensory perception of mechan… BP 27 5 2.13e-5 2.01e-2
## 14 GO:00019… endochondral ossification BP 5 3 2.04e-5 2.01e-2
## 15 GO:00360… replacement ossification BP 5 3 2.04e-5 2.01e-2
## 16 GO:00485… animal organ development BP 544 20 2.54e-5 2.15e-2
## 17 GO:00098… tissue development BP 331 15 2.59e-5 2.15e-2
## 18 GO:00511… striated muscle cell differe… BP 49 6 3.68e-5 2.60e-2
## 19 GO:00302… bone mineralization BP 30 5 3.65e-5 2.60e-2
## 20 GO:00614… connective tissue development BP 49 6 3.68e-5 2.60e-2write.csv(MBvsFS_goRes%>%as.data.frame(),file="midbrainvsfrontal_skinForpolled-ARS_jo-GOresults_p0.05_lfc1.csv")1.4.5.2.2 KEGG entrichment analysis
MBvsFS_keggRes <- kegga(MBvsFS_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
MBvsFS_keggRes <- MBvsFS_keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
MBvsFS_keggRes%>%print(n= dim(MBvsFS_keggRes)[1])## # A tibble: 4 x 6
## keggid Pathway N DE P.DE adjP
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 path:bta04740 Olfactory transduction 983 59 5.26e-24 1.75e-21
## 2 path:bta04512 ECM-receptor interaction 80 10 7.43e- 8 1.24e- 5
## 3 path:bta04974 Protein digestion and absorption 113 9 1.50e- 5 1.67e- 3
## 4 path:bta04510 Focal adhesion 192 11 3.89e- 5 3.24e- 3write.csv(MBvsFS_keggRes%>%as.data.frame(),file="midbrainvsfrontal_skinForpolled-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.4.6 midbrain vs hornbud
1.4.6.1 For horned
MBvsHB_DGE_result <- read_csv("midbrainvshornbudForhorned-ARS_jo-DEG_all.csv")MBvsHB_DGE_result[!((MBvsHB_DGE_result$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
MBvsHB_entrezID_DEG <- MBvsHB_DGE_result%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(MBvsHB_entrezID_DEG) ## [1] 45851.4.6.1.1 Go entrichment analysis
MBvsHB_goRes <- goana(MBvsHB_entrezID_DEG, ALL_entrezID, species = "Bt")
MBvsHB_goRes <- MBvsHB_goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
MBvsHB_goRes%>%print(n= dim(MBvsHB_goRes)[1])## # A tibble: 11 x 7
## goid Term Ont N DE P.DE adjP
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 GO:0005576 extracellular region CC 546 221 4.89e-17 6.92e-13
## 2 GO:0005882 intermediate filament CC 45 32 5.75e-11 4.07e- 7
## 3 GO:0005615 extracellular space CC 243 105 1.01e-10 4.77e- 7
## 4 GO:0045111 intermediate filament cytoske… CC 49 33 2.86e-10 1.01e- 6
## 5 GO:0007601 visual perception BP 63 35 1.26e- 7 2.97e- 4
## 6 GO:0050953 sensory perception of light s… BP 63 35 1.26e- 7 2.97e- 4
## 7 GO:0007600 sensory perception BP 114 53 2.52e- 7 5.10e- 4
## 8 GO:0019835 cytolysis BP 17 14 8.64e- 7 1.53e- 3
## 9 GO:0050877 nervous system process BP 169 69 1.94e- 6 3.04e- 3
## 10 GO:0003796 lysozyme activity MF 12 10 3.08e- 5 3.97e- 2
## 11 GO:0061783 peptidoglycan muralytic activ… MF 12 10 3.08e- 5 3.97e- 2write.csv(MBvsHB_goRes%>%as.data.frame(),file="midbrainvshornbudForhorned-ARS_jo-GOresults_p0.05_lfc1.csv")1.4.6.1.2 KEGG entrichment analysis
MBvsHB_keggRes <- kegga(MBvsHB_entrezID_DEG, ALL_entrezID, species = "Bt")
# for KEGG pathway
MBvsHB_keggRes <- MBvsHB_keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
MBvsHB_keggRes%>%print(n= dim(MBvsHB_keggRes)[1])## # A tibble: 51 x 6
## keggid Pathway N DE P.DE adjP
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 path:bta0… Neuroactive ligand-receptor interac… 345 175 2.64e-26 8.79e-24
## 2 path:bta0… Staphylococcus aureus infection 101 68 1.03e-19 1.71e-17
## 3 path:bta0… Cytokine-cytokine receptor interact… 296 135 1.23e-15 1.37e-13
## 4 path:bta0… Steroid hormone biosynthesis 70 47 5.49e-14 4.57e-12
## 5 path:bta0… Protein digestion and absorption 113 60 5.48e-11 3.04e- 9
## 6 path:bta0… Chemical carcinogenesis 77 46 4.66e-11 3.04e- 9
## 7 path:bta0… Retinol metabolism 65 40 2.38e-10 1.13e- 8
## 8 path:bta0… Natural killer cell mediated cytoto… 133 66 2.90e-10 1.18e- 8
## 9 path:bta0… Graft-versus-host disease 56 36 3.20e-10 1.18e- 8
## 10 path:bta0… Drug metabolism - cytochrome P450 60 37 1.03e- 9 3.43e- 8
## 11 path:bta0… Nicotine addiction 36 26 2.14e- 9 6.49e- 8
## 12 path:bta0… Arachidonic acid metabolism 78 43 6.32e- 9 1.76e- 7
## 13 path:bta0… Complement and coagulation cascades 83 44 2.09e- 8 5.34e- 7
## 14 path:bta0… Olfactory transduction 983 314 3.26e- 8 7.75e- 7
## 15 path:bta0… Metabolism of xenobiotics by cytoch… 66 37 3.93e- 8 8.72e- 7
## 16 path:bta0… Linoleic acid metabolism 33 23 5.47e- 8 1.14e- 6
## 17 path:bta0… Fat digestion and absorption 49 29 2.36e- 7 4.63e- 6
## 18 path:bta0… Type I diabetes mellitus 54 29 3.72e- 6 6.88e- 5
## 19 path:bta0… Hematopoietic cell lineage 103 46 5.83e- 6 1.02e- 4
## 20 path:bta0… Inflammatory bowel disease 69 34 7.04e- 6 1.17e- 4
## 21 path:bta0… JAK-STAT signaling pathway 179 70 9.58e- 6 1.52e- 4
## 22 path:bta0… Tyrosine metabolism 38 22 1.12e- 5 1.69e- 4
## 23 path:bta0… Viral protein interaction with cyto… 91 41 1.42e- 5 2.00e- 4
## 24 path:bta0… Bile secretion 94 42 1.44e- 5 2.00e- 4
## 25 path:bta0… Autoimmune thyroid disease 60 30 1.67e- 5 2.22e- 4
## 26 path:bta0… Drug metabolism - other enzymes 75 35 2.34e- 5 3.00e- 4
## 27 path:bta0… Allograft rejection 50 26 2.50e- 5 3.09e- 4
## 28 path:bta0… Salivary secretion 91 40 3.62e- 5 4.31e- 4
## 29 path:bta0… Maturity onset diabetes of the young 26 16 6.45e- 5 7.41e- 4
## 30 path:bta0… Amoebiasis 107 44 1.06e- 4 1.18e- 3
## 31 path:bta0… Ascorbate and aldarate metabolism 23 14 2.21e- 4 2.37e- 3
## 32 path:bta0… Rheumatoid arthritis 98 40 2.58e- 4 2.68e- 3
## 33 path:bta0… alpha-Linolenic acid metabolism 26 15 3.01e- 4 3.04e- 3
## 34 path:bta0… Antigen processing and presentation 80 33 6.88e- 4 6.74e- 3
## 35 path:bta0… ABC transporters 56 25 7.60e- 4 7.23e- 3
## 36 path:bta0… Estrogen signaling pathway 129 48 8.43e- 4 7.80e- 3
## 37 path:bta0… Pancreatic secretion 100 39 8.89e- 4 8.00e- 3
## 38 path:bta0… Intestinal immune network for IgA p… 54 24 1.04e- 3 8.89e- 3
## 39 path:bta0… Ovarian steroidogenesis 57 25 1.04e- 3 8.89e- 3
## 40 path:bta0… Asthma 34 17 1.12e- 3 9.37e- 3
## 41 path:bta0… Ether lipid metabolism 50 22 1.94e- 3 1.52e- 2
## 42 path:bta0… Calcium signaling pathway 195 66 1.97e- 3 1.52e- 2
## 43 path:bta0… Tuberculosis 188 64 1.94e- 3 1.52e- 2
## 44 path:bta0… Pentose and glucuronate interconver… 30 15 2.18e- 3 1.65e- 2
## 45 path:bta0… Primary immunodeficiency 39 18 2.57e- 3 1.90e- 2
## 46 path:bta0… Cell adhesion molecules 152 52 4.34e- 3 3.08e- 2
## 47 path:bta0… Phototransduction 26 13 4.27e- 3 3.08e- 2
## 48 path:bta0… cAMP signaling pathway 219 71 4.64e- 3 3.22e- 2
## 49 path:bta0… Phenylalanine metabolism 24 12 5.98e- 3 4.05e- 2
## 50 path:bta0… Insulin secretion 83 31 6.08e- 3 4.05e- 2
## 51 path:bta0… Carbohydrate digestion and absorpti… 42 18 6.77e- 3 4.42e- 2write.csv(MBvsHB_keggRes%>%as.data.frame(),file="midbrainvshornbudForhorned-ARS_jo-KEGGresults_p0.05_lfc1.csv")1.4.6.2 For polled
MBvsHB_DGE_result <- read_csv("midbrainvshornbudForpolled-ARS_jo-DEG_all.csv")MBvsHB_DGE_result[!((MBvsHB_DGE_result$logFC)%>%is.na()),]%>%
mutate(Sig = adj.P.Val < 0.05 & abs(logFC) > 1,
Direction = case_when(
Sig & logFC > 1 ~ "Up",
Sig & logFC < -1 ~ "Down",
!Sig ~ "Unchanged"))%>%
ggplot(aes(logFC, -log10(P.Value), colour = Direction)) +
geom_point() +
scale_colour_manual(values = c("blue", "grey", "red")) +
labs(y = "-log10(p)") +
theme_bw()
MBvsHB_entrezID_DEG <- MBvsHB_DGE_result%>%
dplyr::filter(adj.P.Val < 0.05 & abs(logFC) > 1)%>%
as.data.frame()%>%
left_join(data_genesGR, by = c("gene_id" = "gene_id"))%>%
.[["entrezid"]] %>%
unlist() %>%
unique() %>%
subset(!is.na(.))
length(MBvsHB_entrezID_DEG) ## [1] 1931.4.6.2.1 Go entrichment analysis
MBvsHB_goRes <- goana(MBvsHB_entrezID_DEG, ALL_entrezID, species = "Bt")
MBvsHB_goRes <- MBvsHB_goRes %>%
rownames_to_column("goid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
MBvsHB_goRes%>%print(n= dim(MBvsHB_goRes)[1])## # A tibble: 36 x 7
## goid Term Ont N DE P.DE adjP
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 GO:0001… ossification BP 80 11 5.94e-10 2.80e-6
## 2 GO:0031… biomineral tissue development BP 44 9 5.68e-10 2.80e-6
## 3 GO:0110… biomineralization BP 44 9 5.68e-10 2.80e-6
## 4 GO:0005… extracellular region CC 546 24 2.29e- 9 8.11e-6
## 5 GO:0009… tissue development BP 331 16 3.53e- 7 9.99e-4
## 6 GO:0030… bone mineralization BP 30 6 5.37e- 7 1.27e-3
## 7 GO:0032… developmental process BP 1123 29 4.49e- 6 5.29e-3
## 8 GO:0050… epithelial cell proliferation BP 89 8 3.82e- 6 5.29e-3
## 9 GO:0048… animal organ development BP 544 19 3.76e- 6 5.29e-3
## 10 GO:0009… animal organ morphogenesis BP 145 10 2.66e- 6 5.29e-3
## 11 GO:0030… regulation of ossification BP 41 6 3.70e- 6 5.29e-3
## 12 GO:0048… platelet-derived growth factor … MF 4 3 4.29e- 6 5.29e-3
## 13 GO:0048… system development BP 788 23 6.94e- 6 7.41e-3
## 14 GO:0048… anatomical structure development BP 1028 27 7.33e- 6 7.41e-3
## 15 GO:0007… multicellular organism developm… BP 920 25 9.28e- 6 8.29e-3
## 16 GO:0040… growth BP 167 10 9.37e- 6 8.29e-3
## 17 GO:0001… endochondral ossification BP 5 3 1.06e- 5 8.36e-3
## 18 GO:0036… replacement ossification BP 5 3 1.06e- 5 8.36e-3
## 19 GO:0032… multicellular organismal process BP 1326 31 1.49e- 5 1.11e-2
## 20 GO:0009… anatomical structure morphogene… BP 445 16 1.59e- 5 1.12e-2
## 21 GO:0042… regulation of cell population p… BP 267 12 2.22e- 5 1.44e-2
## 22 GO:0061… muscle structure development BP 113 8 2.25e- 5 1.44e-2
## 23 GO:0030… regulation of bone mineralizati… BP 19 4 3.76e- 5 2.31e-2
## 24 GO:0030… extracellular matrix organizati… BP 39 5 4.79e- 5 2.51e-2
## 25 GO:0070… regulation of biomineral tissue… BP 20 4 4.67e- 5 2.51e-2
## 26 GO:0110… regulation of biomineralization BP 20 4 4.67e- 5 2.51e-2
## 27 GO:0005… extracellular space CC 243 11 4.59e- 5 2.51e-2
## 28 GO:0043… extracellular structure organiz… BP 40 5 5.43e- 5 2.74e-2
## 29 GO:0008… cell population proliferation BP 344 13 6.13e- 5 2.99e-2
## 30 GO:0032… collagen metabolic process BP 23 4 8.32e- 5 3.92e-2
## 31 GO:0060… endochondral bone morphogenesis BP 9 3 8.66e- 5 3.95e-2
## 32 GO:0045… positive regulation of growth BP 45 5 9.68e- 5 4.03e-2
## 33 GO:0008… negative regulation of cell pop… BP 103 7 9.56e- 5 4.03e-2
## 34 GO:0048… developmental growth BP 103 7 9.56e- 5 4.03e-2
## 35 GO:0070… protein heterotrimerization BP 2 2 1.06e- 4 4.16e-2
## 36 GO:0045… gamma-catenin binding MF 2 2 1.06e- 4 4.16e-2write.csv(MBvsHB_goRes%>%as.data.frame(),file="midbrainvshornbudForpolled-ARS_jo-GOresults_p0.05_lfc1.csv")1.4.6.2.2 KEGG entrichment analysis
MBvsHB_keggRes <- kegga(MBvsHB_entrezID_DEG, ALL_entrezID, species = "Bt") # for KEGG pathway
MBvsHB_keggRes <- MBvsHB_keggRes %>%
rownames_to_column("keggid") %>%
mutate(adjP = p.adjust(P.DE, method = "fdr"))%>%
subset(adjP < 0.05) %>%
arrange(adjP) %>%
as_tibble()
MBvsHB_keggRes%>%print(n= dim(MBvsHB_keggRes)[1])## # A tibble: 6 x 6
## keggid Pathway N DE P.DE adjP
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 path:bta04512 ECM-receptor interaction 80 11 5.90e-10 1.96e-7
## 2 path:bta04974 Protein digestion and absorption 113 10 2.65e- 7 4.41e-5
## 3 path:bta04510 Focal adhesion 192 11 5.06e- 6 5.62e-4
## 4 path:bta04151 PI3K-Akt signaling pathway 354 13 8.13e- 5 6.77e-3
## 5 path:bta05165 Human papillomavirus infection 329 12 1.63e- 4 1.09e-2
## 6 path:bta04740 Olfactory transduction 983 23 2.06e- 4 1.14e-2write.csv(MBvsHB_keggRes%>%as.data.frame(),file="midbrainvshornbudForpolled-ARS_jo-KEGGresults_p0.05_lfc1.csv")2 Appendix
sessionInfo()## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.4
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ensembldb_2.12.1 AnnotationFilter_1.12.0 GenomicFeatures_1.40.1
## [4] AnnotationDbi_1.50.3 Biobase_2.48.0 GenomicRanges_1.40.0
## [7] GenomeInfoDb_1.24.2 IRanges_2.22.2 S4Vectors_0.26.1
## [10] tidyr_1.1.1 GOfuncR_1.8.0 vioplot_0.3.5
## [13] zoo_1.8-8 sm_2.2-5.6 biomaRt_2.44.1
## [16] magrittr_1.5 AnnotationHub_2.20.2 BiocFileCache_1.12.1
## [19] dbplyr_1.4.4 BiocGenerics_0.34.0 tibble_3.0.3
## [22] limma_3.44.3 ggrepel_0.8.2 pander_0.6.3
## [25] dplyr_1.0.2 ggplot2_3.3.2 readr_1.3.1
##
## loaded via a namespace (and not attached):
## [1] ProtGenerics_1.20.0 matrixStats_0.56.0
## [3] bitops_1.0-6 bit64_4.0.2
## [5] progress_1.2.2 httr_1.4.2
## [7] tools_4.0.2 utf8_1.1.4
## [9] R6_2.4.1 lazyeval_0.2.2
## [11] DBI_1.1.0 colorspace_1.4-1
## [13] withr_2.2.0 tidyselect_1.1.0
## [15] prettyunits_1.1.1 bit_4.0.4
## [17] curl_4.3 compiler_4.0.2
## [19] cli_2.0.2 DelayedArray_0.14.1
## [21] labeling_0.3 rtracklayer_1.48.0
## [23] bookdown_0.20 org.Bt.eg.db_3.11.4
## [25] scales_1.1.1 askpass_1.1
## [27] rappdirs_0.3.1 stringr_1.4.0
## [29] digest_0.6.25 Rsamtools_2.4.0
## [31] rmarkdown_2.3 XVector_0.28.0
## [33] pkgconfig_2.0.3 htmltools_0.5.0
## [35] fastmap_1.0.1 rlang_0.4.7
## [37] RSQLite_2.2.0 shiny_1.5.0
## [39] farver_2.0.3 mapplots_1.5.1
## [41] generics_0.0.2 BiocParallel_1.22.0
## [43] gtools_3.8.2 RCurl_1.98-1.2
## [45] GO.db_3.11.4 GenomeInfoDbData_1.2.3
## [47] Matrix_1.2-18 fansi_0.4.1
## [49] Rcpp_1.0.5 munsell_0.5.0
## [51] lifecycle_0.2.0 stringi_1.4.6
## [53] yaml_2.2.1 SummarizedExperiment_1.18.2
## [55] zlibbioc_1.34.0 grid_4.0.2
## [57] blob_1.2.1 promises_1.1.1
## [59] crayon_1.3.4 lattice_0.20-41
## [61] Biostrings_2.56.0 hms_0.5.3
## [63] knitr_1.29 pillar_1.4.6
## [65] tcltk_4.0.2 XML_3.99-0.5
## [67] glue_1.4.1 BiocVersion_3.11.1
## [69] evaluate_0.14 BiocManager_1.30.10
## [71] vctrs_0.3.2 httpuv_1.5.4
## [73] gtable_0.3.0 openssl_1.4.2
## [75] purrr_0.3.4 assertthat_0.2.1
## [77] xfun_0.16 mime_0.9
## [79] xtable_1.8-4 later_1.1.0.1
## [81] GenomicAlignments_1.24.0 memoise_1.1.0
## [83] ellipsis_0.3.1 interactiveDisplayBase_1.26.3