差分

このページの2つのバージョン間の差分を表示します。

--- cds領域からfastaファイルを作成する [2018/05/30 05:32] – [BED形式について] 133.11.222.89
+++ cds領域からfastaファイルを作成する [Unknown date] (現在) – 削除 - 外部編集 (Unknown date) 127.0.0.1
@@ 行 1: / 行 1: @@
-====== bedファイルの分割 ======
-魚種ごとのミトコンドリアCDS領域の位置情報は
-<code>
-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
-</code>
-という形式で保存した。しかし、それぞれの魚種のミトコンドリアCDSの配列を取ってきたいのでこのbed形式のファイルをそれぞれの魚で分けて保存したほうが都合がいい。そこで次のシェルスクリプトを書いた。（シェルスクリプトで書く必要はあまりなさそう）
-<code bash>
-#!/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
-</code>
-====== 分割したbedの編集 ======
-ここでかなり引っかかったが、''seqkit''を用いてシーケンスを抜き出す際、クエリーとしてのbedファイルの1列目に記載されたシーケンス名が検索対象のfastaのヘッダーと一致していなければならない。そこで各魚のbedに対応するfastaのヘッダーを抜き出してきてbedの1列目を書き換える。また、どの魚か分かるように後ろの列に魚種情報を足しておく。
-<code bash>
-#!/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
-</code>
-====== 配列抜き出し ======
-編集したbedを用いて配列を抜き出す。
-<code bash>
-#!/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
-</code>
-====== BLASTでフレームがずれてないか確認 ======
-bed形式は座標を0始まりにするとか少しややこしいので、抜き出した配列がフレームシフトしてないか心配。そこでblastxでアミノ酸レベルで相同性が取れているか確認する。
-<code bash>
-[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
-[kijima.yusuke@m48 splitBLAST]$ cat ../splitFasta/mitCDSSeq.fa_chunk_0000* | grep ">" | wc -l
-</code>
-クエリーの7~8割は相同性が取れてるから一応大丈夫かな。一応transeqでアミノ酸配列に変換してからもう少しだけ確認する。
-↓
-追記：blastxってframe3まではずらしてくれるんですね。アミノ酸に変換したらBLASTが当たらなくなったのでもう一度確認。
-====== BED形式について ======
-BED形式について調べたところこの記事が分かりやすかった→[[https://qiita.com/suimye/items/e4d067d04d891462207e]]
-{{:bed_intro.png?600|}}
-要するに開始位置の座標のみ-1すればよい。
-というわけで
-<code bash>
-[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
-</code>
-これでseqkitで切り出すスクリプトを少し改変してfastaが入手できた。