- 1 0. Install the package
- 2 1. Simplified GEO microarray workflow
- 3 2. Common small utilities
- 4 3. Download GEO-provided RNA-seq count data
- 5 4. Expression matrix conversion
- 6 5. Common plots
- 7 6. Survival analysis
- 8 7. Network and statistics utilities
- 9 8. Enrichment analysis
- 10 9. Updates
✍️ Author: Xiaojie Sun
🌟 Role: Full-time lecturer at “Bioinformatics Skill Tree” for eight years, co-founder of the WeChat public account “Bioinformatics Planet”, and co-author of an SCI-indexed paper with an impact factor of approximately 8. Co-editor of the second edition of Decoding Life; upcoming book: Practical Bioinformatics Data Analysis with R.
🎓 Focused on bioinformatics data analysis training for beginners and hands-on teaching.
tinyarray is an R package designed to simplify downstream GEO/TCGA analysis. It integrates data download, ID conversion, differential expression analysis, survival analysis, and commonly used bioinformatics visualization. Many of the functions in this package were originally written to meet my own daily analysis, teaching, and Q&A needs, and to solve practical pain points.
tinyarray depends on several Bioconductor packages, and
install.packages() cannot install all of them directly. You can
install those dependencies first and then install tinyarray.
if(!requireNamespace("BiocManager", quietly = TRUE))install.packages("BiocManager")
pkg_list <- c("limma", "GEOquery", "Biobase", "ComplexHeatmap", "AnnotationDbi", "clusterProfiler", "org.Rn.eg.db", "org.Mm.eg.db", "org.Hs.eg.db")
BiocManager::install(pkg_list,update = F,ask = F)if(!require(tinyarray))install.packages("tinyarray")Here we use GSE7305 as an example and walk through downloading,
grouping, annotating, and differential expression analysis in order.
library(tinyarray)
library(stringr)# Download data; use the local cached example when available so the tutorial
# can still render in offline environments.
if (file.exists("GSE7305_eSet.rds")) {
geo_cached <- readRDS("GSE7305_eSet.rds")
if (is.list(geo_cached) && length(geo_cached) == 1L && inherits(geo_cached[[1]], "ExpressionSet")) {
geo_eset <- geo_cached[[1]]
geo1 <- list(
exp = Biobase::exprs(geo_eset),
pd = Biobase::pData(geo_eset),
gpl = Biobase::annotation(geo_eset)
)
} else {
geo1 <- geo_cached
}
} else {
geo1 <- geo_download("GSE7305", destdir = tempdir())
}## you can also use getGEO from GEOquery, by
## getGEO("GSE7305", destdir=".", AnnotGPL = F, getGPL = F)
## 54675 probes, 20 samples from 5.020951 to 22011.934
if (is.null(geo1$group_candidates)) {
geo1$group_candidates <- auto_geo_group(geo1$pd)
}geo_download() returns a list with the following elements:
exp: expression matrixpd: phenotype datagpl: platform accessiongroup_candidates: automatic grouping candidates
This is also the most common entry point for downstream differential analysis, annotation conversion, and automatic grouping.
# Extract expression matrix
exp1 = geo1$exp
# Decide whether to log-transform based on the data range; values between 0 and 20 are usually already logged
if(max(exp1) >50)exp1 <- log2(exp1 + 1)
# Extract clinical information
pd1 <- geo1$pd
# Get grouping information; geo1 contains a group_candidates subset that stores candidate groups
geo1$group_candidates## $`title_Disease-Normal`
## [1] Disease Disease Disease Disease Disease Disease Disease Disease Disease
## [10] Disease Normal Normal Normal Normal Normal Normal Normal Normal
## [19] Normal Normal
## Levels: Disease Normal
##
## $`title_Ovary-Normal`
## [1] Ovary Ovary Ovary Ovary Ovary Ovary Ovary Ovary Ovary Ovary
## [11] Normal Normal Normal Normal Normal Normal Normal Normal Normal Normal
## Levels: Ovary Normal
##
## $`source_name_ch1_G-E`
## [1] E E E E E G G G G G E E E E E G G G G G
## Levels: G E
##
## $`source_name_ch1_M-E`
## [1] M M M M M E E E E E M M M M M E E E E E
## Levels: M E
##
## $`source_name_ch1_N-E`
## [1] E E E E E E E E E E N N N N N N N N N N
## Levels: N E
##
## $`description_Luteal-Follicular`
## [1] Follicular Follicular Luteal Luteal Follicular Follicular
## [7] Follicular Follicular Follicular Follicular Follicular Follicular
## [13] Luteal Luteal Follicular Follicular Follicular Follicular
## [19] Follicular Follicular
## Levels: Luteal Follicular
# For example, select the first candidate as the grouping code
Group1 = geo1$group_candidates[[1]]
# If there is no suitable candidate group, create one manually
k = str_detect(pd1$title, "Disease");table(k)## k
## FALSE TRUE
## 10 10
Group1 = ifelse(k, "Disease", "Normal")
Group1 = factor(Group1, levels = c("Normal", "Disease"))
table(Group1)## Group1
## Normal Disease
## 10 10
# Get probe annotation; for some platforms this cannot be obtained automatically, and the GPL page needs to be downloaded and parsed manually
ids1 <- get_ids(geo1$gpl)##
## The following code was used to obtain the annotation table:
## library(hgu133plus2.db); ids <- toTable(hgu133plus2SYMBOL)
# Check whether the probe annotation is correct; if most values below are TRUE, the annotation is usually correct
table(ids1$probe_id %in% rownames(exp1))##
## TRUE
## 43086
head(ids1)## probe_id symbol
## 1 1007_s_at DDR1
## 2 1053_at RFC2
## 3 117_at HSPA6
## 4 121_at PAX8
## 5 1255_g_at GUCA1A
## 6 1294_at UBA7
get_ids() returns a cleaned two-column table:
- first column:
probe_id - second column:
symbol
It automatically does several things:
- Prefer an available Bioconductor annotation package;
- If none is available, try
AnnoProbe::idmap(); - Remove probes whose symbol is
NAor an empty string""; - Drop extra columns and keep only
probe_idandsymbol; - Convert
probe_idto character; - Reset row names.
If get_ids() gives a warning like:
No annotation table is available for GPL5555. Please extract it manually
from the GPL platform table. You can use get_gpl_txt("GPL5555") to get
the download link for the GPL platform table.
it means that neither a Bioconductor package nor AnnoProbe::idmap()
could provide a valid annotation, and you need to extract it manually
from the GPL platform table.
You can use get_gpl_txt() to get the platform-table download link, for
example:
get_gpl_txt(geo1$gpl)## https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570&targ=self&form=text&view=data
get_gpl_txt("GPL6244")## https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL6244&targ=self&form=text&view=data
Download and read that file, then extract the probe_id and symbol
columns and rename them to "probe_id" and "symbol", and they can be
used in downstream analysis.
# Differential expression analysis
deg1 <- get_deg(exp1, Group1, ids1, logFC_cutoff = 1, entriz = FALSE)
head(deg1)## logFC AveExpr t P.Value adj.P.Val B probe_id
## 1 6.270309 8.436140 45.39552 9.106509e-24 2.489492e-19 41.58809 202992_at
## 2 3.943359 7.351799 35.25755 2.600407e-21 4.739242e-17 37.34483 204971_at
## 3 2.318498 6.631187 32.33367 1.785551e-20 1.952500e-16 35.77275 228564_at
## 4 4.905540 8.140399 30.78154 5.320673e-20 4.155826e-16 34.85700 208131_s_at
## 5 4.878195 6.815838 29.02740 1.950627e-19 1.333131e-15 33.74579 210002_at
## 6 4.106051 9.045949 28.82714 2.273193e-19 1.380965e-15 33.61341 212190_at
## symbol change
## 1 C7 up
## 2 CSTA up
## 3 LINC01116 up
## 4 PTGIS up
## 5 GATA6 up
## 6 SERPINE2 up
# Check the number of DE genes
table(deg1$change)##
## down stable up
## 576 19613 623
# Differential expression analysis + plots
dcp1 <- get_deg_all(exp1, Group1, ids1, logFC_cutoff = 2, entriz = TRUE)##
##
## 'select()' returned 1:many mapping between keys and columns
## 138 down genes,130 up genes
dcp1$plotsThis group of functions mainly solves the most basic data-cleaning problems. If the input is handled correctly first, the downstream analysis becomes much easier.
dumd()gives a quick look at an object’s structure.intersect_all()andunion_all()handle intersections and unions of multiple sets.interaction_to_edges()andedges_to_nodes()convert between a “relation table” and a “node table”.cor.full()andcor.one()are used for correlation analysis.
dumd(iris)## colname count
## 3 Petal.Length 43
## 1 Sepal.Length 35
## 2 Sepal.Width 23
## 4 Petal.Width 22
## 5 Species 3
intersect_all(1:5, 3:8, 4:10)## [1] 4 5
union_all(1:5, 3:8, 4:10)## [1] 1 2 3 4 5 6 7 8 9 10
df <- data.frame(
a = c("gene1", "gene2", "gene3"),
b = c("d,f,a,b", "c,e,g", "a,b,d")
)
df## a b
## 1 gene1 d,f,a,b
## 2 gene2 c,e,g
## 3 gene3 a,b,d
interaction_to_edges(df)## a1 a2
## 1 gene1 d
## 2 gene1 f
## 3 gene1 a
## 4 gene1 b
## 5 gene2 c
## 6 gene2 e
## 7 gene2 g
## 8 gene3 a
## 9 gene3 b
## 10 gene3 d
edges_to_nodes(interaction_to_edges(df))## gene type
## 1 gene1 a1
## 2 gene2 a1
## 3 gene3 a1
## 4 d a2
## 5 f a2
## 6 a a2
## 7 b a2
## 8 c a2
## 9 e a2
## 10 g a2
x <- iris[, -5]
cor.full(x)## p.value cor obsnumber
## Sepal.Length:Sepal.Width 1.518983e-01 -0.1175698 150
## Sepal.Length:Petal.Length 1.038667e-47 0.8717538 150
## Sepal.Length:Petal.Width 2.325498e-37 0.8179411 150
## Sepal.Width:Petal.Length 4.513314e-08 -0.4284401 150
## Sepal.Width:Petal.Width 4.073229e-06 -0.3661259 150
## Petal.Length:Petal.Width 4.675004e-86 0.9628654 150
cor.one(x, "Sepal.Width")## p.value cor obsnumber
## Sepal.Length 1.518983e-01 -0.1175698 150
## Petal.Length 4.513314e-08 -0.4284401 150
## Petal.Width 4.073229e-06 -0.3661259 150
The outputs are straightforward:
dumd()gives a quick summary of the object so you can confirm the data type and dimensions first.intersect_all()andunion_all()return the set-operation results.interaction_to_edges()outputs an edge table, andedges_to_nodes()extracts nodes from the edge table.cor.full()returns the overall correlation result, andcor.one()returns the correlation between one variable and the others.
For selected human RNA-seq datasets, GEO provides reprocessed count matrices with Entrez IDs as row names. These matrices can be used directly for transcriptome-wide differential expression analysis.
get_count_txt("GSE162550")## https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=GSE162550&format=file&file=GSE162550_raw_counts_GRCh38.p13_NCBI.tsv.gz
By default, this function only returns the download URL for the count
matrix provided by GEO. You can copy the link into a browser to download
it. If you want to save it locally, set download = TRUE. If the
network is unstable, the download may fail. If the given GSE accession
does not have an official count file, the URL is invalid and the
download will fail.
This section turns an “expression matrix with original IDs” into a matrix that can be used directly in downstream analysis.
trans_array()converts a probe-level expression matrix into a gene-level matrix, and it can also rename rows based on a custom old-to-new ID mapping table.trans_exp_new()can convert matrices with Ensembl or Entrez row names.trans_entrezexp()specifically converts Entrez ID matrices to symbols.trans_ensembl_exp()specifically converts Ensembl ID matrices to gene symbols.make_tcga_group()is used for TCGA sample grouping.
We first create a few example matrices.
# Example array expression matrix
exp_array <- matrix(1:20, nrow = 5)
rownames(exp_array) <- paste0("probe", 1:5)
colnames(exp_array) <- paste0("sample", 1:4)
show_mat(exp_array)## sample1 sample2 sample3 sample4
## probe1 1 6 11 16
## probe2 2 7 12 17
## probe3 3 8 13 18
## probe4 4 9 14 19
ids_array <- data.frame(
probe_id = c("probe1", "probe2", "probe3", "probe4", "probe5"),
symbol = c("A", "B", "B", "C", "C")
)
head(ids_array)## probe_id symbol
## 1 probe1 A
## 2 probe2 B
## 3 probe3 B
## 4 probe4 C
## 5 probe5 C
# Example Ensembl expression matrix
exp_ensembl <- matrix(rnorm(200), nrow = 20)
rownames(exp_ensembl) <- c(
sample(mRNA_anno$gene_id, 10, replace = TRUE),
sample(lnc_anno$gene_id, 10, replace = TRUE)
)
colnames(exp_ensembl) <- paste0("s", 1:10)
show_mat(exp_ensembl)## s1 s2 s3 s4
## ENSG00000105519.12 -0.76670821 -0.6530894 2.1291767 2.64591554
## ENSG00000075043.17 0.89973925 -1.0277558 0.1924601 -0.48781338
## ENSG00000121289.17 0.31039150 -0.2447213 0.1996285 0.35487582
## ENSG00000235631.1 0.00426799 0.5441694 2.3501687 0.08757415
# Example Entrez expression matrix
gene_demo <- sample(genes, 20, replace = TRUE)
exp_entrez <- matrix(rnorm(length(gene_demo) * 6), ncol = 6)
rownames(exp_entrez) <- gene_demo
colnames(exp_entrez) <- paste0("s", 1:6)
show_mat(exp_entrez)## s1 s2 s3 s4
## 140679 0.3831703 -0.1183658 0.24566685 -0.5400305
## 51046 0.9888129 1.0435844 1.07239869 0.1997114
## 56853 -0.1362858 0.5830176 1.19305547 0.2940977
## 166752 -1.3732949 -0.2314533 0.07126605 -0.2882327
trans_array() is suitable when you already have a probe_id -> symbol
mapping table. It rewrites the row names of the expression matrix
according to that table and returns a new matrix. If one symbol
corresponds to multiple probes, only one row is kept by default; if you
set average = TRUE, those duplicated entries are merged by mean.
trans_array(exp_array, ids_array)## 3 rownames transformed after duplicate rows removed
## sample1 sample2 sample3 sample4
## A 1 6 11 16
## B 2 7 12 17
## C 4 9 14 19
trans_array(exp_array, ids_array, average = TRUE)## 3 rownames transformed after averaging duplicate rows
## sample1 sample2 sample3 sample4
## A 1.0 6.0 11.0 16.0
## B 2.5 7.5 12.5 17.5
## C 4.5 9.5 14.5 19.5
The first result is a matrix that keeps only unique symbols, and the
second result is a matrix that averages duplicate mappings. Both can be
used for downstream differential expression analysis or plotting; they only differ
in how repeated information is handled. The other trans_*() functions
described below also support the average argument.
s1 = trans_exp_new(exp_ensembl)## 'select()' returned 1:many mapping between keys and columns
## 16 rownames transformed after duplicate rows removed
show_mat(s1)## s1 s2 s3 s4
## CAPS -0.76670821 -0.6530894 2.1291767 2.64591554
## KCNQ2 0.89973925 -1.0277558 0.1924601 -0.48781338
## CEP89 0.31039150 -0.2447213 0.1996285 0.35487582
## RNF148 0.00426799 0.5441694 2.3501687 0.08757415
It also supports the species argument, which can be set to human,
mouse, or rat.
It also supports the type argument, which can handle Ensembl (default,
equivalent to trans_ensembl_exp) and Entrez (equivalent to
trans_entrezexp) IDs.
- type = “entrez”
s2 = trans_exp_new(exp_entrez,type = "entrez")## 'select()' returned 1:1 mapping between keys and columns
## 19 rownames transformed after duplicate rows removed
show_mat(s2)## s1 s2 s3 s4
## SLC32A1 0.3831703 -0.1183658 0.24566685 -0.5400305
## ST8SIA3 0.9888129 1.0435844 1.07239869 0.1997114
## CELF4 -0.1362858 0.5830176 1.19305547 0.2940977
## FREM3 -1.3732949 -0.2314533 0.07126605 -0.2882327
- type = “array” with an explicit
idstable, which converts IDs just liketrans_array().
s3 = trans_exp_new(exp_array,type = "array",ids = ids_array)## 3 rownames transformed after duplicate rows removed
show_mat(s3)## sample1 sample2 sample3 sample4
## A 1 6 11 16
## B 2 7 12 17
## C 4 9 14 19
s4 = trans_ensembl_exp(exp_ensembl)## 'select()' returned 1:many mapping between keys and columns
## 16 rownames transformed after duplicate rows removed
show_mat(s4)## s1 s2 s3 s4
## CAPS -0.76670821 -0.6530894 2.1291767 2.64591554
## KCNQ2 0.89973925 -1.0277558 0.1924601 -0.48781338
## CEP89 0.31039150 -0.2447213 0.1996285 0.35487582
## RNF148 0.00426799 0.5441694 2.3501687 0.08757415
s5 = trans_entrezexp(exp_entrez)## 'select()' returned 1:1 mapping between keys and columns
## 19 rownames transformed after duplicate rows removed
show_mat(s5)## s1 s2 s3 s4
## SLC32A1 0.3831703 -0.1183658 0.24566685 -0.5400305
## ST8SIA3 0.9888129 1.0435844 1.07239869 0.1997114
## CELF4 -0.1362858 0.5830176 1.19305547 0.2940977
## FREM3 -1.3732949 -0.2314533 0.07126605 -0.2882327
make_tcga_group() can accept either an expression matrix or a vector
of sample names. It splits TCGA samples into normal and tumor
according to TCGA sample barcodes and returns a factor.
show_mat(exp_hub1)## GTEX-S33H-1226-SM-4AD69 GTEX-VJYA-0826-SM-4KL1M TCGA-FB-A545-01
## CXCL8 34.00041 50.99857 1341.029
## FN1 591.01914 960.00492 137998.164
## COL3A1 2226.94824 3682.08761 177553.205
## ISG15 125.00175 135.00350 1784.980
## GTEX-ZF3C-2026-SM-4WWB5
## CXCL8 388.9918
## FN1 5218.1531
## COL3A1 15621.0916
## ISG15 513.9898
make_tcga_group(exp_hub1[,1:10])## [1] normal normal tumor normal normal normal tumor tumor tumor tumor
## Levels: normal tumor
make_tcga_group(colnames(exp_hub1)[1:10])## [1] normal normal tumor normal normal normal tumor tumor tumor tumor
## Levels: normal tumor
The output of make_tcga_group() is a normal / tumor factor that
can be used directly for differential expression analysis.
This chapter collects the most common visualization functions in one place so that you can find the right function by plot type.
draw_volcano() accepts a differential expression results table and
uses adjusted P values (padj) by default.
draw_volcano(deseq_data)
These functions all belong to the category of “turn a matrix into a plot”:
draw_pca()performs principal component analysis.draw_tsne()performs tSNE dimensionality reduction.draw_heatmap()draws a heatmap.draw_venn()draws a Venn diagram.
draw_pca(t(iris[, 1:4]), iris$Species)
draw_pca(t(iris[, 1:4]), iris$Species, style = "ggplot2")
draw_pca(t(iris[, 1:4]), iris$Species, style = "3D")
exp_tsne <- matrix(rnorm(50 * 60), nrow = 50)
group_tsne <- factor(rep(c("A", "B"), each = 30))
draw_tsne(exp_tsne, group_tsne, perplexity = 15)
exp_heat <- matrix(abs(rnorm(60, sd = 16)), nrow = 10)
exp_heat[, 4:6] <- exp_heat[, 4:6] + 20
colnames(exp_heat) <- paste0("sample", 1:6)
rownames(exp_heat) <- paste0("gene", 1:10)
group_heat <- factor(rep(c("A", "B"), each = 3))
draw_heatmap(exp_heat, group_heat)
draw_heatmap(exp_heat, group_heat, split_column = TRUE, show_column_title = TRUE)
xset <- list(
Deseq2 = sample(1:100, 30),
edgeR = sample(1:100, 30),
limma = sample(1:100, 30)
)
draw_venn(xset, "test")
draw_venn(xset, "test", color = c("darkgreen", "darkblue", "#B2182B"))This group is better for looking at the distribution of one or a few
genes across different groups. draw_boxplot() and exp_boxplot()
mainly draw grouped boxplots, corheatmap() and corscatterplot() are
for correlation analysis, and ggheat() is a ggplot-style heatmap.
draw_boxplot(t(iris[, 1:4]), iris$Species)
draw_boxplot(exp_heat, group_heat)
exp_boxplot(exprSet_hub1)[[1]]
xg <- rownames(exprSet_hub1)[1:3]
yg <- rownames(exprSet_hub1)[4:7]
corheatmap(exprSet_hub1, xg, yg)
corscatterplot(iris, "Sepal.Length", "Sepal.Width")
exp_heat_gg <- matrix(abs(rnorm(60, sd = 16)), nrow = 10)
exp_heat_gg[, 4:6] <- exp_heat_gg[, 4:6] + 20
rownames(exp_heat_gg) <- paste0("sample", 1:10)
colnames(exp_heat_gg) <- paste0("gene", 1:6)
group_heat_gg <- factor(rep(c("A", "B"), each = 5))
ggheat(exp_heat_gg, group_heat_gg)
draw_boxplot() returns a boxplot for quickly comparing the expression
of a few genes across groups. corheatmap()
returns a correlation heatmap, and corscatterplot() returns a pairwise
scatter correlation plot.
box_surv() combines an expression boxplot and a survival plot, which
is useful for single-gene visualization. Rather than producing a
standalone plot, it places differences in expression and survival on the
same page for direct comparison.
box_surv(
log2(exp_hub1+1),
exprSet_hub1,
meta1
)[[1]]
This section revolves around KM curves, Cox regression, and optimal
cut-point estimation. point_cut() is used to find the cut point first,
surv_KM() and surv_cox() are the two most commonly used
survival-analysis entry points, exp_surv() bundles survival plots for
multiple genes, and risk_plot() is often used for risk-score display.
point_cut(exprSet_hub1, meta1)## 8 rows with less than 5 values were ignored in cut.point calculations
## $CXCL8
## [1] 9.63068
##
## $FN1
## [1] 18.12339
##
## $COL3A1
## [1] 18.78103
##
## $ISG15
## [1] 11.28338
##
## $COL1A2
## [1] 18.63055
##
## $CXCL10
## [1] 6.884411
##
## $ICAM1
## [1] 12.27501
##
## $MMP9
## [1] 10.58174
surv_KM(exprSet_hub1, meta1)## ICAM1
## 0.0218004
surv_KM(exprSet_hub1, meta1, cut.point = TRUE)## 8 rows with less than 5 values were ignored in cut.point calculations
## CXCL8 ICAM1 FN1
## 0.005832445 0.010023406 0.038235032
surv_cox(exprSet_hub1, meta1)## coef se z p HR HRse HRz
## ICAM1 -0.4813703 0.2119765 -2.270866 0.02315511 0.6179361 0.130988 -2.916787
## HRp HRCILL HRCIUL
## ICAM1 0.003536574 0.4078578 0.9362209
surv_cox(exprSet_hub1, meta1, cut.point = TRUE)## 8 rows with less than 5 values were ignored in cut.point calculations
## coef se z p HR HRse HRz
## CXCL8 1.1142476 0.4246572 2.623876 0.008693553 3.0472747 1.2940471 1.582071
## FN1 0.4730967 0.2304361 2.053049 0.040067814 1.6049566 0.3698400 1.635725
## ICAM1 -0.5388259 0.2118283 -2.543692 0.010968763 0.5834328 0.1235876 -3.370624
## HRp HRCILL HRCIUL
## CXCL8 0.1136333125 1.3256923 7.0045535
## FN1 0.1018971198 1.0216819 2.5212208
## ICAM1 0.0007499825 0.3851965 0.8836889
tmp <- exp_surv(exprSet_hub1, meta1)
patchwork::wrap_plots(tmp) + patchwork::plot_layout(guides = "collect")
box_surv(exprSet_hub1, exprSet_hub1, meta1)[[1]]
risk_plot(exprSet_hub1, meta1, riskscore = rnorm(nrow(meta1)))
lung_demo <- survival::lung
lung_demo$event <- lung_demo$status == 2
lung_demo$sex <- factor(lung_demo$sex, levels = c(1, 2), labels = c("male", "female"))
draw_KM(meta = lung_demo, group_list = lung_demo$sex, event_col = "event")
Most outputs here are plot objects. surv_KM() returns a Kaplan-Meier
survival curve, surv_cox() returns a Cox-analysis plot or table, and
exp_surv() returns a plot list that can be combined. draw_KM() is
more of a general-purpose survival curve function: as long as you
provide an event column and a grouping column, it can draw the plot.
This group of functions mainly serves lncRNA-miRNA-mRNA network
analysis. plcortest() first screens for correlated lncRNA-mRNA pairs,
after which hypertest() uses a hypergeometric test to assess shared
miRNAs. The former returns a list, and the latter returns a result table.
lnc_exp <- matrix(rnorm(15), nrow = 3)
colnames(lnc_exp) <- paste0("S", 1:5)
rownames(lnc_exp) <- paste0("lnc", 1:3)
mRNA_exp <- matrix(rnorm(25), nrow = 5)
colnames(mRNA_exp) <- paste0("S", 1:5)
rownames(mRNA_exp) <- paste0("gene", 1:5)
plcortest(lnc_exp, mRNA_exp, cor_cutoff = 0)## $lnc1
## logical(0)
##
## $lnc2
## logical(0)
##
## $lnc3
## logical(0)
lnctarget <- data.frame(
lncRNA = c("lnc1", "lnc1", "lnc2"),
miRNA = c("miR1", "miR2", "miR2")
)
pctarget <- data.frame(
mRNA = c("gene1", "gene2", "gene2"),
miRNA = c("miR1", "miR2", "miR3")
)
hypertest(c("lnc1", "lnc2"), c("gene1", "gene2"), deMIR = c("miR2"), lnctarget = lnctarget, pctarget = pctarget)## lncRNAs Genes Counts listTotal popHits popTotal foldEnrichment
## 1 lnc1 gene1 1 2 1 3 1.5
## 4 lnc2 gene2 1 1 2 3 1.5
## 2 lnc1 gene2 1 2 2 3 0.75
## hyperPValue miRNAs deMIRCounts deMIRs
## 1 0.666666666666667 miR1 0
## 4 0.666666666666667 miR2 1 miR2
## 2 1 miR2 1 miR2
plcortest() returns a list named by lncRNA, and each element contains
the mRNA names correlated with it. hypertest() returns a data.frame
containing the lncRNA and mRNA identifiers, the number of shared miRNAs,
fold enrichment, and the P value from the hypergeometric test.
This section mainly covers KEGG and GO enrichment. quick_enrich()
returns a list that usually contains kk, go, and the corresponding
dotplot objects; if an enrichment analysis returns no enriched terms,
the function reports this instead of attempting to draw an empty plot.
double_enrich() is suitable for enriching up- and down-regulated genes
separately and comparing them in one figure.
g <- quick_enrich(genes)## Reading KEGG annotation online: "https://rest.kegg.jp/link/hsa/pathway"...
## Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/hsa"...
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
names(g)## [1] "kk" "go" "kk.dot" "go.dot"
g$kk## #
## # over-representation test
## #
## #...@organism hsa
## #...@ontology KEGG
## #...@keytype ENTREZID
## #...@gene chr [1:511] "960" "3371" "10763" "51029" "2146" "2149" "2316" "7153" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...54 enriched terms found
## 'data.frame': 54 obs. of 14 variables:
## $ category : chr NA "Human Diseases" "Organismal Systems" "Human Diseases" ...
## $ subcategory : chr NA "Substance dependence" "Nervous system" "Substance dependence" ...
## $ ID : chr "hsa04082" "hsa05032" "hsa04727" "hsa05033" ...
## $ Description : chr "Neuroactive ligand signaling" "Morphine addiction" "GABAergic synapse" "Nicotine addiction" ...
## $ GeneRatio : chr "30/236" "20/236" "19/236" "14/236" ...
## $ BgRatio : chr "199/9421" "91/9421" "89/9421" "41/9421" ...
## $ RichFactor : num 0.151 0.22 0.213 0.341 0.114 ...
## $ FoldEnrichment: num 6.02 8.77 8.52 13.63 4.56 ...
## $ zScore : num 11.47 11.94 11.43 12.99 9.21 ...
## $ pvalue : num 1.27e-15 5.72e-14 4.38e-13 5.04e-13 5.68e-12 ...
## $ p.adjust : num 3.25e-13 7.29e-12 3.22e-11 3.22e-11 2.90e-10 ...
## $ qvalue : num 2.39e-13 5.36e-12 2.36e-11 2.36e-11 2.13e-10 ...
## $ geneID : chr "OPRK1/DRD1/ADCY2/GABRB2/GRM5/HTR2C/CHRM3/GABRG1/GRIN2A/GABRB3/PRKCB/HRH3/GNG3/SLC32A1/GNB5/GRIN3A/HTR2A/SLC17A6"| __truncated__ "DRD1/KCNJ3/ADCY2/GABRB2/GABRG1/GABRB3/PRKCB/PDE1A/GNG3/SLC32A1/GNB5/GNAO1/GABRD/PDE2A/KCNJ9/GABRA5/ADCY5/GABRA1/CACNA1B/GABRG2" "ADCY2/GABRB2/GABRG1/GABRB3/PRKCB/SLC12A5/GNG3/SLC32A1/GNB5/GNAO1/GABRD/GLS2/GABRA5/ADCY5/GAD1/GABRA1/CACNA1B/GABRG2/GAD2" "GABRB2/GABRG1/GRIN2A/GABRB3/SLC32A1/GRIN3A/SLC17A6/SLC17A7/GABRD/GABRA5/GRIN2B/GABRA1/CACNA1B/GABRG2" ...
## $ Count : int 30 20 19 14 29 22 18 15 15 15 ...
## #...Citation
## S Xu, E Hu, Y Cai, Z Xie, X Luo, L Zhan, W Tang, Q Wang, B Liu, R Wang, W Xie, T Wu, L Xie, G Yu. Using clusterProfiler to characterize multiomics data. Nature Protocols. 2024, 19(11):3292-3320
g$go## #
## # over-representation test
## #
## #...@organism Homo sapiens
## #...@ontology GOALL
## #...@keytype ENTREZID
## #...@gene chr [1:511] "960" "3371" "10763" "51029" "2146" "2149" "2316" "7153" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...564 enriched terms found
## 'data.frame': 564 obs. of 13 variables:
## $ ONTOLOGY : chr "BP" "BP" "BP" "BP" ...
## $ ID : chr "GO:0006836" "GO:0099504" "GO:0099003" "GO:0007269" ...
## $ Description : chr "neurotransmitter transport" "synaptic vesicle cycle" "vesicle-mediated transport in synapse" "neurotransmitter secretion" ...
## $ GeneRatio : chr "38/466" "38/466" "40/466" "31/466" ...
## $ BgRatio : chr "216/18860" "220/18860" "259/18860" "154/18860" ...
## $ RichFactor : num 0.176 0.173 0.154 0.201 0.201 ...
## $ FoldEnrichment: num 7.12 6.99 6.25 8.15 8.15 ...
## $ zScore : num 14.4 14.2 13.5 14.2 14.2 ...
## $ pvalue : num 1.05e-21 2.06e-21 1.23e-20 1.06e-19 1.06e-19 ...
## $ p.adjust : num 3.98e-18 3.98e-18 1.59e-17 8.20e-17 8.20e-17 ...
## $ qvalue : num 3.30e-18 3.30e-18 1.32e-17 6.79e-17 6.79e-17 ...
## $ geneID : chr "GFAP/SLC6A15/DRD1/RIMS3/CAMK2A/MEF2C/SYT4/SV2B/SNCG/NAPB/PRKCB/SYT13/SYT7/HRH3/PCLO/RIMS1/SYN2/SYT1/STXBP5/PPFI"| __truncated__ "RIMS3/ATP6V1G2/SYT4/SV2B/SNCG/NAPB/PRKCB/SYT13/DNM3/SYT7/PCLO/RIMS1/SYN2/SYT1/STXBP5/PPFIA2/RAPGEF4/SLC32A1/SLC"| __truncated__ "RIMS3/ATP6V1G2/SYT4/SV2B/SNCG/NAPB/AAK1/PRKCB/SYT13/DNM3/SYT7/PCLO/RIMS1/SYN2/SYT1/STXBP5/PPFIA2/RAPGEF4/SLC32A"| __truncated__ "RIMS3/CAMK2A/MEF2C/SYT4/SV2B/SNCG/NAPB/PRKCB/SYT13/SYT7/HRH3/PCLO/RIMS1/SYN2/SYT1/STXBP5/PPFIA2/SLC32A1/UNC13C/"| __truncated__ ...
## $ Count : int 38 38 40 31 31 39 41 47 27 23 ...
## #...Citation
## S Xu, E Hu, Y Cai, Z Xie, X Luo, L Zhan, W Tang, Q Wang, B Liu, R Wang, W Xie, T Wu, L Xie, G Yu. Using clusterProfiler to characterize multiomics data. Nature Protocols. 2024, 19(11):3292-3320
g$kk.dot
g$go.dot
deg_enrich <- data.frame(
ENTREZID = genes[1:20],
change = rep(c("up", "down"), 10)
)
double_enrich(deg_enrich)## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
## <tinyarray_double_enrich>
## components: kp, gp
You can first use names(g) to inspect the structure of
quick_enrich() output, and then check g$kk and g$go separately.
double_enrich() returns an object with a custom class, and its print
method summarizes the available components.
This version focuses on reducing manual data wrangling, improving the reliability of annotation conversion, and making the plotting workflow more convenient. The new functions are introduced first, followed by changes to existing functions.
get_ids()retrieves the platform annotation table and returns a standardized two-column table with the fixed column namesprobe_idandsymbol. The result can be used directly for expression-matrix conversion, differential expression analysis, and enrichment analysis.
-
geo_download()now returnsgroup_candidatesin addition toexp,pd, andgpl. This makes it easier to perform automatic grouping or manually select grouping information after downloading GEO data. -
get_deg()andget_deg_all()now supportids = NULL. If the input is already a gene-level expression matrix, you can run differential expression analysis directly without providing a probe-to-gene mapping table. -
trans_array(),trans_exp_new(), andtrans_entrezexp()now supportaverage = TRUE. When multiple probes or IDs map to the samesymbol, their expression values can be averaged so that the resulting gene symbols remain unique. -
draw_volcano()uses adjusted P values (padj) by default. To use raw P values (P.value) instead, explicitly setadjust = FALSE. -
get_deg_all()now issues a warning and skips the heatmap when there are no differential genes available for plotting, instead of stopping with an error. -
make_tcga_group()can now accept either an expression matrix or a vector of sample names, which makes it more flexible to use.