kuroの覚え書き

96の個人的覚え書き

multi fastaファイルを1遺伝子ごとのファイルに分割するには

multi fastaファイルを1個ずつのfastaに分割したい。
まずはfastaのseq部分の改行をなくす

$ awk -v ORS= '/^>/ { $0 = (NR==1 ? "" : RS) $0 RS } END { printf RS }1' fasta.txt > fasta_awk.txt

次にfastaを2行ごとに分割。多数のファイルが同じ階層にできるので、対象のmultifastaファイルと共に新しいフォルダに入れておく。

$ split -l 2 fasta_awk.txt gene_

これでgene_aaのようなファイルが出力されるので

#!/bin/bash

for file in gene_*; do
  firstline=$(cat "$file" | head -1 | sed -e 's/\//-/g')
  filename=${firstline#>}
  [ -n "$filename" ] && [ ! -e "${filename}.txt" ] &&
  mv "$file" "${filename}.fasta"
done

こういうスクリプトで1行目の>{trasscript_ID}から>を除いたIDを取得して、{transcript_ID}.fastaというファイル名に置き換えてやる。


これを全部まとめて対象ファイルを引数にするスクリプト

#!/bin/sh

fpath=$1
echo "input file: ${fpath}"
faname=`basename $1`
fdir="${fpath%/*}"
fname="${faname%.*}"
fext="${faname#*.}"
opath1="${fdir}/${fname}_1.${fext}"
odir="${fdir}/fastas"
opath="${odir}/gene_"
mkdir ${odir}
awk -v ORS= '/^>/ { $0 = (NR==1 ? "" : RS) $0 RS } END { printf RS }1' ${fpath} > ${opath1}
split -l 2 ${opath1} ${opath}

for file in ${odir}/gene_*; do
  firstline=$(cat "$file" | head -1 | sed -e 's/\//-/g')
  filename=${firstline#>}
  [ -n "$filename" ] && [ ! -e "${filename}.txt" ] &&
  mv "$file" "${odir}/${filename}.${fext}"
done

こんな感じ。