Analyzing gene expression data using iDEP
iDEP 0.69 http://ge-lab.org/idep/, by Steven Xijin.Ge@sdstate.edu
Mon Mar 12 20:23:55 2018
1. Read data
setwd('C:/Users/Xijin.Ge/Downloads') # Needs to be changed
source('iDEP_core_functions.R') ##
## Installing lots of R packages, this may take several hours ...
##
## Each of these packages took years to develop.
## So be patient while we steal a ton of code.
##
##
## Note 1: Sometimes dependencies needs to be installed manually.
##
## Note 2: If you are using an older version of R, and having trouble with package installation.
##
## Sometimes, it is easier to install the new version of R and remove all old packages, and start fresh.
## inputFile <- 'Downloaded_Converted_Data.csv' # Expression matrix
sampleInfoFile <- 'Downloaded_sampleInfoFile.csv' # Experiment design file
geneInfoFile <- 'Mouse__mmusculus_gene_ensembl_GeneInfo.csv' #Gene symbols, location etc.
geneSetFile <- 'Mouse__mmusculus_gene_ensembl.db' # pathway database in SQL; can be GMT format
STRING10_speciesFile <- 'https://raw.githubusercontent.com/iDEP-SDSU/idep/master/shinyapps/idep/STRING10_species.csv'
input_missingValue <- 'geneMedian' #Missing values imputation method
input_dataFileFormat <- 1 #1- read counts, 2 FKPM/RPKM or DNA microarray
input_minCounts <- 0.5 #Min counts
input_NminSamples <- 1 #Minimum number of samples
input_countsLogStart <- 4 #Pseudo count for log CPM
input_CountsTransform <- 1 #Methods for data transformation of counts. 1-VST, 2-rlog, 3-EdgeR readData.out <- readData(inputFile)
head(readData.out$data) ## p53_mock_1 p53_mock_2 p53_mock_3 p53_mock_4 p53_IR_1
## ENSMUSG00000075014 5.986190 10.245204 11.622005 12.445437 19.58986
## ENSMUSG00000076617 18.174331 18.020764 17.901045 16.892389 18.06544
## ENSMUSG00000015656 15.349228 15.086274 15.392536 14.817324 15.95674
## ENSMUSG00000052305 13.766147 11.436670 10.523028 15.927563 10.57175
## ENSMUSG00000069917 13.511114 11.202619 10.240337 15.915347 10.33765
## ENSMUSG00000075015 4.023733 7.791965 9.161421 9.809036 16.98655
## p53_IR_2 p53_IR_3 p53_IR_4 null_mock_1 null_mock_2
## ENSMUSG00000075014 21.07925 8.279954 17.40506 13.23741 12.91427
## ENSMUSG00000076617 18.52019 16.933927 17.68636 17.18909 17.39665
## ENSMUSG00000015656 15.75896 16.233039 16.04712 15.03582 14.98886
## ENSMUSG00000052305 11.22625 10.125088 10.20851 11.81668 17.59129
## ENSMUSG00000069917 11.20676 10.280270 10.32910 11.51218 17.41403
## ENSMUSG00000075015 18.47457 5.767006 14.96939 10.59269 10.25800
## null_IR_1 null_IR_2
## ENSMUSG00000075014 6.949322 14.57512
## ENSMUSG00000076617 17.663021 17.21374
## ENSMUSG00000015656 18.603855 17.00075
## ENSMUSG00000052305 10.311098 15.82093
## ENSMUSG00000069917 10.236792 15.73450
## ENSMUSG00000075015 4.742948 12.00540 readSampleInfo.out <- readSampleInfo(sampleInfoFile)
readSampleInfo.out ## p53 Treatment
## p53_mock_1 "wt" "mock"
## p53_mock_2 "wt" "mock"
## p53_mock_3 "wt" "mock"
## p53_mock_4 "wt" "mock"
## p53_IR_1 "wt" "IR"
## p53_IR_2 "wt" "IR"
## p53_IR_3 "wt" "IR"
## p53_IR_4 "wt" "IR"
## null_mock_1 "null" "mock"
## null_mock_2 "null" "mock"
## null_IR_1 "null" "IR"
## null_IR_2 "null" "IR" input_selectOrg ="NEW"
input_selectGO <- 'GOBP' #Gene set category
input_noIDConversion = TRUE
allGeneInfo.out <- geneInfo(geneInfoFile)
converted.out = NULL
convertedData.out <- convertedData()
nGenesFilter() ## [1] "6882 genes in 12 samples. 6882 genes passed filter (see above). Original gene IDs used." convertedCounts.out <- convertedCounts() # converted counts 2. Pre-process
parDefault = par()
par(mar=c(12,4,2,2))
# barplot of total read counts
barplot( colSums(readData.out$rawCounts)/1e6,
col="green",las=3, main="Total read counts (millions)") 
### Box plot
x = readData.out$data
boxplot(x, las = 2,
ylab='Transformed expression levels',
main='Distribution of transformed data') 
#### Density plot
par(parDefault) ## Warning in par(parDefault): graphical parameter "cin" cannot be set## Warning in par(parDefault): graphical parameter "cra" cannot be set## Warning in par(parDefault): graphical parameter "csi" cannot be set## Warning in par(parDefault): graphical parameter "cxy" cannot be set## Warning in par(parDefault): graphical parameter "din" cannot be set## Warning in par(parDefault): graphical parameter "page" cannot be set densityPlot() 
### Scatter plot of the first two samples
plot(x[,1:2],xlab=colnames(x)[1],ylab=colnames(x)[2],
main='Scatter plot of first two samples') 
####plot gene or gene family
input_selectOrg ="BestMatch"
input_geneSearch <- 'HOXA' #Gene ID for searching
genePlot() 
input_useSD <- 'FALSE' #Use standard deviation instead of standard error in error bar?
geneBarPlotError() 
## 3. Heatmap
# hierarchical clustering tree
x <- readData.out$data
maxGene <- apply(x,1,max)
# remove bottom 25% lowly expressed genes, which inflate the PPC
x <- x[which(maxGene > quantile(maxGene)[1] ) ,]
plot(as.dendrogram(hclust2( dist2(t(x)))), ylab="1 - Pearson C.C.", type = "rectangle") 
#Correlation matrix
input_labelPCC <- TRUE #Show correlation coefficient?
correlationMatrix() 
# Parameters for heatmap
input_nGenes <- 1000 #Top genes for heatmap
input_geneCentering <- TRUE #centering genes ?
input_sampleCentering <- FALSE #Center by sample?
input_geneNormalize <- FALSE #Normalize by gene?
input_sampleNormalize <- FALSE #Normalize by sample?
input_noSampleClustering <- FALSE #Use original sample order
input_heatmapCutoff <- 4 #Remove outliers beyond number of SDs
input_distFunctions <- 1 #which distant funciton to use
input_hclustFunctions <- 1 #Linkage type
input_heatColors1 <- 1 #Colors
input_selectFactorsHeatmap <- 'Sample_Name' #Sample coloring factors png('heatmap.png', width = 10, height = 15, units = 'in', res = 300)
staticHeatmap()
dev.off() ## png
## 2[heatmap] (heatmap.png)
heatmapPlotly() # interactive heatmap using Plotly input_nGenesKNN <- 2000 #Number of genes fro k-Means
input_nClusters <- 4 #Number of clusters
maxGeneClustering = 12000
input_kmeansNormalization <- 'geneMean' #Normalization
input_KmeansReRun <- 0 #Random seed
distributionSD() #Distribution of standard deviations 
KmeansNclusters() #Number of clusters 
Kmeans.out = Kmeans() #Running K-means
KmeansHeatmap() #Heatmap for k-Means 
#Read gene sets for enrichment analysis
sqlite <- dbDriver('SQLite')
input_selectGO3 <- 'GOBP' #Gene set category
input_minSetSize <- 15 #Min gene set size
input_maxSetSize <- 2000 #Max gene set size
GeneSets.out <-readGeneSets( geneSetFile,
convertedData.out, input_selectGO3,input_selectOrg,
c(input_minSetSize, input_maxSetSize) )
# Alternatively, users can use their own GMT files by
#GeneSets.out <- readGMTRobust('somefile.GMT')
KmeansGO() #Enrichment analysis for k-Means clusters## Cluster
## RNA processing A
## Ribonucleoprotein complex biogenesis
## NcRNA metabolic process
## Ribosome biogenesis
## NcRNA processing
## Organonitrogen compound metabolic process
## Organonitrogen compound biosynthetic process
## RRNA processing
## RRNA metabolic process
## Translation
## Cellular response to stress B
## Apoptotic process
## Cell death
## Programmed cell death
## Apoptotic signaling pathway
## Positive regulation of cell death
## Cellular response to DNA damage stimulus
## Intrinsic apoptotic signaling pathway
## Positive regulation of programmed cell death
## Intrinsic apoptotic signaling pathway in response to DNA damage
## Immune system process C
## Positive regulation of response to stimulus
## Immune response
## Establishment of protein localization
## Macromolecular complex assembly
## Regulation of immune system process
## Protein complex assembly
## Protein transport
## Protein complex biogenesis
## Protein complex subunit organization
## Immune system process1 D
## Positive regulation of response to stimulus1
## Positive regulation of molecular function
## Protein phosphorylation
## Defense response
## Regulation of intracellular signal transduction
## Cell activation
## Regulation of transport
## Response to external stimulus
## Vesicle-mediated transport
## adj.Pval
## RNA processing 2.685004e-144
## Ribonucleoprotein complex biogenesis 1.992995e-140
## NcRNA metabolic process 2.289446e-134
## Ribosome biogenesis 3.496906e-117
## NcRNA processing 6.343278e-115
## Organonitrogen compound metabolic process 3.482536e-107
## Organonitrogen compound biosynthetic process 1.711515e-101
## RRNA processing 6.861087e-93
## RRNA metabolic process 2.888202e-92
## Translation 5.045222e-75
## Cellular response to stress 2.075602e-25
## Apoptotic process 4.559668e-23
## Cell death 4.559668e-23
## Programmed cell death 5.845102e-23
## Apoptotic signaling pathway 5.845102e-23
## Positive regulation of cell death 1.232970e-20
## Cellular response to DNA damage stimulus 1.439213e-20
## Intrinsic apoptotic signaling pathway 3.577557e-19
## Positive regulation of programmed cell death 1.011040e-18
## Intrinsic apoptotic signaling pathway in response to DNA damage 1.179594e-18
## Immune system process 2.394955e-63
## Positive regulation of response to stimulus 8.229461e-45
## Immune response 1.180337e-44
## Establishment of protein localization 1.514473e-42
## Macromolecular complex assembly 7.228082e-42
## Regulation of immune system process 2.287631e-40
## Protein complex assembly 2.287631e-40
## Protein transport 2.287631e-40
## Protein complex biogenesis 2.287631e-40
## Protein complex subunit organization 2.287631e-40
## Immune system process1 9.615041e-94
## Positive regulation of response to stimulus1 6.158140e-55
## Positive regulation of molecular function 5.378810e-52
## Protein phosphorylation 6.500699e-52
## Defense response 1.857527e-51
## Regulation of intracellular signal transduction 7.246982e-51
## Cell activation 1.023829e-49
## Regulation of transport 2.242061e-47
## Response to external stimulus 1.004958e-46
## Vesicle-mediated transport 5.131263e-46
## Genes
## RNA processing 121
## Ribonucleoprotein complex biogenesis 106
## NcRNA metabolic process 103
## Ribosome biogenesis 84
## NcRNA processing 84
## Organonitrogen compound metabolic process 144
## Organonitrogen compound biosynthetic process 119
## RRNA processing 65
## RRNA metabolic process 65
## Translation 74
## Cellular response to stress 36
## Apoptotic process 36
## Cell death 37
## Programmed cell death 36
## Apoptotic signaling pathway 25
## Positive regulation of cell death 24
## Cellular response to DNA damage stimulus 23
## Intrinsic apoptotic signaling pathway 18
## Positive regulation of programmed cell death 22
## Intrinsic apoptotic signaling pathway in response to DNA damage 14
## Immune system process 115
## Positive regulation of response to stimulus 86
## Immune response 76
## Establishment of protein localization 81
## Macromolecular complex assembly 76
## Regulation of immune system process 69
## Protein complex assembly 69
## Protein transport 75
## Protein complex biogenesis 69
## Protein complex subunit organization 72
## Immune system process1 147
## Positive regulation of response to stimulus1 99
## Positive regulation of molecular function 93
## Protein phosphorylation 93
## Defense response 84
## Regulation of intracellular signal transduction 86
## Cell activation 72
## Regulation of transport 88
## Response to external stimulus 95
## Vesicle-mediated transport 73
## Pathways
## RNA processing RNA processing
## Ribonucleoprotein complex biogenesis Ribonucleoprotein complex biogenesis
## NcRNA metabolic process NcRNA metabolic process
## Ribosome biogenesis Ribosome biogenesis
## NcRNA processing NcRNA processing
## Organonitrogen compound metabolic process Organonitrogen compound metabolic process
## Organonitrogen compound biosynthetic process Organonitrogen compound biosynthetic process
## RRNA processing RRNA processing
## RRNA metabolic process RRNA metabolic process
## Translation Translation
## Cellular response to stress Cellular response to stress
## Apoptotic process Apoptotic process
## Cell death Cell death
## Programmed cell death Programmed cell death
## Apoptotic signaling pathway Apoptotic signaling pathway
## Positive regulation of cell death Positive regulation of cell death
## Cellular response to DNA damage stimulus Cellular response to DNA damage stimulus
## Intrinsic apoptotic signaling pathway Intrinsic apoptotic signaling pathway
## Positive regulation of programmed cell death Positive regulation of programmed cell death
## Intrinsic apoptotic signaling pathway in response to DNA damage Intrinsic apoptotic signaling pathway in response to DNA damage
## Immune system process Immune system process
## Positive regulation of response to stimulus Positive regulation of response to stimulus
## Immune response Immune response
## Establishment of protein localization Establishment of protein localization
## Macromolecular complex assembly Macromolecular complex assembly
## Regulation of immune system process Regulation of immune system process
## Protein complex assembly Protein complex assembly
## Protein transport Protein transport
## Protein complex biogenesis Protein complex biogenesis
## Protein complex subunit organization Protein complex subunit organization
## Immune system process1 Immune system process
## Positive regulation of response to stimulus1 Positive regulation of response to stimulus
## Positive regulation of molecular function Positive regulation of molecular function
## Protein phosphorylation Protein phosphorylation
## Defense response Defense response
## Regulation of intracellular signal transduction Regulation of intracellular signal transduction
## Cell activation Cell activation
## Regulation of transport Regulation of transport
## Response to external stimulus Response to external stimulus
## Vesicle-mediated transport Vesicle-mediated transport input_seedTSNE <- 0 #Random seed for t-SNE
input_colorGenes <- TRUE #Color genes in t-SNE plot?
tSNEgenePlot() #Plot genes using t-SNE 
## 5. PCA and beyond
input_selectFactors <- 'Sample_Name' #Factor coded by color
input_selectFactors2 <- 'Sample_Name' #Factor coded by shape
input_tsneSeed2 <- 0 #Random seed for t-SNE
#PCA, MDS and t-SNE plots
PCAplot() 
MDSplot() 
tSNEplot() 
#Read gene sets for pathway analysis using PGSEA on principal components
input_selectGO6 <- 'GOBP'
GeneSets.out <-readGeneSets( geneSetFile,
convertedData.out, input_selectGO6,input_selectOrg,
c(input_minSetSize, input_maxSetSize) )
PCApathway() # Run PGSEA analysis 
cat( PCA2factor() ) #The correlation between PCs with factors ##
## Correlation between Principal Components (PCs) with factors
## PC1 is correlated with Treatment (p=6.56e-05).6. DEG1
input_CountsDEGMethod <- 3 #DESeq2= 3,limma-voom=2,limma-trend=1
input_limmaPval <- 0.1 #FDR cutoff
input_limmaFC <- 2 #Fold-change cutoff
input_selectModelComprions <- NULL #Selected comparisons
input_selectFactorsModel <- NULL #Selected comparisons
input_selectInteractions <- NULL #Selected comparisons
input_selectBlockFactorsModel <- NULL #Selected comparisons
factorReferenceLevels.out <- NULL
limma.out <- limma()
DEG.data.out <- DEG.data()
limma.out$comparisons ## [1] "null_mock-null_IR" "p53_IR-null_IR" "p53_mock-null_IR"
## [4] "p53_IR-null_mock" "p53_mock-null_mock" "p53_mock-p53_IR" input_selectComparisonsVenn = limma.out$comparisons[1:3] # use first three comparisons
input_UpDownRegulated <- FALSE #Split up and down regulated genes
vennPlot() # Venn diagram 
sigGeneStats() # number of DEGs as figure 
sigGeneStatsTable() # number of DEGs as table ## Comparisons Up Down
## null_mock-null_IR null_mock-null_IR 702 1532
## p53_IR-null_IR p53_IR-null_IR 1118 577
## p53_mock-null_IR p53_mock-null_IR 910 1687
## p53_IR-null_mock p53_IR-null_mock 1741 540
## p53_mock-null_mock p53_mock-null_mock 20 16
## p53_mock-p53_IR p53_mock-p53_IR 487 18567. DEG2
##########################
# 7. DEG2
##########################
input_selectContrast <- 'null_mock-null_IR' #Selected comparisons
selectedHeatmap.data.out <- selectedHeatmap.data()
selectedHeatmap() # heatmap for DEGs in selected comparison
# Save gene lists and data into files
write.csv( selectedHeatmap.data()$genes, 'heatmap.data.csv')
write.csv(DEG.data(),'DEG.data.csv' )
write(AllGeneListsGMT() ,'AllGeneListsGMT.gmt') input_selectGO2 <- 'GOBP' #Gene set category
geneListData.out <- geneListData()
volcanoPlot() 
scatterPlot() 
MAplot() 
geneListGOTable.out <- geneListGOTable()
# Read pathway data again
GeneSets.out <-readGeneSets( geneSetFile,
convertedData.out, input_selectGO2,input_selectOrg,
c(input_minSetSize, input_maxSetSize) )
input_removeRedudantSets <- TRUE #Remove highly redundant gene sets?
geneListGO() ## Direction adj.Pval
## Organonitrogen compound metabolic process Down regulated 6.1e-157
## RNA processing 9.4e-155
## NcRNA metabolic process 4.1e-151
## Ribonucleoprotein complex biogenesis 7.9e-143
## Organonitrogen compound biosynthetic process 3.3e-136
## NcRNA processing 2.7e-120
## Ribosome biogenesis 6.8e-120
## Small molecule metabolic process 4.3e-114
## RRNA processing 6.5e-90
## RRNA metabolic process 7.6e-89
## Cellular response to organic substance Up regulated 3.1e-34
## Immune system process 3.0e-32
## Response to oxygen-containing compound 1.7e-31
## Response to endogenous stimulus 2.2e-30
## Regulation of transport 1.7e-26
## Response to nitrogen compound 1.9e-26
## Lipid metabolic process 2.8e-26
## Cellular response to oxygen-containing compound 3.6e-24
## Response to organonitrogen compound 5.5e-24
## Transmembrane transport 5.5e-24
## nGenes
## Organonitrogen compound metabolic process 251
## RNA processing 159
## NcRNA metabolic process 135
## Ribonucleoprotein complex biogenesis 131
## Organonitrogen compound biosynthetic process 192
## NcRNA processing 104
## Ribosome biogenesis 102
## Small molecule metabolic process 202
## RRNA processing 76
## RRNA metabolic process 76
## Cellular response to organic substance 86
## Immune system process 89
## Response to oxygen-containing compound 72
## Response to endogenous stimulus 70
## Regulation of transport 69
## Response to nitrogen compound 52
## Lipid metabolic process 57
## Cellular response to oxygen-containing compound 51
## Response to organonitrogen compound 47
## Transmembrane transport 57
## Pathways
## Organonitrogen compound metabolic process Organonitrogen compound metabolic process
## RNA processing RNA processing
## NcRNA metabolic process NcRNA metabolic process
## Ribonucleoprotein complex biogenesis Ribonucleoprotein complex biogenesis
## Organonitrogen compound biosynthetic process Organonitrogen compound biosynthetic process
## NcRNA processing NcRNA processing
## Ribosome biogenesis Ribosome biogenesis
## Small molecule metabolic process Small molecule metabolic process
## RRNA processing RRNA processing
## RRNA metabolic process RRNA metabolic process
## Cellular response to organic substance Cellular response to organic substance
## Immune system process Immune system process
## Response to oxygen-containing compound Response to oxygen-containing compound
## Response to endogenous stimulus Response to endogenous stimulus
## Regulation of transport Regulation of transport
## Response to nitrogen compound Response to nitrogen compound
## Lipid metabolic process Lipid metabolic process
## Cellular response to oxygen-containing compound Cellular response to oxygen-containing compound
## Response to organonitrogen compound Response to organonitrogen compound
## Transmembrane transport Transmembrane transport # STRING-db API access
STRING10_species = read.csv(STRING10_speciesFile)
ix = grep('Mus musculus', STRING10_species$official_name )
findTaxonomyID.out <- STRING10_species[ix,1] # find taxonomyID
findTaxonomyID.out ## [1] 10090 # users can also skip the above and assign NCBI taxonomy id directly by
# findTaxonomyID.out = 10090 # mouse 10090, human 9606 etc. STRINGdb_geneList.out <- STRINGdb_geneList() #convert gene lists## Warning: we couldn't map to STRING 9% of your identifiers input_STRINGdbGO <- 'Process' #'Process', 'Component', 'Function', 'KEGG', 'Pfam', 'InterPro'
stringDB_GO_enrichmentData() ## Direction adj.Pval nGenes
## 1 Up regulated 7.3e-08 49
## 2 5.1e-07 41
## 3 8.4e-07 65
## 4 8.4e-07 46
## 5 8.4e-07 53
## 6 9.5e-07 37
## 7 1.0e-06 49
## 8 1.3e-06 60
## 9 2.3e-06 60
## 10 9.5e-06 40
## 11 9.5e-06 36
## 12 9.5e-06 35
## 13 5.0e-05 26
## 14 5.4e-05 46
## 15 6.6e-05 51
## 16 Down regulated 7.5e-62 107
## 17 2.8e-55 106
## 18 2.1e-52 81
## 19 1.6e-51 136
## 20 4.2e-50 146
## 21 2.3e-46 85
## 22 1.3e-43 181
## 23 9.8e-41 61
## 24 7.4e-40 61
## 25 8.8e-32 84
## 26 1.9e-29 72
## 27 1.9e-28 130
## 28 1.1e-27 72
## 29 2.2e-27 91
## 30 8.4e-27 156
## Pathways
## 1 lipid metabolic process
## 2 response to nitrogen compound
## 3 regulation of transport
## 4 transmembrane transport
## 5 response to oxygen-containing compound
## 6 response to organonitrogen compound
## 7 ion transport
## 8 positive regulation of molecular function
## 9 immune system process
## 10 positive regulation of transport
## 11 positive regulation of hydrolase activity
## 12 cellular response to oxygen-containing compound
## 13 positive regulation of GTPase activity
## 14 regulation of hydrolase activity
## 15 intracellular signal transduction
## 16 ribonucleoprotein complex biogenesis
## 17 ncRNA metabolic process
## 18 ribosome biogenesis
## 19 RNA processing
## 20 organonitrogen compound biosynthetic process
## 21 ncRNA processing
## 22 organonitrogen compound metabolic process
## 23 rRNA processing
## 24 rRNA metabolic process
## 25 amide biosynthetic process
## 26 translation
## 27 single-organism biosynthetic process
## 28 peptide biosynthetic process
## 29 cellular amide metabolic process
## 30 small molecule metabolic process # PPI network retrieval and analysis
input_nGenesPPI <- 100 #Number of top genes for PPI retrieval and analysis
stringDB_network1(1) #Show PPI network 
write(stringDB_network_link(), 'PPI_results.html') # write results to html file ## Warning: we couldn't map to STRING 9% of your identifiers browseURL('PPI_results.html') # open in browser 8. Pathway analysis
##########################
# 8. Pathway analysis
##########################
input_selectContrast1 <- 'null_mock-null_IR' #select Comparison
#input_selectContrast1 = limma.out$comparisons[3] # manually set
input_selectGO <- 'GOBP' #Gene set category
#input_selectGO='custom' # if custom gmt file
input_minSetSize <- 15 #Min size for gene set
input_maxSetSize <- 2000 #Max size for gene set
# Read pathway data again
GeneSets.out <-readGeneSets( geneSetFile,
convertedData.out, input_selectGO,input_selectOrg,
c(input_minSetSize, input_maxSetSize) )
input_pathwayPvalCutoff <- 0.2 #FDR cutoff
input_nPathwayShow <- 30 #Top pathways to show
input_absoluteFold <- FALSE #Use absolute values of fold-change?
input_GenePvalCutoff <- 1 #FDR to remove genes
input_pathwayMethod = 1 # 1 GAGE
gagePathwayData.out <- gagePathwayData() # pathway analysis using GAGE
gagePathwayData.out ## Direction
## NcRNA metabolic process Down
## Ribosome biogenesis
## NcRNA processing
## Ribonucleoprotein complex biogenesis
## RRNA processing
## RRNA metabolic process
## RNA processing
## RNA modification
## TRNA metabolic process
## Ribosomal large subunit biogenesis
## TRNA processing
## TRNA modification
## Maturation of SSU-rRNA
## Protein folding
## Translation
## Organonitrogen compound biosynthetic process
## Peptide biosynthetic process
## Ribonucleoprotein complex subunit organization
## Ribonucleoprotein complex assembly
## Amide biosynthetic process
## Cellular amino acid metabolic process
## RNA localization
## Ribosomal small subunit biogenesis
## MRNA metabolic process
## Translational initiation
## DNA metabolic process
## Nucleoside metabolic process
## Ribonucleoside metabolic process
## Peptide metabolic process
## MRNA processing
## GAGE analysis: null_mock vs null_IR
## NcRNA metabolic process NcRNA metabolic process
## Ribosome biogenesis Ribosome biogenesis
## NcRNA processing NcRNA processing
## Ribonucleoprotein complex biogenesis Ribonucleoprotein complex biogenesis
## RRNA processing RRNA processing
## RRNA metabolic process RRNA metabolic process
## RNA processing RNA processing
## RNA modification RNA modification
## TRNA metabolic process TRNA metabolic process
## Ribosomal large subunit biogenesis Ribosomal large subunit biogenesis
## TRNA processing TRNA processing
## TRNA modification TRNA modification
## Maturation of SSU-rRNA Maturation of SSU-rRNA
## Protein folding Protein folding
## Translation Translation
## Organonitrogen compound biosynthetic process Organonitrogen compound biosynthetic process
## Peptide biosynthetic process Peptide biosynthetic process
## Ribonucleoprotein complex subunit organization Ribonucleoprotein complex subunit organization
## Ribonucleoprotein complex assembly Ribonucleoprotein complex assembly
## Amide biosynthetic process Amide biosynthetic process
## Cellular amino acid metabolic process Cellular amino acid metabolic process
## RNA localization RNA localization
## Ribosomal small subunit biogenesis Ribosomal small subunit biogenesis
## MRNA metabolic process MRNA metabolic process
## Translational initiation Translational initiation
## DNA metabolic process DNA metabolic process
## Nucleoside metabolic process Nucleoside metabolic process
## Ribonucleoside metabolic process Ribonucleoside metabolic process
## Peptide metabolic process Peptide metabolic process
## MRNA processing MRNA processing
## statistic Genes adj.Pval
## NcRNA metabolic process -11.8689 154 3.7e-23
## Ribosome biogenesis -10.8929 110 1.1e-18
## NcRNA processing -10.5686 114 2.0e-18
## Ribonucleoprotein complex biogenesis -10.562 154 2.9e-19
## RRNA processing -9.7918 81 1.3e-14
## RRNA metabolic process -9.6836 82 1.6e-14
## RNA processing -9.4752 216 5.3e-17
## RNA modification -7.8231 42 3.5e-08
## TRNA metabolic process -7.5376 51 1.0e-08
## Ribosomal large subunit biogenesis -6.7767 31 5.7e-06
## TRNA processing -6.1338 30 1.2e-05
## TRNA modification -5.7857 20 4.8e-04
## Maturation of SSU-rRNA -5.6703 21 4.8e-04
## Protein folding -5.6113 72 1.1e-05
## Translation -5.6082 172 5.6e-06
## Organonitrogen compound biosynthetic process -5.4954 397 5.7e-06
## Peptide biosynthetic process -5.1789 178 3.3e-05
## Ribonucleoprotein complex subunit organization -5.1249 70 7.8e-05
## Ribonucleoprotein complex assembly -4.9952 66 1.3e-04
## Amide biosynthetic process -4.9881 203 7.4e-05
## Cellular amino acid metabolic process -4.9853 70 1.3e-04
## RNA localization -4.9661 67 1.4e-04
## Ribosomal small subunit biogenesis -4.8412 27 7.3e-04
## MRNA metabolic process -4.6528 120 3.4e-04
## Translational initiation -4.6098 33 1.3e-03
## DNA metabolic process -4.4818 288 4.9e-04
## Nucleoside metabolic process -4.4411 132 6.9e-04
## Ribonucleoside metabolic process -4.2843 127 1.2e-03
## Peptide metabolic process -4.247 227 1.2e-03
## MRNA processing -4.2087 94 1.7e-03 pathwayListData.out = pathwayListData()
enrichmentPlot(pathwayListData.out, 25 ) 
enrichmentNetwork(pathwayListData.out ) 
enrichmentNetworkPlotly(pathwayListData.out) 
input_pathwayMethod = 3 # 1 fgsea
fgseaPathwayData.out <- fgseaPathwayData() #Pathway analysis using fgsea
fgseaPathwayData.out ## Direction GSEA analysis: null_mock vs null_IR NES
## 1247 Down NcRNA metabolic process -2.636
## 1234 NcRNA processing -2.6069
## 1339 Ribosome biogenesis -2.5705
## 923 Ribonucleoprotein complex biogenesis -2.5522
## 289 RRNA processing -2.511
## 804 RRNA metabolic process -2.5083
## 294 TRNA metabolic process -2.438
## 525 TRNA processing -2.4074
## 292 RNA processing -2.3922
## 619 RNA modification -2.3841
## 1343 Ribosomal large subunit biogenesis -2.2866
## 295 TRNA modification -2.2736
## 322 Cellular amino acid metabolic process -2.2331
## 304 Protein folding -2.2127
## 2352 RNA phosphodiester bond hydrolysis endonucleolytic -2.1758
## 992 Maturation of SSU-rRNA -2.1705
## 300 Translation -2.1502
## 2262 Ribonucleoprotein complex subunit organization -2.1476
## 925 Ribonucleoprotein complex assembly -2.1474
## 32 Telomere maintenance -2.1332
## 1099 Telomere organization -2.1332
## 301 Translational initiation -2.106
## 303 TRNA aminoacylation for protein translation -2.0919
## 1402 Peptide biosynthetic process -2.0795
## 298 RNA localization -2.0793
## 731 Telomere maintenance via telomere lengthening -2.0743
## 1344 Ribosomal small subunit biogenesis -2.0724
## 1400 Amino acid activation -2.0694
## 1401 TRNA aminoacylation -2.0694
## 2124 Chaperone-mediated protein folding -2.0685
## Genes adj.Pval
## 1247 154 3.6e-03
## 1234 114 3.6e-03
## 1339 110 3.6e-03
## 923 154 3.6e-03
## 289 81 3.6e-03
## 804 82 3.6e-03
## 294 51 3.6e-03
## 525 30 3.6e-03
## 292 216 3.6e-03
## 619 42 3.6e-03
## 1343 31 3.6e-03
## 295 20 3.6e-03
## 322 70 3.6e-03
## 304 72 3.6e-03
## 2352 19 3.6e-03
## 992 21 3.6e-03
## 300 172 3.6e-03
## 2262 70 3.6e-03
## 925 66 3.6e-03
## 32 35 3.6e-03
## 1099 35 3.6e-03
## 301 33 3.6e-03
## 303 19 3.6e-03
## 1402 178 3.6e-03
## 298 67 3.6e-03
## 731 24 3.6e-03
## 1344 27 3.6e-03
## 1400 20 3.6e-03
## 1401 20 3.6e-03
## 2124 29 3.6e-03 pathwayListData.out = pathwayListData()
enrichmentPlot(pathwayListData.out, 25 ) 
enrichmentNetwork(pathwayListData.out ) 
enrichmentNetworkPlotly(pathwayListData.out) 
PGSEAplot() # pathway analysis using PGSEA ##
## Computing P values using ANOVA
9. Chromosome
input_selectContrast2 <- 'null_mock-null_IR' #select Comparison
#input_selectContrast2 = limma.out$comparisons[3] # manually set
input_limmaPvalViz <- 0.1 #FDR to filter genes
input_limmaFCViz <- 2 #FDR to filter genes
genomePlotly() # shows fold-changes on the genome ## Warning in eval(quote(list(...)), env): NAs introduced by coercion## Warning in genomePlotly(): NAs introduced by coercion input_nGenesBiclust <- 1000 #Top genes for biclustering
input_biclustMethod <- 'BCCC()' #Method: 'BCCC', 'QUBIC', 'runibic' ...
biclustering.out = biclustering() # run analysis
input_selectBicluster <- 1 #select a cluster
biclustHeatmap() # heatmap for selected cluster 
input_selectGO4 <- 'GOBP' #Gene set
# Read pathway data again
GeneSets.out <-readGeneSets( geneSetFile,
convertedData.out, input_selectGO4,input_selectOrg,
c(input_minSetSize, input_maxSetSize) )
geneListBclustGO() # enrichment analysis ## adj.Pval Genes
## Immune system process 3.1e-124 207
## Organonitrogen compound metabolic process 1.5e-89 157
## Positive regulation of molecular function 4.7e-86 148
## Cell death 5.4e-86 155
## Homeostatic process 2.0e-83 139
## Macromolecular complex assembly 3.2e-81 136
## Programmed cell death 3.5e-80 146
## Apoptotic process 2.9e-79 144
## Regulation of cell death 1.1e-75 132
## Cell-cell adhesion 2.2e-75 117
## Pathways
## Immune system process Immune system process
## Organonitrogen compound metabolic process Organonitrogen compound metabolic process
## Positive regulation of molecular function Positive regulation of molecular function
## Cell death Cell death
## Homeostatic process Homeostatic process
## Macromolecular complex assembly Macromolecular complex assembly
## Programmed cell death Programmed cell death
## Apoptotic process Apoptotic process
## Regulation of cell death Regulation of cell death
## Cell-cell adhesion Cell-cell adhesion11. Co-expression network
##########################
# 11. Co-expression network
##########################
input_mySoftPower <- 5 #SoftPower to cutoff
input_nGenesNetwork <- 1000 #Number of top genes
input_minModuleSize <- 20 #Module size minimum
wgcna.out = wgcna() # run WGCNA## ==========================================================================
## *
## * Package WGCNA 1.63 loaded.
## *
## * Important note: It appears that your system supports multi-threading,
## * but it is not enabled within WGCNA in R.
## * To allow multi-threading within WGCNA with all available cores, use
## *
## * allowWGCNAThreads()
## *
## * within R. Use disableWGCNAThreads() to disable threading if necessary.
## * Alternatively, set the following environment variable on your system:
## *
## * ALLOW_WGCNA_THREADS=<number_of_processors>
## *
## * for example
## *
## * ALLOW_WGCNA_THREADS=8
## *
## * To set the environment variable in linux bash shell, type
## *
## * export ALLOW_WGCNA_THREADS=8
## *
## * before running R. Other operating systems or shells will
## * have a similar command to achieve the same aim.
## *
## ==========================================================================## Warning: executing %dopar% sequentially: no parallel backend registered## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1 0.9720 2.7600 0.966 572.0 612.0 711
## 2 2 0.9490 1.2200 0.956 412.0 449.0 578
## 3 3 0.9050 0.7070 0.923 325.0 353.0 495
## 4 4 0.7480 0.4190 0.797 269.0 287.0 439
## 5 5 0.5350 0.2490 0.615 230.0 239.0 398
## 6 6 0.1820 0.1090 0.255 200.0 204.0 367
## 7 7 0.0040 0.0131 0.185 178.0 175.0 342
## 8 8 0.0965 -0.0738 0.144 160.0 152.0 321
## 9 9 0.2500 -0.1240 0.238 145.0 133.0 303
## 10 10 0.4570 -0.1930 0.431 132.0 118.0 287
## 11 12 0.5670 -0.2940 0.514 113.0 98.9 263
## 12 14 0.6340 -0.3830 0.587 97.8 84.7 243
## 13 16 0.6890 -0.4420 0.655 86.2 73.5 227
## 14 18 0.7120 -0.4930 0.673 76.9 63.8 213
## 15 20 0.7290 -0.5420 0.683 69.3 55.7 201
## TOM calculation: adjacency..
## ..will not use multithreading.
## Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done. softPower() # soft power curve 
modulePlot() # plot modules 
listWGCNA.Modules.out = listWGCNA.Modules() #modules input_selectGO5 <- 'GOBP' #Gene set
# Read pathway data again
GeneSets.out <-readGeneSets( geneSetFile,
convertedData.out, input_selectGO5,input_selectOrg,
c(input_minSetSize, input_maxSetSize) )
input_selectWGCNA.Module <- '1. turquoise (495 genes)' #Select a module
input_topGenesNetwork <- 10 #SoftPower to cutoff
input_edgeThreshold <- 0.4 #Number of top genes
moduleNetwork() # show network of top genes in selected module## softConnectivity: FYI: connecitivty of genes with less than 4 valid samples will be returned as NA.
## ..calculating connectivities..
input_removeRedudantSets <- TRUE #Remove redundant gene sets
networkModuleGO() # Enrichment analysis of selected module ## adj.Pval Genes
## Organonitrogen compound metabolic process 6.2e-68 106
## Ribonucleoprotein complex biogenesis 3.0e-62 61
## RNA processing 6.9e-62 68
## Organonitrogen compound biosynthetic process 4.0e-61 84
## NcRNA metabolic process 2.4e-56 57
## Immune system process 4.2e-53 101
## Macromolecular complex assembly 2.0e-49 81
## Ribosome biogenesis 3.0e-48 46
## Establishment of protein localization 3.4e-48 84
## Cellular amide metabolic process 4.7e-48 63
## Pathways
## Organonitrogen compound metabolic process Organonitrogen compound metabolic process
## Ribonucleoprotein complex biogenesis Ribonucleoprotein complex biogenesis
## RNA processing RNA processing
## Organonitrogen compound biosynthetic process Organonitrogen compound biosynthetic process
## NcRNA metabolic process NcRNA metabolic process
## Immune system process Immune system process
## Macromolecular complex assembly Macromolecular complex assembly
## Ribosome biogenesis Ribosome biogenesis
## Establishment of protein localization Establishment of protein localization
## Cellular amide metabolic process Cellular amide metabolic process