library(yulab.utils)
pload(airway)
pload(pheatmap)
pload(DESeq2)
pload(SummarizedExperiment)
pload(org.Hs.eg.db)
pload(clusterProfiler)
pload(enrichplot)
pload(ggplot2)
pload(dplyr)转录组分析
R 环境
环境准备
yulab.utils::pload() 会检查 R 包是否已安装；缺少时提示安装，安装后自动加载。
glmpca、tximeta、AnnotationHub 和翻译/交互扩展为可选包；未安装 glmpca 时，PCA 小节使用 DESeq2 的 vst() 备用路线。
示例数据
library(yulab.utils)
pload(airway)── loading the package: airway ──dir <- system.file("extdata", package="airway", mustWork=TRUE)
list.files(dir) [1] "GSE52778_series_matrix.txt" "Homo_sapiens.GRCh37.75_subset.gtf"
[3] "quants" "sample_table.csv"
[5] "SraRunInfo_SRP033351.csv" "SRR1039508_subset.bam"
[7] "SRR1039509_subset.bam" "SRR1039512_subset.bam"
[9] "SRR1039513_subset.bam" "SRR1039516_subset.bam"
[11] "SRR1039517_subset.bam" "SRR1039520_subset.bam"
[13] "SRR1039521_subset.bam" list.files(file.path(dir, "quants"))[1] "SRR1039508" "SRR1039509"csvfile <- file.path(dir, "sample_table.csv")
coldata <- read.csv(csvfile, row.names=1, stringsAsFactors=FALSE)
coldata SampleName cell dex albut Run avgLength Experiment
SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126 SRX384345
SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126 SRX384346
SRR1039512 GSM1275866 N052611 untrt untrt SRR1039512 126 SRX384349
SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87 SRX384350
SRR1039516 GSM1275870 N080611 untrt untrt SRR1039516 120 SRX384353
SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126 SRX384354
SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520 101 SRX384357
SRR1039521 GSM1275875 N061011 trt untrt SRR1039521 98 SRX384358
Sample BioSample
SRR1039508 SRS508568 SAMN02422669
SRR1039509 SRS508567 SAMN02422675
SRR1039512 SRS508571 SAMN02422678
SRR1039513 SRS508572 SAMN02422670
SRR1039516 SRS508575 SAMN02422682
SRR1039517 SRS508576 SAMN02422673
SRR1039520 SRS508579 SAMN02422683
SRR1039521 SRS508580 SAMN02422677coldata <- coldata[1:2,]
coldata$names <- coldata$Run
coldata$files <- file.path(dir, "quants", coldata$names, "quant.sf.gz")
file.exists(coldata$files)[1] TRUE TRUEcoldata SampleName cell dex albut Run avgLength Experiment
SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126 SRX384345
SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126 SRX384346
Sample BioSample names
SRR1039508 SRS508568 SAMN02422669 SRR1039508
SRR1039509 SRS508567 SAMN02422675 SRR1039509
files
SRR1039508 /home/wang/Rlib/airway/extdata/quants/SRR1039508/quant.sf.gz
SRR1039509 /home/wang/Rlib/airway/extdata/quants/SRR1039509/quant.sf.gz定量
Code
pload("tximeta")
se <- tximeta(coldata)
dim(se)
head(rownames(se))
gse <- summarizeToGene(se)
dim(gse)上面折叠代码运行时间较长，可以直接加载airway包中的示例数据gse。
data(gse, package="airway")
gseclass: RangedSummarizedExperiment
dim: 58294 8
metadata(6): tximetaInfo quantInfo ... txomeInfo txdbInfo
assays(3): counts abundance length
rownames(58294): ENSG00000000003.14 ENSG00000000005.5 ...
ENSG00000285993.1 ENSG00000285994.1
rowData names(1): gene_id
colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
colData names(3): names donor conditiongse是一个存储为SummarizedExperiment结构的数据，其中的基因定量丰度表, 样本信息以及基因信息的访问，具体操作如下。
访问相应的信息
pload(SummarizedExperiment)── loading the package: SummarizedExperiment ──head(rownames(gse))[1] "ENSG00000000003.14" "ENSG00000000005.5" "ENSG00000000419.12"
[4] "ENSG00000000457.13" "ENSG00000000460.16" "ENSG00000000938.12"assayNames(gse)[1] "counts" "abundance" "length" head(assay(gse), 3) SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
ENSG00000000003.14 708.164 467.962 900.992 424.368 1188.295
ENSG00000000005.5 0.000 0.000 0.000 0.000 0.000
ENSG00000000419.12 455.000 510.000 604.000 352.000 583.000
SRR1039517 SRR1039520 SRR1039521
ENSG00000000003.14 1090.668 805.929 599.337
ENSG00000000005.5 0.000 0.000 0.000
ENSG00000000419.12 773.999 409.999 499.000rowRanges(gse)GRanges object with 58294 ranges and 1 metadata column:
seqnames ranges strand | gene_id
<Rle> <IRanges> <Rle> | <character>
ENSG00000000003.14 chrX 100627109-100639991 - | ENSG00000000003.14
ENSG00000000005.5 chrX 100584802-100599885 + | ENSG00000000005.5
ENSG00000000419.12 chr20 50934867-50958555 - | ENSG00000000419.12
ENSG00000000457.13 chr1 169849631-169894267 - | ENSG00000000457.13
ENSG00000000460.16 chr1 169662007-169854080 + | ENSG00000000460.16
... ... ... ... . ...
ENSG00000285990.1 chr14 19244904-19269380 - | ENSG00000285990.1
ENSG00000285991.1 chr6 149817937-149896011 - | ENSG00000285991.1
ENSG00000285992.1 chr8 47129262-47132628 + | ENSG00000285992.1
ENSG00000285993.1 chr18 46409197-46410645 - | ENSG00000285993.1
ENSG00000285994.1 chr10 12563151-12567351 + | ENSG00000285994.1
-------
seqinfo: 25 sequences (1 circular) from hg38 genomeseqinfo(rowRanges(gse))Seqinfo object with 25 sequences (1 circular) from hg38 genome:
seqnames seqlengths isCircular genome
chr1 248956422 FALSE hg38
chr2 242193529 FALSE hg38
chr3 198295559 FALSE hg38
chr4 190214555 FALSE hg38
chr5 181538259 FALSE hg38
... ... ... ...
chr21 46709983 FALSE hg38
chr22 50818468 FALSE hg38
chrX 156040895 FALSE hg38
chrY 57227415 FALSE hg38
chrM 16569 TRUE hg38colData(gse)DataFrame with 8 rows and 3 columns
names donor condition
<factor> <factor> <factor>
SRR1039508 SRR1039508 N61311 Untreated
SRR1039509 SRR1039509 N61311 Dexamethasone
SRR1039512 SRR1039512 N052611 Untreated
SRR1039513 SRR1039513 N052611 Dexamethasone
SRR1039516 SRR1039516 N080611 Untreated
SRR1039517 SRR1039517 N080611 Dexamethasone
SRR1039520 SRR1039520 N061011 Untreated
SRR1039521 SRR1039521 N061011 Dexamethasone差异分析
创建DESeqDataSet对象
pload(DESeq2)── loading the package: DESeq2 ──dds <- DESeqDataSet(gse, design = ~ donor + condition)using counts and average transcript lengths from tximetamodel.matrix(~ donor + condition, data = as.data.frame(colData(dds))) (Intercept) donorN061011 donorN080611 donorN61311
SRR1039508 1 0 0 1
SRR1039509 1 0 0 1
SRR1039512 1 0 0 0
SRR1039513 1 0 0 0
SRR1039516 1 0 1 0
SRR1039517 1 0 1 0
SRR1039520 1 1 0 0
SRR1039521 1 1 0 0
conditionDexamethasone
SRR1039508 0
SRR1039509 1
SRR1039512 0
SRR1039513 1
SRR1039516 0
SRR1039517 1
SRR1039520 0
SRR1039521 1
attr(,"assign")
[1] 0 1 1 1 2
attr(,"contrasts")
attr(,"contrasts")$donor
[1] "contr.treatment"
attr(,"contrasts")$condition
[1] "contr.treatment"hist(log(rowSums(counts(dds))))
nrow(dds)[1] 58294keep <- rowSums(counts(dds)) > 1
dds <- dds[keep, ]
nrow(dds)[1] 31604先检查 design formula 和样本元数据是否一致：donor 用来吸收配对/个体差异，condition 才是本次比较的主效应。若实验设计改变，不能只改 contrast 名称而不重新检查模型矩阵。
PCA
if (requireNamespace("glmpca", quietly = TRUE)) {
pload(glmpca)
gpca <- glmpca(counts(dds), L = 2)
gpca.dat <- gpca$factors
gpca.dat$condition <- dds$condition
gpca.dat$donor <- dds$donor
head(gpca.dat)
pload(ggplot2)
ggplot(gpca.dat, aes(x = dim1, y = dim2,
color = condition, shape = donor)) +
geom_point(size = 3) +
ggtitle("glmpca - Generalized PCA")
} else {
message("未安装 glmpca，使用 DESeq2::vst() + prcomp() 完成 PCA。")
vsd <- vst(dds, blind = TRUE)
pca <- prcomp(t(assay(vsd)))
pca.dat <- as.data.frame(pca$x[, 1:2])
pca.dat$condition <- colData(vsd)$condition
pca.dat$donor <- colData(vsd)$donor
pload(ggplot2)
ggplot(pca.dat, aes(x = PC1, y = PC2,
color = condition, shape = donor)) +
geom_point(size = 3) +
ggtitle("PCA after variance-stabilizing transformation")
}未安装 glmpca，使用 DESeq2::vst() + prcomp() 完成 PCA。using 'avgTxLength' from assays(dds), correcting for library size── loading the package: ggplot2 ──
样本相似性
PCA 只把主要变异压到两个轴上，样本之间的整体相似性更适合用相关系数矩阵看：同一 donor 的两次处理应该比跨 donor 更接近，condition 的效应则应该在同一 donor 内部体现。这里用 vst() 之后的表达矩阵算 Pearson 相关（vst 让不同表达量级的基因可比，直接对 counts 求相关会被高表达基因主导）。
pload(pheatmap)── loading the package: pheatmap ──vsd <- vst(dds, blind = TRUE)using 'avgTxLength' from assays(dds), correcting for library sizesample_cor <- cor(assay(vsd))
anno_col <- as.data.frame(colData(dds)[, c("donor", "condition")])
pheatmap(sample_cor,
annotation_col = anno_col,
annotation_row = anno_col,
display_numbers = TRUE,
number_format = "%.3f",
main = "样本间相关性（vst 表达矩阵）")
读这张图的顺序：先看对角线上的 1，再看组内与组间的差别。如果某个样本与所有样本的相关都明显偏低，先回去查 QC 与建库记录，不要直接删（第 4 讲的样本审问一节讲的就是这件事）。
差异表达分析
dds <- DESeq(dds)estimating size factorsusing 'avgTxLength' from assays(dds), correcting for library sizeestimating dispersionsgene-wise dispersion estimatesmean-dispersion relationshipfinal dispersion estimatesfitting model and testingcoef_name <- grep("^condition_", resultsNames(dds), value = TRUE)
coef_name[1] "condition_Dexamethasone_vs_Untreated"res <- results(dds, name = coef_name)
res_shrunk <- lfcShrink(dds, coef = coef_name, type = "normal")using 'normal' for LFC shrinkage, the Normal prior from Love et al (2014).
Note that type='apeglm' and type='ashr' have shown to have less bias than type='normal'.
See ?lfcShrink for more details on shrinkage type, and the DESeq2 vignette.
Reference: https://doi.org/10.1093/bioinformatics/bty895reslog2 fold change (MLE): condition Dexamethasone vs Untreated
Wald test p-value: condition Dexamethasone vs Untreated
DataFrame with 31604 rows and 6 columns
baseMean log2FoldChange lfcSE stat pvalue
<numeric> <numeric> <numeric> <numeric> <numeric>
ENSG00000000003.14 739.940717 -0.3611537 0.106869 -3.379419 0.000726392
ENSG00000000419.12 511.735722 0.2063147 0.128665 1.603509 0.108822318
ENSG00000000457.13 314.194855 0.0378308 0.158633 0.238479 0.811509461
ENSG00000000460.16 79.793622 -0.1152590 0.314991 -0.365912 0.714430444
ENSG00000000938.12 0.307267 -1.3691185 3.503764 -0.390757 0.695977205
... ... ... ... ... ...
ENSG00000285979.1 38.353886 0.3423657 0.359511 0.952310 0.340940
ENSG00000285987.1 1.562508 0.7064145 1.547295 0.456548 0.647996
ENSG00000285990.1 0.642315 0.3647333 3.433276 0.106235 0.915396
ENSG00000285991.1 11.276284 -0.1165515 0.748601 -0.155692 0.876275
ENSG00000285994.1 3.651041 -0.0960094 1.068660 -0.089841 0.928414
padj
<numeric>
ENSG00000000003.14 0.00531137
ENSG00000000419.12 0.29318870
ENSG00000000457.13 0.92255697
ENSG00000000460.16 0.87298038
ENSG00000000938.12 NA
... ...
ENSG00000285979.1 0.600750
ENSG00000285987.1 NA
ENSG00000285990.1 NA
ENSG00000285991.1 0.952921
ENSG00000285994.1 NAsummary(res)
out of 31604 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 2373, 7.5%
LFC < 0 (down) : 1949, 6.2%
outliers [1] : 0, 0%
low counts [2] : 14706, 47%
(mean count < 9)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?resultshead(res_shrunk[, c("log2FoldChange", "lfcSE", "padj")])log2 fold change (MAP): condition Dexamethasone vs Untreated
DataFrame with 6 rows and 3 columns
log2FoldChange lfcSE padj
<numeric> <numeric> <numeric>
ENSG00000000003.14 -0.3469688 0.1027268 5.31137e-03
ENSG00000000419.12 0.1946659 0.1215961 2.93189e-01
ENSG00000000457.13 0.0344560 0.1457579 9.22557e-01
ENSG00000000460.16 -0.0971219 0.2334461 8.72980e-01
ENSG00000000938.12 -0.0361966 0.0791913 NA
ENSG00000000971.15 0.4357043 0.0869471 9.08589e-06这里显式记录 contrast 名称，并同时保留原始估计与收缩后的效应量。lfcShrink() 的目标是让低计数或高方差基因的 log2 fold change 更稳定，不是把结果“变显著”。正式分析应在报告中说明使用的 shrinkage 方法。
MA 图
MA 图把效应量与平均表达放在同一张图上：横轴是标准化后的平均计数，纵轴是 log2 fold change。它有两个用处——一是看效应量是否随表达量变化（低表达端通常更发散），二是直接看出 lfcShrink() 把低计数基因的夸张倍数压回了多少。左右并排画原始与收缩后的结果。
par(mfrow = c(1, 2))
plotMA(res, alpha = 0.01, ylim = c(-6, 6), main = "原始 log2FC")
plotMA(res_shrunk, alpha = 0.01, ylim = c(-6, 6), main = "lfcShrink 之后")
注意 plotMA() 的灰色点是 padj 不显著、蓝色点是显著；alpha 只影响着色，不影响模型。也可以自己用 ggplot2 画同一个东西（横轴建议 scale_x_log10()），但要记得 MA 图的横轴是表达量、纵轴是效应量，和火山图（横轴效应量、纵轴显著性）正好换了轴。
火山图
pload(org.Hs.eg.db)── loading the package: org.Hs.eg.db ──df <- as.data.frame(res)
df$gene_id <- sub("\\.\\w+$", "", rownames(df))
df$symbol <- mapIds(org.Hs.eg.db,
keys = df$gene_id,
column = "SYMBOL",
keytype = "ENSEMBL",
multiVals = "first")'select()' returned 1:many mapping between keys and columnshead(df) baseMean log2FoldChange lfcSE stat
ENSG00000000003.14 739.9407174 -0.36115368 0.10686857 -3.3794190
ENSG00000000419.12 511.7357220 0.20631474 0.12866453 1.6035090
ENSG00000000457.13 314.1948551 0.03783077 0.15863341 0.2384792
ENSG00000000460.16 79.7936215 -0.11525899 0.31499067 -0.3659124
ENSG00000000938.12 0.3072668 -1.36911852 3.50376359 -0.3907565
ENSG00000000971.15 5716.7968219 0.44708633 0.08938729 5.0016766
pvalue padj gene_id symbol
ENSG00000000003.14 7.263921e-04 5.311369e-03 ENSG00000000003 TSPAN6
ENSG00000000419.12 1.088223e-01 2.931887e-01 ENSG00000000419 DPM1
ENSG00000000457.13 8.115095e-01 9.225570e-01 ENSG00000000457 SCYL3
ENSG00000000460.16 7.144304e-01 8.729804e-01 ENSG00000000460 FIRRM
ENSG00000000938.12 6.959772e-01 NA ENSG00000000938 FGR
ENSG00000000971.15 5.683387e-07 9.085892e-06 ENSG00000000971 CFHpload(ivolcano)── loading the package: ivolcano ──ivolcano(df,
logFC_col = "log2FoldChange",
pval_col = "padj",
gene_col = "symbol",
top_n = 5,
onclick_fun=onclick_genecards)基因表达水平
resSig <- subset(res, padj < 0.01)
topGene <- rownames(res)[which.min(res$padj)]
topGene[1] "ENSG00000189221.9"geneCounts <- plotCounts(dds, gene = topGene,
intgroup = c("donor", "condition"),
returnData = TRUE)
geneCounts count donor condition
SRR1039508 445.0772 N61311 Untreated
SRR1039509 4570.3817 N61311 Dexamethasone
SRR1039512 468.9803 N052611 Untreated
SRR1039513 4398.0477 N052611 Dexamethasone
SRR1039516 517.5761 N080611 Untreated
SRR1039517 4607.8726 N080611 Dexamethasone
SRR1039520 253.5867 N061011 Untreated
SRR1039521 3732.9201 N061011 Dexamethasoneggplot(geneCounts, aes(condition, count)) +
geom_boxplot(aes(fill=condition), alpha=.3) +
geom_point(aes(color=donor)) +
geom_line(aes(group=donor, color=donor)) +
scale_y_log10()
差异基因热图
rlog_expr <- assay(rlog(dds))
mat1 <- rlog_expr[order(res$padj)[1:20], ]
pload(dplyr)── loading the package: dplyr ──res2 <- as.data.frame(res)
res2$ID <- rownames(res2)
de <- arrange(res2, padj) |>
group_by(sign(log2FoldChange)) |>
slice(1:20) |>
pull(ID)
mat2 <- rlog_expr[de, ]
anno <- as.data.frame(colData(dds)[, c("donor", "condition")])
plot_heatmap <- function(mat, anno) {
mat <- mat - rowMeans(mat)
p <- pheatmap::pheatmap(mat, annotation_col = anno)
return(p)
}
p1 <- plot_heatmap(mat1, anno)p2 <- plot_heatmap(mat2, anno)aplot::plot_list(p1, p2)
富集分析
ORA：先设阈值，再看名单
df$entrez <- mapIds(org.Hs.eg.db,
keys = df$gene_id,
column = "ENTREZID",
keytype = "ENSEMBL",
multiVals = "first")'select()' returned 1:many mapping between keys and columns# resSig <- subset(df, padj < 0.01)
resSig <- df[order(df$padj, na.last = NA)[1:100], ]
gene_set <- unique(na.omit(resSig$entrez))
universe <- unique(na.omit(df$entrez))
pload(clusterProfiler)── loading the package: clusterProfiler ──x <- enrichGO(
gene = gene_set,
universe = universe,
OrgDb = org.Hs.eg.db,
ont = "BP",
readable = TRUE
)
dotplot(x)
fc <- resSig$log2FoldChange
names(fc) <- resSig$symbol
cnetplot(x, foldChange = fc)
富集分析同时报告 gene_set 和 universe。universe 应代表本次实验中真正有机会被检测、注释并进入检验的基因集合；默认使用全基因组可能系统性放大富集显著性。正式报告还应记录注释包版本和数据库日期。
GSEA：不设阈值，用整张排序表
ORA 必须先把基因切成”显著 / 不显著”，阈值一动结果就换一批；GSEA 换一个问法——把全部基因按效应量排成一张表，看某条通路的基因是不是挤在排序的两端。这样做的好处是那些单个变化不大、但方向一致的通路不会被阈值切掉。
跑 GSEA 之前要先处理两件事：把 ENSEMBL 映射到 ENTREZ 并去重（一个 ENTREZ 只能出现一次，否则排序表不合法），以及按 log2 fold change 从大到小排序。
gl <- df[!is.na(df$log2FoldChange) & !is.na(df$entrez), c("entrez", "log2FoldChange")]
gl <- gl[order(-abs(gl$log2FoldChange)), ]
gl <- gl[!duplicated(gl$entrez), ]
geneList <- sort(setNames(gl$log2FoldChange, gl$entrez), decreasing = TRUE)
length(geneList)[1] 21940head(geneList) 247 100862671 642778 7704 101928819 140711
9.505148 7.363608 7.166855 6.377496 5.923731 5.616437 pload(enrichplot)── loading the package: enrichplot ──gse <- gseGO(geneList = geneList,
OrgDb = org.Hs.eg.db,
keyType = "ENTREZID",
ont = "BP",
minGSSize = 10,
pvalueCutoff = 0.05)Warning in prepare_gsea_inputs(geneList, scoreType, exponent): There are ties
in the preranked stats (1.22% of the list). The order of those tied genes will
be arbitrary, which may produce unexpected results.gse_df <- as.data.frame(gse)
head(gse_df[, c("ID", "Description", "setSize", "NES", "p.adjust")]) ID Description setSize NES p.adjust
GO:0014888 GO:0014888 striated muscle adaptation 30 2.117744 0.02500666
GO:0010273 GO:0010273 detoxification of copper ion 11 1.998013 0.03081380
GO:1990169 GO:1990169 stress response to copper ion 11 1.998013 0.03081380
GO:0014896 GO:0014896 muscle hypertrophy 76 1.955855 0.03081380
GO:0043500 GO:0043500 muscle adaptation 88 1.941661 0.02500666
GO:0003300 GO:0003300 cardiac muscle hypertrophy 72 1.937692 0.03081380dotplot(gse, showCategory = 10, split = ".sign") + facet_grid(. ~ sign(NES))
NES（normalized enrichment score）的符号说明方向：正值表示这条通路的基因在排序表前端富集（在处理组里更高），负值相反。setSize 是通路大小，通路太小时结果不稳定，所以上面设了 minGSSize。
top_id <- gse_df$ID[1]
gseaplot2(gse, geneSetID = top_id, title = gse_df$Description[1])
最后是 leading edge——驱动这条通路信号的那批核心基因。它比整条通路更值得拿去解释：
le <- strsplit(gse_df$core_enrichment[1], "/")[[1]]
le_symbol <- mapIds(org.Hs.eg.db, keys = le, column = "SYMBOL",
keytype = "ENTREZID", multiVals = "first")'select()' returned 1:1 mapping between keys and columnshead(le_symbol, 20) 28999 84676 51778 54206 2308 90523 58529 8837
"KLF15" "TRIM63" "MYOZ2" "ERRFI1" "FOXO1" "MLIP" "MYOZ1" "CFLAR"
2660 5468 51676 8557
"MSTN" "PPARG" "ASB2" "TCAP" 报告 GSEA 结果时必须一起写清：用的是哪张排序表（哪个统计量）、置换次数（nPermSimple，默认 1000）、minGSSize 与 pvalueCutoff、以及注释库版本。和 ORA 一样，GSEA 给的是假说，不是机制。
结果审计任务
结果审计：
- 比较原始
log2FoldChange与res_shrunk的效应量：在 MA 图上挑一个低计数、收缩前后差别最大的基因，说明 shrinkage 为什么会改变它的排序； - 在 PCA、样本相关热图、MA 图、火山图和 top-gene 表中，至少指出一个“统计显著但效应量很小”或“效应量大但不显著”的基因；
- 重新运行一次富集分析，把
universe改成全基因组，记录结果变化并解释背景偏差； - 对比 ORA 与 GSEA 的结果：找出两条路线都支持的通路，以及只在其中一条里出现的通路，并解释它们为什么不一致；
- 取 GSEA 首条通路的 leading edge 基因，说明这批基因与 ORA 的基因名单重叠多少；
- 在报告末尾写清 contrast、design formula、padj 阈值、lfcShrink 方法、GSEA 的排序统计量与置换次数、注释版本和随机性/软件版本。
if (requireNamespace("fanyi", quietly = TRUE)) {
pload(fanyi)
message("翻译扩展需要先配置 fanyi 的 appid/key；未配置时跳过在线翻译。")
# 配置后再取消下面代码的注释：
# set_translate_option(appid = Sys.getenv("FANYI_APPID"),
# key = Sys.getenv("FANYI_KEY"), source = "baidu")
# gs <- gene_summary(df$entrez[1:10])
# gs$summary2 <- translate(gs$summary)
# head(gs)
# translate_ggplot(dotplot(x), axis = "y")
} else {
message("未安装可选包 fanyi，跳过在线翻译扩展。")
}── loading the package: fanyi ──fanyi v0.1.1 Learn more at https://yulab-smu.top/
Please cite:
D Wang, G Chen, L Li, S Wen, Z Xie, X Luo, L Zhan, S Xu, J Li, R
Wang, Q Wang, G Yu. Reducing language barriers, promoting information
absorption, and communication using fanyi. Chinese Medical Journal. 2024,
137(16):1950-1956. doi: 10.1097/CM9.0000000000003242翻译扩展需要先配置 fanyi 的 appid/key；未配置时跳过在线翻译。在线翻译为可选扩展。使用 translate_ggplot 前，请先按 fanyi 包 README 申请并配置百度翻译 appid 与 key，再运行上一个代码块中的配置示例。