こんにちは，h--ishiです．
夜はすっかり肌寒い季節になってきましたね．
日中との寒暖差で体調を崩さないようにご自愛ください．
今日は，bedtools intersectの出力結果がおかしくなった話をしたいと思います．
「bedtools intersectとはなんぞや？」という方のために，弊社ブログ過去記事リンクを貼っておきます．
staffblog.amelieff.jp
一言でいうと，gff, bed, vcfなどの二つのファイルを比べてオーバーラップしている場所を抽出したり，
あるファイル内の各ゲノム領域内において，もう片方のファイルに記述のある領域（例えば，変異など）がどれくらいあるかを数えたりできるというツールです．
先日，bedtools intersectを使っていて，出力結果の列が元のファイルの列と融合するという現象に遭遇しました．
bedtoolsのバージョンは，2.27.1でした．
これからその時の再現を行いたいと思います．
テストデータセット
- unicorn_no_genes.gtf: ユニコーンの遺伝子情報が記載されているgtfファイル
- unicornA_no_henni.vcf: ユニコーン（個体A）の変異情報が記載されているvcfファイル
- unicornB_no_henni.vcf: ユニコーン（個体B）の変異情報が記載されているvcfファイル
unicorn_no_genes.gtfの中身．
ユニコーンの遺伝子が二つあり，それぞれ1番染色体の1-250 bp，251-500 bpにあるようです．
chr1 RefSeq region 1 500 . + . ID=Uni_001.1 chr1 Gnomon gene 1 250 . - . ID=gene-Unicorn_no_idenshi1 chr1 Gnomon gene 251 500 . - . ID=gene_Unicorn_no_idenshi2
unicornA_no_henni.vcfの中身
こちらは，ユニコーンAの変異情報が書かれているファイルです．1番染色体の50 bp，255 bp，472 bpに一塩基置換があるようです．
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT chr1 50 . A T 476 . . . chr1 255 . C G 409 . . . chr1 472 . T C 601 . . .
unicornB_no_henni.vcfの中身
こちらは，ユニコーンBの変異情報が書かれているファイルです．1番染色体の50 bp，255 bpに一塩基置換があるようです．
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT chr1 50 . A T 476 . . . chr1 255 . C G 409 . . .
・やろうとしたこと
unicorn_no_genes.gtfファイル中の各行（feature領域）に，unicornA_no_henni.vcf, unicornB_no_henni.vcfに記載のある変異がそれぞれいくつ存在するかをbedtools intersect を用いて集計しようとした．gffは9列ありますが，10列目にユニコーンAの変異数，11列目にユニコーンBの変異数が記載されているタブ区切りテキストファイルが出力されるのを期待していました．
・やったこと
私は，
まず，unicornA_no_henni.vcf 中の変異がunicorn_no_genes.gtf中の領域にいくつあるかを数えて結果output.txt を作成しました．
実行コマンドはこちらです．
bedtools intersect -header -c -a unicorn_no_genes.gtf -b unicornA_no_henni.vcf > output.txt
そうすると，10列目に集計結果が追加されます．
次に，このoutput.txt ファイルはgffフォーマット+10列目に集計結果というファイルなので，ほとんどgffフォーマットです．
そこで，output.txtを入力ファイルとして，11列目にunicornB_no_henni.vcfについての集計結果を吐き出してくれるだろうと横着をしました．
実行コマンドです．
bedtools intersect -c -header -a output.txt -b unicornB_no_henni.vcf > output2.txt
・結果
output.txt の中身はこちらです．
ユニコーンAは，1番染色体の50 bp，255 bp，472 bpに一塩基置換がありました．
そのため，Unicorn_no_idenshi1 には1箇所変異があり，Unicorn_no_idenshi2 には2箇所変異があるという結果が10列目に出力されました．
##gff-version 3 #!gff-spec-version 1.21 #!processor NCBI annotwriter #!genome-build Unicorn_2.0 #!genome-build-accession NCBI_Assembly:Uni_001.1 ### chr1 RefSeq region 1 500 . + . ID=Uni_001.1 3 chr1 Gnomon gene 1 250 . - . ID=gene-Unicorn_no_idenshi1 1 chr1 Gnomon gene 251 500 . - . ID=gene_Unicorn_no_idenshi2 2
私が期待したoutput2.txt はgffファイルの10列目にunicornA_no_henni.vcf，11列目にunicornB_no_henni.vcfの集計結果が追加されているファイル．
しかし，
できあがったものは...
##gff-version 3 #!gff-spec-version 1.21 #!processor NCBI annotwriter #!genome-build Unicorn_2.0 #!genome-build-accession NCBI_Assembly:Uni_001.1 ### chr1 RefSeq region 1 500 . + . ID=Uni_001.13 2 chr1 Gnomon gene 1 250 . - . ID=gene-Unicorn_no_idenshi11 1 chr1 Gnomon gene 251 500 . - . ID=gene_Unicorn_no_idenshi22 1
おわかりいただけただろうか？（杉本○みさん風の声で）
10列目のデータが9列目と融合してしまっています．
2回目のintersectの入力ファイルが10列（9列目まではgffフォーマット）になっていて，
そのファイルを読み込むとき，あるいは出力ファイルを出力する時に不具合が起きているのかなと想像しています．
（bedtools intersect のコードを見ているわけではないので，想像の域を出ません）
・まとめ
横着せず，プログラムが指定するちゃんとしたgff, vcf, bedフォーマットのファイルを入力ファイルとして使用しましょう（自戒）．
