多个表达矩阵文件合并

今天群主给了我们学徒一个任务,下载数据集:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84073

我想看到,HCC,CHC,CC这3组,跟healthy的分开比较,然后出3个火山图,3个热图。那么首先需要下载counts值矩阵,样本信息如下:

GSM2653819 Healthy1-Tissue_1 [RNA-Seq]
GSM2653820 Healthy1-Tissue_2 [RNA-Seq]
GSM2653821 Healthy2-Tissue [RNA-Seq]
GSM2653822 Healthy3-Tissue [RNA-Seq]
GSM2653823 HCC1-Tissue_1 [RNA-Seq]
GSM2653824 HCC1-Tissue_2 [RNA-Seq]
GSM2653825 HCC3-Tissue [RNA-Seq]
GSM2653826 CHC1-Tissue_1 [RNA-Seq]
GSM2653827 CHC1-Tissue_2 [RNA-Seq]
GSM2653828 CHC2-Tissue [RNA-Seq]
GSM2653829 CC1-Tissue_1 [RNA-Seq]
GSM2653830 CC1-Tissue_2 [RNA-Seq]
GSM2653831 CC2-Tissue [RNA-Seq]
GSM2653832 CC3-Tissue [RNA-Seq]

肯定是不能一个个手动点击样本信息进入寻找文件下载链接,那样低效。

查看具体的每个文件

压缩包解压的方式下载表达矩阵后,发现,每个样本都是一个文本文件

GSM2653819_Counts_notmergedTR_Healthy1_Tissue_1.txt.gz
GSM2653820_Counts_notmergedTR_Healthy1_Tissue_2.txt.gz
GSM2653821_Counts_notmergedTR_Healthy2_Tissue.txt.gz
GSM2653822_Counts_notmergedTR_Healthy3_Tissue.txt.gz
GSM2653823_Counts_notmergedTR_HCC1_Tissue_1.txt.gz
GSM2653824_Counts_notmergedTR_HCC1_Tissue_2.txt.gz
GSM2653825_Counts_notmergedTR_HCC3_Tissue.txt.gz
GSM2653826_Counts_notmergedTR_CHC1_Tissue_1.txt.gz
GSM2653827_Counts_notmergedTR_CHC1_Tissue_2.txt.gz
GSM2653828_Counts_notmergedTR_CHC2_Tissue.txt.gz
GSM2653829_Counts_notmergedTR_CC1_Tissue_1.txt.gz
GSM2653830_Counts_notmergedTR_CC1_Tissue_2.txt.gz
GSM2653831_Counts_notmergedTR_CC2_Tissue.txt.gz
GSM2653832_Counts_notmergedTR_CC3_Tissue.txt.gz
GSM2653833_Counts_notmergedTR_Healthy1_Organoid_1.txt.gz
GSM2653834_Counts_notmergedTR_Healthy1_Organoid_2.txt.gz
GSM2653835_Counts_notmergedTR_Healthy2_Organoid.txt.gz
GSM2653836_Counts_notmergedTR_Healthy2_Organoid_DM.txt.gz
GSM2653837_Counts_notmergedTR_Healthy3_Organoid.txt.gz
GSM2653838_Counts_notmergedTR_Healthy3_Organoid_DM.txt.gz
GSM2653839_Counts_notmergedTR_HCC1_Organoid_1.txt.gz
GSM2653840_Counts_notmergedTR_HCC1_Organoid_2.txt.gz
GSM2653841_Counts_notmergedTR_HCC3_Organoid.txt.gz
GSM2653842_Counts_notmergedTR_CHC1_Organoid_1.txt.gz
GSM2653843_Counts_notmergedTR_CHC1_Organoid_2.txt.gz
GSM2653844_Counts_notmergedTR_CHC1_Organoid_a.txt.gz
GSM2653845_Counts_notmergedTR_CHC1_Organoid_b.txt.gz
GSM2653846_Counts_notmergedTR_CHC2_Organoid.txt.gz
GSM2653847_Counts_notmergedTR_CC1_Organoid_1.txt.gz
GSM2653848_Counts_notmergedTR_CC1_Organoid_2.txt.gz
GSM2653849_Counts_notmergedTR_CC1_Organoid_a.txt.gz
GSM2653850_Counts_notmergedTR_CC1_Organoid_b.txt.gz
GSM2653851_Counts_notmergedTR_CC1_Organoid_c.txt.gz
GSM2653852_Counts_notmergedTR_CC2_Organoid.txt.gz
GSM2653853_Counts_notmergedTR_CC3_Organoid.txt.gz

格式很统一,如下:

Ensembl Symbol Counts
ENSG00000278566 OR4F29 0
ENSG00000273547 OR4F16 0
ENSG00000187634 SAMD11 33
ENSG00000188976 NOC2L 160
ENSG00000187961 KLHL17 13
ENSG00000187583 PLEKHN1 0
ENSG00000187642 PERM1 0
ENSG00000188290 HES4 1
ENSG00000187608 ISG15 5
ENSG00000188157 AGRN 187
ENSG00000237330 RNF223 5
ENSG00000131591 C1orf159 0
ENSG00000162571 TTLL10 8

现在就需要批量依次读取这些文件,然后合并成为表达矩阵!

