相信用過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

先給一下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)