Counting.genes.all.data

# CHECK R VERSION
 stopifnot(R.version.string == "R version 3.5.3 (2019-03-11)")

# CLEANUP
# Clear all variables:
rm(list=ls(all=TRUE))

# Unload current packages:
if (!is.null(names(utils::sessionInfo()[["otherPkgs"]]))) pacman::p_unload("all") 

# KNITR DEFAULTS
knitr::opts_chunk$set(eval = TRUE, fig.height = 10)

Load packages

library(magrittr)
library(dplyr)
library(rstudioapi)

# GET CURRENT SCRIPT NAME
(this.script <- rstudioapi::getActiveDocumentContext() %>% .$path %>% basename)
getwd()
list.files()
stopifnot(this.script != "")

Data files

Upregulated genes

data.files.up <- list.files(path="data.out", pattern="genes\\.up\\.uniq.*", full.names = TRUE)
data.files.up

[1] "data.out/A-BUGS-23-DC.analysis.rmd.genes.up.uniq.txt"            
[2] "data.out/A-BUGS-23-MP.analysis.rmd.genes.up.uniq.txt"            
[3] "data.out/GSE103092.analysis.rmd.genes.up.uniq.txt"               
[4] "data.out/GSE108363.bcg.discovery.analysis.Rmd.genes.up.uniq.txt" 
[5] "data.out/GSE108363.bcg.validation.analysis.Rmd.genes.up.uniq.txt"
[6] "data.out/GSE108363.mtb.discovery.analysis.Rmd.genes.up.uniq.txt" 
[7] "data.out/GSE108363.mtb.validation.analysis.Rmd.genes.up.uniq.txt"
[8] "data.out/GSE34151.analysis.rmd.genes.up.uniq.txt"                

basename(data.files.up)

[1] "A-BUGS-23-DC.analysis.rmd.genes.up.uniq.txt"            
[2] "A-BUGS-23-MP.analysis.rmd.genes.up.uniq.txt"            
[3] "GSE103092.analysis.rmd.genes.up.uniq.txt"               
[4] "GSE108363.bcg.discovery.analysis.Rmd.genes.up.uniq.txt" 
[5] "GSE108363.bcg.validation.analysis.Rmd.genes.up.uniq.txt"
[6] "GSE108363.mtb.discovery.analysis.Rmd.genes.up.uniq.txt" 
[7] "GSE108363.mtb.validation.analysis.Rmd.genes.up.uniq.txt"
[8] "GSE34151.analysis.rmd.genes.up.uniq.txt"                

Downregulated genes

data.files.down <- list.files(path="data.out", pattern="genes\\.down\\.uniq.*", full.names = T)
basename(data.files.down)

[1] "A-BUGS-23-DC.analysis.rmd.genes.down.uniq.txt"            
[2] "A-BUGS-23-MP.analysis.rmd.genes.down.uniq.txt"            
[3] "GSE103092.analysis.rmd.genes.down.uniq.txt"               
[4] "GSE108363.bcg.discovery.analysis.Rmd.genes.down.uniq.txt" 
[5] "GSE108363.bcg.validation.analysis.Rmd.genes.down.uniq.txt"
[6] "GSE108363.mtb.discovery.analysis.Rmd.genes.down.uniq.txt" 
[7] "GSE108363.mtb.validation.analysis.Rmd.genes.down.uniq.txt"
[8] "GSE34151.analysis.rmd.genes.down.uniq.txt"                

Load data: upregulated genes

genes.up <- lapply(data.files.up, function(f) {
  scan(file = f, what = "character")
}) %>% 
  magrittr::set_names(basename(data.files.up))

Read 646 items
Read 515 items
Read 8 items
Read 193 items
Read 379 items
Read 57 items
Read 265 items
Read 665 items

# class(genes.up) # "list"
# length(genes.up) # 9
# sapply(genes.up, class)
# sapply(genes.up, length)

data.frame(File=names(genes.up),
            Gene.count=sapply(genes.up, length)) %>% 
  knitr::kable(., row.names = F)

File	Gene.count
A-BUGS-23-DC.analysis.rmd.genes.up.uniq.txt	646
A-BUGS-23-MP.analysis.rmd.genes.up.uniq.txt	515
GSE103092.analysis.rmd.genes.up.uniq.txt	8
GSE108363.bcg.discovery.analysis.Rmd.genes.up.uniq.txt	193
GSE108363.bcg.validation.analysis.Rmd.genes.up.uniq.txt	379
GSE108363.mtb.discovery.analysis.Rmd.genes.up.uniq.txt	57
GSE108363.mtb.validation.analysis.Rmd.genes.up.uniq.txt	265
GSE34151.analysis.rmd.genes.up.uniq.txt	665

