bedファイルの分割

魚種ごとのミトコンドリアCDS領域の位置情報は

NC_000860_Salvelinus_fontinalis 2843    3817
NC_000860_Salvelinus_fontinalis 4032    5080
NC_000860_Salvelinus_fontinalis 5472    7022
NC_000860_Salvelinus_fontinalis 7186    7876
NC_000860_Salvelinus_fontinalis 7952    8119
NC_000860_Salvelinus_fontinalis 8110    8792

という形式で保存した。しかし、それぞれの魚種のミトコンドリアCDSの配列を取ってきたいのでこのbed形式のファイルをそれぞれの魚で分けて保存したほうが都合がいい。そこで次のシェルスクリプトを書いた。（シェルスクリプトで書く必要はあまりなさそう）

#!/bin/bash
 
cat ~/work/mitCysFish/annotation/mitoCDSStartEnd | cut -f1 | sort | uniq > ~/work/mitCysFish/CDSseq/tmp/speciesName.tmp
cd ~/work/mitCysFish/CDSseq/tmp/
 
mkdir bed_raw
 
cat speciesName.tmp | while read line
do
    cat ~/work/mitCysFish/annotation/mitoCDSStartEnd | awk -v species=$line 'match($1,species){print $0}' > bed_raw/$line.bed
done

ここでかなり引っかかったが、seqkitを用いてシーケンスを抜き出す際、クエリーとしてのbedファイルの1列目に記載されたシーケンス名が検索対象のfastaのヘッダーと一致していなければならない。そこで各魚のbedに対応するfastaのヘッダーを抜き出してきてbedの1列目を書き換える。また、どの魚か分かるように後ろの列に魚種情報を足しておく。

#!/bin/bash
 
num=1
 
cd ~/work/mitCysFish/CDSseq/tmp/
mkdir bed_edit
 
for i in bed_raw/*
do
    name=$(cat speciesName.tmp | awk -v rowN=$num 'NR==rowN{print}')
    queryName=$(cat ~/work/mitCysFish/mitoGenome/raw/$name.fa | grep ">gi" | sed -e 's/>//g')
    spTmp=$(cat <(echo ${i##*/} | sed -e 's/.bed//g'))
    ((num++)) #インクリメント
    cat $i | awk -v q="$queryName" -v sp="$spTmp" 'BEGIN{OFS="\t"}{print q,$2,$3,sp,NR}' > bed_edit/"$spTmp"_edit.bed 
    #シェル変数をawkに導入する際、""で囲った方が無難な模様。特に今回は魚種名にスペースが入っていたためエラーが起きた。
done

編集したbedを用いて配列を抜き出す。

#!/bin/bash
 
cd ~/work/mitCysFish/CDSseq
 
for i in tmp/bed_edit/NC*
    do
        nameTmp=$(echo ${i##*/} | sed -e 's/_edit.bed//g')
        seqkit subseq --id-ncbi --bed $i ~/work/mitCysFish/mitoGenome/raw/"$nameTmp".fa
    done > mitCDSSeq.fa

bed形式は座標を0始まりにするとか少しややこしいので、抜き出した配列がフレームシフトしてないか心配。そこでblastxでアミノ酸レベルで相同性が取れているか確認する。

[kijima.yusuke@m48 CDSseq]$ for i in ~/files/m48/kijima.yusuke/work/mitCysFish/CDSseq/BLAST/splitFasta/mitCDSSeq.fa_chunk_0000*; do runGE-1cpu-4gb "blastx -query $i -db ~/files/m48/kijima.yusuke/work/DBforBLAST/mitoCDS/mito.aa -outfmt 6 -evalue 1e-50 -max_target_seqs 1 -out ${i##*/}.blast"; done
 
[kijima.yusuke@m48 splitBLAST]$ cat *.blast | wc -l
29355
[kijima.yusuke@m48 splitBLAST]$ cat ../splitFasta/mitCDSSeq.fa_chunk_0000* | grep ">" | wc -l
36894

クエリーの7~8割は相同性が取れてるから一応大丈夫かな。一応transeqでアミノ酸配列に変換してからもう少しだけ確認する。 ↓ 追記：blastxってframe3まではずらしてくれるんですね。アミノ酸に変換したらBLASTが当たらなくなったのでもう一度確認。

BED形式について調べたところこの記事が分かりやすかった→https://qiita.com/suimye/items/e4d067d04d891462207e 要するに開始位置の座標のみ-1すればよい。というわけで

[kijima.yusuke@m48 BLAST]$ for i in ~/work/mitCysFish/CDSseq/tmp/bed_edit/NC_0*; do cat $i | awk -F'\t' 'BEGIN{OFS="\t"}{print $1,$2-1,$3,$4,$5}'>${i##*/};done