✍️ 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.934if (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 candidatesThis 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 10Group1 = 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
## 43086head(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 UBA7get_ids() returns a cleaned two-column table:
probe_idsymbolIt automatically does several things:
AnnoProbe::idmap();NA or an empty string
"";probe_id and
symbol;probe_id to character;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=dataget_gpl_txt("GPL6244")## https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL6244&targ=self&form=text&view=dataDownload 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 genesdcp1$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() and union_all() handle
intersections and unions of multiple sets.interaction_to_edges() and
edges_to_nodes() convert between a “relation table” and a
“node table”.cor.full() and cor.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 3intersect_all(1:5, 3:8, 4:10)## [1] 4 5union_all(1:5, 3:8, 4:10)## [1] 1 2 3 4 5 6 7 8 9 10df <- 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,dinteraction_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 dedges_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 a2x <- 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 150cor.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 150The outputs are straightforward:
dumd() gives a quick summary of the object so you can
confirm the data type and dimensions first.intersect_all() and union_all() return the
set-operation results.interaction_to_edges() outputs an edge table, and
edges_to_nodes() extracts nodes from the edge table.cor.full() returns the overall correlation result, and
cor.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.gzBy 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 19ids_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.2882327trans_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 19trans_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.5The 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 removedshow_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.08757415It 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.
s2 = trans_exp_new(exp_entrez,type = "entrez")## 'select()' returned 1:1 mapping between keys and columns
## 19 rownames transformed after duplicate rows removedshow_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.2882327ids table, which
converts IDs just like trans_array().s3 = trans_exp_new(exp_array,type = "array",ids = ids_array)## 3 rownames transformed after duplicate rows removedshow_mat(s3)## sample1 sample2 sample3 sample4
## A 1 6 11 16
## B 2 7 12 17
## C 4 9 14 19s4 = trans_ensembl_exp(exp_ensembl)## 'select()' returned 1:many mapping between keys and columns
## 16 rownames transformed after duplicate rows removedshow_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.08757415s5 = trans_entrezexp(exp_entrez)## 'select()' returned 1:1 mapping between keys and columns
## 19 rownames transformed after duplicate rows removedshow_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.2882327make_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.9898make_tcga_group(exp_hub1[,1:10])## [1] normal normal tumor normal normal normal tumor tumor tumor tumor
## Levels: normal tumormake_tcga_group(colnames(exp_hub1)[1:10])## [1] normal normal tumor normal normal normal tumor tumor tumor tumor
## Levels: normal tumorThe 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.58174surv_KM(exprSet_hub1, meta1)## ICAM1
## 0.0218004surv_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.038235032surv_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.9362209surv_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.8836889tmp <- 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 miR2plcortest() 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 columnsnames(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-3320g$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-3320g$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, gpYou 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 names
probe_id and symbol. The result can be used
directly for expression-matrix conversion, differential expression
analysis, and enrichment analysis.geo_download() now returns
group_candidates in addition to exp,
pd, and gpl. This makes it easier to perform
automatic grouping or manually select grouping information after
downloading GEO data.
get_deg() and get_deg_all() now support
ids = 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(), and
trans_entrezexp() now support average = TRUE.
When multiple probes or IDs map to the same symbol, 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 set
adjust = 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.