首先参考群主的WGCNA教程的合并方法

当时群主的代码是linux的shell脚本+R里面的dcast函数,如果大家感兴趣群主的WGCNA教程,见:

我仔细看了看代码,其实就是首先在linux是把多个文件合并成为 tmp.txt 文本。

 ## https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48213
 #wget -c ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE48nnn/GSE48213/suppl/GSE48213_RAW.tar
 #tar -xf GSE48213_RAW.tar
 #gzip -d *.gz
 ## 首先在GSE48213_RAW目录里面生成tmp.txt文件,使用shell脚本
 # awk '{print FILENAME"\t"$0}' GSM*.txt |grep -v EnsEMBL_Gene_ID >tmp.txt
 # 其实也可以直接使用R来读取GSE48213_RAW.tar里面的gz文件,这里就不演示了
 # 可以参考:https://mp.weixin.qq.com/s/OLc9QmfN0YcT548VAYgOPA 里面的教程
 ## 然后把tmp.txt导入R语言里面用reshape2处理即可
 # 这个 tmp.txt 文件应该是100M左右大小哦。

这个文本有点特殊,其实就是把每个txt文件夹,按照行的方式首尾连接起来成为一个大文本,但是第一列加上了样本信息!

然后在R里面读取后,使用reshape2包的dcast函数即可,如下所示,一句话搞定!

a=read.table('GSE48213_RAW/tmp.txt',sep = '\t',stringsAsFactors = F)
 library(reshape2)
 fpkm <- dcast(a,formula = V2~V1)

上面的方法当然是可行的,但是依赖于linux环境,在mac下面稍微有点不一样,在Windows就需要借助于git等软件来使用shell脚本。我猜想应该是那个WGCNA教程已经是四年前的啦,当时群主的主要编程语言并不是R,所以这样的文本合并需求,会采取LINUX+R的方式搞定!

第二种方法是lapply循环读取文件

这个是纯粹的R语言解决方案,我也是在群主的指点下完成的,可以看到里面使用了 do.calllapply 函数批量读取txt文本文件:

rm(list = ls())
options(stringsAsFactors = F)
fs=list.files('countsFiles/')
a=do.call(cbind,lapply(list.files('countsFiles/'), function(x){
 read.table(file.path('countsFiles',x),
 header = T,sep = '\t',row.names = 1)[,2]
}))
a[1:4,1:4]

# 下面是合并后的表达矩阵添加行名和列名
rownames(a)=rownames(read.table('countsFiles/GSM2653820_Counts_notmergedTR_Healthy1_Tissue_2.txt.gz',
 header = T,sep = '\t',row.names = 1))
colnames(a)=gsub('.txt.gz','',substring(fs,nchar("GSM2653853_Counts_notmergedTR_")+1,1000))
a[1:4,1:4]

我不知道什么样的函数叫做优雅,但是看起来这个就有点高大上!

第3种方法你来写吧

反正数据集就是GSE84073,进入就看到了可以下载的txt文件,自行摸索合并!

最后当然是拿到表达矩阵后做差异分析了

这个群主的教程已经足够多了,走标准分析流程,火山图,热图,GO/KEGG数据库注释等等。这些流程的视频教程都在B站和GitHub上,目录如下:

  • 第一讲:GEO,表达芯片与R
  • 第二讲:从GEO下载数据得到表达量矩阵
  • 第三讲:对表达量矩阵用GSEA软件做分析
  • 第四讲:根据分组信息做差异分析
  • 第五讲:对差异基因结果做GO/KEGG超几何分布检验富集分析
  • 第六讲:指定基因分组boxplot指定基因list画热图

仅仅是最后得到的差异分子,并不是以前的mRNA后面的基因名,而是miRNA,lncRNA,甚至circRNA的ID,看起来很陌生罢了。感兴趣可以细读表达芯片的公共数据库挖掘系列推文 ;

如果表达矩阵需要注释探针

也可以看群主在2019年的尾巴推出3个R包:

  • 第一个是整合全部的bioconductor里面的芯片探针注释包。
  • 第二个是整合全部GPL的soft文件里面的芯片探针注释包。
  • 第三个是下载全部的GPL的soft文件里面的探针碱基序列比对后注释包。

配合着详细的介绍:

因为这些包暂时托管在GitHub平台,但是非常多的朋友访问GitHub困难,尤其是我打包了好几百个GPL平台的注释信息后, 我的GitHub包变得非常臃肿,大家下载安装困难,所以我重新写一个精简包。也在:芯片探针ID的基因注释以前很麻烦 和 :芯片探针序列的基因注释已经无需你自己亲自做了, 里面详细介绍了。最重要的是idmap函数,安装方法说到过:芯片探针序列的基因注释已经无需你自己亲自做了, 使用起来也非常简单:

library(AnnoProbe)
ids=idmap('GPL570',type = 'soft')
head(ids)

仅仅是一句话,就拿到了这个平台的探针的注释信息。需要注意的是,这个函数的type参数,其实是有3个选择,这里我演示的是选择soft这个来源的基因注释信息。

并不是所有的平台都是有soft注释,也不是所有的平台都被我的这个工具囊括哦。

Comments are closed.