多个表达矩阵文件合并

Posted on by

今天群主给了我们学徒一个任务，下载数据集：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.call 和 lapply 函数批量读取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.