ixxmu
commented
9 months ago 规范统一格式的GEO RNA-seq count及其标准化数据 by 生信星球
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
参考网址：https://www.ncbi.nlm.nih.gov/geo/info/rnaseqcounts.html#why
0.痛点
GEO网站的芯片数据是有统一格式的，都放在series.matrix.txt.gz文件里面了。但是RNA-seq数据，series.matrix.txt.gz文件里只能拿到临床信息（表型信息），而表达矩阵却是空的。大多数的作者在supplementary file里面提供了表达矩阵，格式不受任何限制，xls格式、csv、tsv、txt都有，最麻烦的是那种每个样本一个文件，需要自行读取并合并成矩阵的。例如：还有的数据没有提供差异分析所需使用的count，而是提供了rpkm、tpm、fpkm等标准化后的数据，这些数据不适合用作差异分析，且不能再逆转回count，这就使得GEO的RNA-seq数据的前期整理过程比较难搞。
1.统一格式的转录组数据来啦
GEO是NCBI旗下的网站，为了提高存量RNA-seq数据的利用率，SRA 和 GEO 团队搭建了一套流程，统一计算 RNA-seq 基因表达矩阵，方便进行差异分析和可视化。
2.范围
提交给 GEO 的人类和小鼠 RNA-seq 数据，人类数据已经可用，小鼠数据官方说2023年秋季上线，但今天（2023.12.27）目前暂时还不能用，估计快了。原有数据已经可用，新提交的数据大约需要等待一周时间
3.可用文件
count、rpkm、tpm、基因注释表格
4.流程
HISAT2→Subread featureCounts
5.如何获取
举个栗子，这是GSE204753的页面，其他数据改编号即可https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=GSE204753https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=GSE164073&format=file&file=GSE164073_raw_counts_GRCh38.p13_NCBI.tsv.gz.
在上面这个链接里替换GSE编号即可。不大喜欢多次复制粘贴改来改去，所以我在2.3.2以上版本的tinyarray中加入了一个函数：get_count_txt，目前这个版本只在github上。预计要2024年1月中旬才能出现在cran上啦。下载方式
devtools::install_github("xjsun1221/tinyarray") #2024年1月中旬以前用这个
install.packages("tinyarray")#2024年1月中旬以后用这个
使用方式
library(tinyarray)
get_count_txt("GSE204753")
## https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=GSE204753&format=file&file=GSE204753_raw_counts_GRCh38.p13_NCBI.tsv.gz
出现的网址就是上面那个网址替换你的编号的结果啦。这样就不用总是复制粘贴了。
我还顺便写了一个转换行名为symbol的函数嘞。是用bitr函数转换的id。
dat = read.delim("GSE204753_raw_counts_GRCh38.p13_NCBI.tsv.gz",
row.names = 1)
dat[1:4,1:4]
## GSM6190621 GSM6190623 GSM6190624 GSM6190625
## 100287102 6 5 7 7
## 653635 10 12 18 17
## 102466751 0 0 0 0
## 107985730 0 0 0 0
exp = trans_entrezexp(dat)
## 'select()' returned 1:1 mapping between keys and columns
## Warning in clusterProfiler::bitr(rownames(entrezexp), fromType = "ENTREZID", :
## 4.17% of input gene IDs are fail to map...
## 37734 rownames transformed after duplicate rows removed
exp[1:4,1:4]
## GSM6190621 GSM6190623 GSM6190624 GSM6190625
## DDX11L1 6 5 7 7
## WASH7P 10 12 18 17
## MIR6859-1 0 0 0 0
## MIR1302-2HG 0 0 0 0
6.注意事项
(1) NCBI生成的表达矩阵未必与数据提交者上传的数据完全一致
因NCBI流程与数据提交者使用的软件、参数的不同，有所差别属于正常现象。
(2) 质控
The only parameter a run must pass for inclusion in the NCBI pipeline is that the run is of type 'transcriptomic' and it has a genome alignment rate over 50%.
type必须是转录组，且比对率超过50%
(3) 检查样本是否具有可比性
提交者经常在同一研究中存放多种类型的数据（例如，RNA-seq和RIP-seq），这意味着即使在同一个矩阵中，RNA counts也无法直接比较。或者，尽管样本属于同一类型，但它们可能仍然不用于比较。需要查看原始记录，以确定研究中的所有样本是否都打算直接比较。
(4)不同研究的count比较时要谨慎
尽管所有NCBIcount数据都是由同一流程生成的，但如果试图比较不同研究的count，应特别小心。NCBI 流程没有矫正不同实验室的差别或者其他混杂因素。
(5)标准化矩阵文件可能未充分标准化
当总 RNA 含量及其分布差异很大时，FPKM （RPKM） 和 TPM 不应用于样品间的定量比较 。正如Zhao等人所讨论的，“一个常见的误解是RPKM和TPM值已经标准化，因此应该在样本或RNA-seq项目中具有可比性。然而，RPKM 和 TPM 代表测序总体中转录本的相对丰度，因此取决于样品中 RNA 总体的组成。通常，可以合理地假设总RNA浓度和分布在比较样品中非常接近。然而，在不同的实验条件下和/或不同测序方案中，测序的RNA库可能存在显着差异;因此，在这种情况下，基因表达的比例不能直接比较”。
(6)样本缺失
没达到50% 比对率的样本会被删掉，或者有的样本会因技术原因处理失败。
如果因为代码看不懂，而跟不上正文的节奏，可以来找我，相当于来一个新手保护期。我的课程都是循环开课。下一期的时间，点进去咨询微信咯 生信分析直播课程 生信新手保护学习小组
https://mp.weixin.qq.com/s/Wq0vlCfzIVoRr0BmZ4dZvg