linux コマンドによる htseq-count の結果の結合

HTSeq などのプログラムを利用して得られたカウントデータは、タブ区切りのテキストファイルに保存される。一列目は遺伝子の ID、二列目はカウントデータの値が記載されている。ターミナルで headtail を利用すると、ファイルの最初と最後の数行のデータを見ることができる。

head SRR119866.count.txt
## ENSG00000000003	26
## ENSG00000000005	0
## ENSG00000000419	25
## ENSG00000000457	25
## ENSG00000000460	13
## ENSG00000000938	18
## ENSG00000000971	15
## ENSG00000001036	126
## ENSG00000001084	253
## ENSG00000001167	90

tail SRR119866.count.txt 
## ENSG00000278854	0
## ENSG00000278855	0
## ENSG00000278856	0
## ENSG00000278857	0
## ENSG00000278858	0
## __no_feature	3676532
## __ambiguous	61377
## __too_low_aQual	1454582
## __not_aligned	2240357
## __alignment_not_unique	0

一般に、HTSeq の htseq-count を実行した結果では 1 FASTQ ファイルから 1 つこのようなカウントデータのファイルが得られる。Perl や Python などで簡単なプログラムを作成して、複数のこのようなファイルを一つの行列にまとめることはできるが、ここでは Linux コマンドを利用した方法を紹介する。

join コマンドによる htseq-count の出力結果の結合

2 つのファイル(a.count.txt と b.count.txt)を 1 つのファイル(all.count.txt)に結合させる例。

# カウントデータを結合して一時ファイル(tmp.count.txt)に保存
join -j 1 a.count.txt b.count.txt > tmp.count.txt

# ヘッダー行をつける
echo "GeneID    A    B" > all.count.txt

# カウントデータを all.count.txt に書き入れる
cat tmp.count.txt >> all.count.txt

# 一時ファイルを削除
rm tmp.count.txt

結合するファイルが 3 つ以上の場合は以下のようにする。

join -j 1 a.count.txt b.count.txt | \
    join -j - c.count.txt | \
    join -j - d.count.txt | \
    join -j - e.count.txt > tmp.count.txt

echo "GeneID    A    B    C    D    E" > all.count.txt

cat tmp.count.txt >> all.count.txt

rm tmp.count.txt

awk による htseq-count の出力結果の結合

まず、結合したい複数の htseq-count の出力結果を一つのディレクトリ(counts_dir)に保存する。htseq-count の出力結果が .count.txt で終わるものとする。次に、そのディレクトリに移動し、awk コマンドを実行し、複数のファイルの結合を行う。

cd counts_dir

# ヘッダー行を作成
ls *.count.txt | \
  awk 'BEGIN{FS="__"}{print $1}' | \
  tr '\n' '\t' | ¥
  sed 's/\t$/\n/' > all.count.txt

# htseq-count の出力結果を結合
awk 'NF > 1 {a[$1] = a[$1]"\t"$2}END{for (i in a) print i a[i]}' *.count.txt >> all.count.txt

awk を利用して結合した結果は、遺伝子の並び順が元の htseq-count の順番と異なってくる。そこで、1 列目に対してソートを行えばよい。

sort all.count.txt > all.count.sorted.txt

一部の UNIX 環境では __no_feature や __ambiguous などのように「__」から始まる単語では「__」を無視してソートを行う。そのため、ソート結果が期待する結果と異なるものになったりする。そこで、「__」を含む単語と含まない単語をそれぞれソートしてから、両ファイルを結合するという手法をとるとよい。

cat <(grep -v "__" all.count.txt | sort ) <(grep "^__"  all.count.txt) > all.count.sorted.txt