Load data: downregulated genes

genes.down <- lapply(data.files.down, function(f) {
  scan(file = f, what = "character")
}) %>% 
  magrittr::set_names(basename(data.files.down))

Read 792 items
Read 802 items
Read 48 items
Read 277 items
Read 225 items
Read 43 items
Read 208 items
Read 741 items

data.frame(File=names(genes.down),
            Gene.count=sapply(genes.down, length)) %>% 
  knitr::kable(., row.names = F)

File	Gene.count
A-BUGS-23-DC.analysis.rmd.genes.down.uniq.txt	792
A-BUGS-23-MP.analysis.rmd.genes.down.uniq.txt	802
GSE103092.analysis.rmd.genes.down.uniq.txt	48
GSE108363.bcg.discovery.analysis.Rmd.genes.down.uniq.txt	277
GSE108363.bcg.validation.analysis.Rmd.genes.down.uniq.txt	225
GSE108363.mtb.discovery.analysis.Rmd.genes.down.uniq.txt	43
GSE108363.mtb.validation.analysis.Rmd.genes.down.uniq.txt	208
GSE34151.analysis.rmd.genes.down.uniq.txt	741

Number of regulated genes:

Number of upregulated genes

up.all <- unlist(genes.up) %>% toupper
# class(up.all) # "character"
length(up.all) # 2728

[1] 2728

table.up   <- table(up.all) %>% as.data.frame(., stringsAsFactors = FALSE) %>%
  set_colnames(c("Gene", "Count")) %>%
  arrange(desc(Count))
dim(table.up) # 1955    2    2

[1] 1955    2

DT::datatable(table.up)

Number of downregulated genes

down.all <- unlist(genes.down) %>% toupper
length(down.all) # 3136

[1] 3136

table.down <- table(down.all) %>% as.data.frame(.,stringsAsFactors = FALSE) %>%
  set_colnames(c("Gene", "Count")) %>%
  arrange(desc(Count))
DT::datatable(table.down)

length(unique(table.up$Gene)) # 1955

[1] 1955

length(unique(table.down$Gene)) # 2215

[1] 2215

t.up.count <- table(table.up$Count) %>%
 as.data.frame(stringsAsFactors  = FALSE) %>%
  set_colnames(c("Count", "Freq")) 
knitr::kable(t.up.count)

Count	Freq
1	1429
2	375
3	82
4	42
5	27

# insert this as a table with all counts for both up and down regulated genes.
# extract gene symbols in 3 or more datasets and perform gene set enrichment analysis; enrichr

t.down.count <- table(table.down$Count) %>%
 as.data.frame(stringsAsFactors  = FALSE) %>%
  set_colnames(c("Count", "Freq")) 
knitr::kable(t.down.count)

Count	Freq
1	1561
2	493
3	93
4	40
5	18
6	10

Creating list of up-regulated genes with a count >= 3

genes.up.min.3 <- table.up %>%
  filter(Count >= 3) %>%
  magrittr::use_series("Gene") %>% 
  unique %>% 
  sort
# genes.up.min.3
length(genes.up.min.3) # 151

[1] 151

head(genes.up.min.3) #"ABCC3"  "ADORA3" "AIF1"   "ALDH2"  "ANG"    "APOC1" 

[1] "ABCC3"  "ADORA3" "AIF1"   "ALDH2"  "ANG"    "APOC1" 

Creating list of down-regulated genes with a count >= 3

genes.down.min.3 <- table.down %>%
  filter(Count >=3) %>%
  # .$Gene %>% 
  magrittr::use_series("Gene") %>% 
  unique %>% 
  sort 
head(genes.down.min.3) #  "ACSL1"   "ACSL5"   "ADM"     "ADORA2A" "ADTRP"   "AK4"    

[1] "ACSL1"   "ACSL5"   "ADM"     "ADORA2A" "ADTRP"   "AK4"    

length(genes.down.min.3) # 161

[1] 161

Save results

out.file.up <- paste0("data.out/", this.script, ".up.min3.txt")
cat(genes.up.min.3, file = out.file.up, sep = "\n" )
cat("File saved:", out.file.up)

File saved: data.out/counting.genes.all.DK.Rmd.up.min3.txt

out.file.down <- paste0("data.out/", this.script, ".down.min3.txt")
cat(genes.down.min.3, file = out.file.down, sep = "\n" )
cat("File saved:", out.file.down)

File saved: data.out/counting.genes.all.DK.Rmd.down.min3.txt