生物信息学技能面试题(第4题)-多个同样的行列式文件合并起来

相信用过htseq-count的朋友都知道，它是分开对每个样本计算所有的基因表达量，所以会生成一个个独立的文件，我用perl脚本模仿它的结果如下：

$ head a.txt gene_1 178 gene_2 692 gene_3 486 gene_4 666 gene_5 395 gene_6 48 gene_7 926 gene_8 733 gene_9 660 gene_10 578

第一列是基因，第二列是该基因的counts值，共有a~z这26个样本的counts文件，需要合并成一个大的行列式，这样才能导入到R里面做差异分析，如果手工用excel表格做，当然是可以的，但是太麻烦，如果有500个样本，正常人都不会去手工做了，需要编程。 生成测试文件的代码如下：

#首先新建文件tmp.sh 输入这个代码：

perl -le '{print "gene_$_\t".int(rand(1000)) foreach 1..99}'

## 然后用perl脚本调用这个tmp.sh文件：

perl -e 'system(" bash tmp.sh >$_.txt") foreach a..z'

##这样就生成了a~z这26个样本的counts文件

用shell或者perl或者python，设置R语言都可以做，但是各有优缺点，而且如果每个样本的基因顺序并不一致，这时候你应该怎么做呢？ 实际需求如下：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48213 里面有56个文件(ftp://ftp.ncbi.nlm.nih.gov/geo/s ... pl/GSE48213_RAW.tar)，需要合并成一个表达矩阵，来根据cell-line的不同，分组做差异分析。 paper是：https://www.ncbi.nlm.nih.gov/pubmed/24176112 输出的表达矩阵，如下所示：

先给一下shell结合R语言的做法：

## 首先在GSE48213_RAW目录里面生成tmp.txt文件

awk '{print FILENAME"\t"$0}' * |grep -v EnsEMBL_Gene_ID >tmp.txt

## 然后把tmp.txt导入R语言里面用reshape2处理即可！

setwd('tmp/GSE48213_RAW/')

a=read.table('tmp.txt',sep = '\t',stringsAsFactors = F)

library(reshape2)

fpkm <- dcast(a,formula = V2~V1)