作者仅提供了fpkm格式表达量矩阵的转录组测序数据集该如何重新分析呢

一个2021的糖尿病转录组文章:《Altered human alveolar bone gene expression in type 2 diabetes—A cross-sectional study》,在线链接是:https://onlinelibrary.wiley.com/doi/abs/10.1111/jre.12947

研究者们在GEO数据库是有数据分享:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE182923

  • RNA samples were further purified using ribosomal RNA depletion technique and processed for RNA sequencing and analysis.
  • Expression levels for mRNAs were performed by calculating FPKM ([total_exon_fragments/mapped reads (millions) × exon length (kB)]),
  • 最后是差异分析后确定统计学阈值: differentially expressed mRNAs were selected with log2 (fold change) >1 or log2 (fold change) ≤1 and with a parametric F test comparing nested linear models.

可以看到的是作者给出来的是57.5 Mb 的矩阵文件 :

GSE182923_genes_fpkm_expression.txt.gz

转录组测序数据的表达量矩阵可以有多种格式,每种格式都有其特定的用途和优势。以下是一些常见的格式:

  1. 计数矩阵(Count Matrix)
    • 这是最基本的格式,通常由比对到参考基因组的读段生成。
    • 每一行代表一个基因或转录本,每一列代表一个样本。
    • 单元格中的值表示该基因在该样本中的读段计数。
  2. FPKM/FPKM-UQ(每千个碱基每百万片段的比率/未量化的FPKM)
    • FPKM是标准化的表达量指标,考虑了基因长度和测序深度。
    • FPKM-UQ是未量化的FPKM,它没有经过标准化处理,通常用于避免引入人为的表达量变化。
  3. TPM(每千个转录本每百万片段的比率)
    • TPM是另一种标准化的表达量指标,它考虑了样本中的总转录本数。
    • TPM使得不同样本间的基因表达量可比。
  4. CPM(每百万计数的比率)
    • CPM是一种简单的标准化方法,将计数除以样本的总计数乘以百万。
    • 它用于归一化数据,使得不同样本间的表达量可比。
  5. RSEM/Cufflinks输出
    • RSEM(RNA-Seq by Expectation-Maximization)和Cufflinks是两种软件工具,它们提供了一种估计基因和转录本表达量的方法。
    • 输出通常包括每个基因的估计表达量(如FPKM)、表达量的不确定性和统计评估。
  6. Salmon输出
    • Salmon是一种用于RNA-Seq数据的无需比对的定量工具,它使用轻量级比对和EM算法来估计表达量。
    • 输出通常包括每个转录本的TPM和预期计数(expected count)。

我们通常是针对转录组测序使用DESeq2/edgeR进行差异分析,而DESeq2/edgeR要求的输入数据是计数矩阵(raw Count Matrix)格式,作者并没有提供,而且我们不可能依据作者提供的FPKM矩阵去反推出来原始的计数矩阵(raw Count Matrix)。

这里我们推荐:https://www.ncbi.nlm.nih.gov/geo/geo2r/?acc=GSE182923

而且这个geo2r网页工具还贴心的给出来了代码,如下所示:

# Version info: R 4.2.2, Biobase 2.58.0, GEOquery 2.66.0, limma 3.54.0
################################################################
# Data plots for selected GEO samples

# load counts table from GEO
urld <- "https://www.ncbi.nlm.nih.gov/geo/download/?format=file&type=rnaseq_counts"
path <- paste(urld, "acc=GSE182923", "file=GSE182923_raw_counts_GRCh38.p13_NCBI.tsv.gz", sep="&");
tbl <- as.matrix(data.table::fread(path, header=T, colClasses="integer"), rownames=1)

# pre-filter low count genes
# keep genes with at least 2 counts > 10
keep <- rowSums( tbl >= 10 ) >= 2
tbl <- tbl[keep, ]

# log transform raw counts
# instead of raw counts can display vst(as.matrix(tbl)) i.e. variance stabilized counts
dat <- log10(tbl + 1)

# box-and-whisker plot
par(mar=c(7,4,2,1))
boxplot(dat, boxwex=0.7, notch=T, main="GSE182923", ylab="lg(cnt + 1)", outline=F, las=2)

# UMAP plot (dimensionality reduction)
library(umap)
dat <- dat[!duplicated(dat), ] # first remove duplicates
ump <- umap(t(dat), n_neighbors = 5, random_state = 123)
plot(ump$layout, main="GSE182923 UMAP plot, nbrs =5", xlab="", ylab="", pch=20, cex=1.5)
library(car)
pointLabel(ump$layout, labels = rownames(ump$layout), method="SANN", cex=0.6)

值得注意的是表达量矩阵并不是每个样品会被纳入

可以看到的是作者虽然是 Cross-sectional, no replicates, 10 healthy, 8 diabetic,但是geo2r仅仅是纳入了4个疾病和7个对照哦,我推测应该是这个数据集的测序质量很差,有一些样品不满足前面的转录组定量要求就被暴力删除了,其实也是合理的选择样品 :

不满足前面的转录组定量要求就被暴力删除了

当然了,就算是我们拿到了DESeq2/edgeR要求的输入数据是计数矩阵(raw Count Matrix)格式的文件,做后面的差异分析也很难,因为文章自己就一个很垃圾的差异分析结果,如下所示:

很垃圾的差异分析结果

GEO数据库的任意转录组测序数据集均可获得count矩阵

虽然说上面的案例(糖尿病数据集GSE182923)是因为作者自己的原因导致我们虽然是获得count矩阵但是差异分析结果也丑爆了。但是这个解决方案是 通用的, 理论上GEO数据库的任意转录组测序数据集均可获得count矩阵。比如这个GSE148241-先兆子痫-数据集,是 placentae from 9 patients with early-onset severe preeclampsia (EOSPE) and 32 normal controls, 同样的方式获取count矩阵和分组信息 :

# load counts table from GEO
urld <- "https://www.ncbi.nlm.nih.gov/geo/download/?format=file&type=rnaseq_counts"
path <- paste(urld, "acc=GSE148241", "file=GSE148241_raw_counts_GRCh38.p13_NCBI.tsv.gz", sep="&");
path
tbl <- as.matrix(data.table::fread(path, header=T, colClasses="integer"), rownames=1)

# data<-data.table::fread("matrix.txt",
# data.table = F)
dim(tbl)
mat=as.data.frame(tbl) 
library(AnnoProbe)
library(GEOquery) 
gset = getGEO("GSE148241", destdir = '.', getGPL = F) 
pd = pData(gset[[1]]) 
table(group_list)
group_list=ifelse(grepl('Normal',pd$source_name_ch1),
 'control','case' )

就可以常规的差异分析,如下所示的火山图和热图:

火山图和热图

可以看到有两个样品是离群点, 其实这个GSE148241-先兆子痫-数据集页面也指出来了,但是数据集配套的文献并没有关心这个差异分析结果,反而是做了一个wgcna分析。如果你恰好是先兆子痫研究方向, 就可以把这个数据集更加细致的解读和挖掘一下,未必不是一个课题哦!

Comments are closed.