多个表达矩阵文件合并
前些天群主给了我们学徒一个任务,下载数据集：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 <- 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，看起来很陌生罢了。感兴趣可以细读表达芯片的公共数据库挖掘系列推文 ；
- 解读GEO数据存放规律及下载，一文就够
- 解读SRA数据库规律一文就够
- 从GEO数据库下载得到表达矩阵 一文就够
- GSEA分析一文就够（单机版+R语言版）
- 根据分组信息做差异分析- 这个一文不够的
- 差异分析得到的结果注释一文就够
如果表达矩阵需要注释探针
也可以看群主在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注释，也不是所有的平台都被我的这个工具囊括哦。