c++ 合并2个txt_多个表达矩阵文件合并

前些天群主给了我们学徒一个任务,下载数据集：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教程，见：

• 一文看懂WGCNA 分析(2019更新版)
• 通过WGCNA作者的测试数据来学习
• 重复一篇WGCNA分析的文章(代码版)
• 重复一篇WGCNA分析的文章(解读版)(逆向收费读文献2019-19)
• 关键问题答疑：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 

上面的方法当然是可行的，但是依赖于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，看起来很陌生罢了。感兴趣可以细读表达芯片的公共数据库挖掘系列推文 ；

• 解读GEO数据存放规律及下载，一文就够

• 解读SRA数据库规律一文就够

• 从GEO数据库下载得到表达矩阵 一文就够

• GSEA分析一文就够(单机版+R语言版)

• 根据分组信息做差异分析- 这个一文不够的

• 差异分析得到的结果注释一文就够

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

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

• 第一个是整合全部的bioconductor里面的芯片探针注释包。

• 第二个是整合全部GPL的soft文件里面的芯片探针注释包。

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

配合着详细的介绍：

• 第三个万能芯片探针ID注释平台R包

• 第二个万能芯片探针ID注释平台R包

• 第一个万能芯片探针ID注释平台R包

• GEO数据库中国区镜像横空出世

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

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

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

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