一个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
可以看到的是作者给出来的是57.5 Mb 的矩阵文件 :
GSE182923_genes_fpkm_expression.txt.gz
转录组测序数据的表达量矩阵可以有多种格式,每种格式都有其特定的用途和优势。以下是一些常见的格式:
我们通常是针对转录组测序使用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)格式的文件,做后面的差异分析也很难,因为文章自己就一个很垃圾的差异分析结果,如下所示:
很垃圾的差异分析结果
虽然说上面的案例(糖尿病数据集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分析。如果你恰好是先兆子痫研究方向, 就可以把这个数据集更加细致的解读和挖掘一下,未必不是一个课题哦!
如果你确实觉得我的教程对你的科研课题有帮助,让你茅塞顿开,或者说你的课题大量使用我的技能,烦请日后在发表自己的成果的时候,加上一个简短的致谢,如下所示:
We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
十年后我环游世界各地的高校以及科研院所(当然包括中国大陆)的时候,如果有这样的情谊,我会优先见你。