「ほっ」と。キャンペーン

<   2013年 12月 ( 4 )   > この月の画像一覧


2013年 12月 20日

マルチファスタをたくさんのシングルファスタに

手元にあるマルチファスタファイルをシングルファスタのたくさんのファイルにします。
今回使うマルチファスタファイルは multi.fasta という名前のファイルです。
それをバラバラにして、遺伝子名.fasta という名前のファイルにします。


#使うものつくる
$ more multi.fasta | grep ">" | tr -d ">" | awk '{print "more multi.fasta | sed -e %/>/s/$/@/g% | tr -d \"\\n\" | tr \"\>\" \"\\n\" | grep "$1" | awk -F \"@\" %{print \"\>\"$1\"\\n\"$2}% > "$1".fasta"}' | tr "%" "'" > m2s.sh
#動くようにする
$ chmod 777 m2s.sh
#動かす
$ sh m2s.sh

使ったもの:more, grep, tr, awk, sed, chmod, sh


#####################################################

#multi.fasta の中を見てみる
$ more multi.fasta
>gene_a
CTAGCTAGTCTAGCTATCTATCGTATCGTAGCTGATGCTAGTCGTAGCTGATCGTAGCTGATGCTGTAGCGATGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTAGTCGATGTGTGTGTGGGGGCGCGCCGGCGCGGCGC
>gene_b
TGACTGATGCTATACGACGATGCTAGCTACGTAGCTAGCTAGCTAGCTGATGAATATATTTATTTCTGCTAGCTATAGCTGATGCTAGTCGTAGTGGGGGCTAGCTAGTCGTAGCTAGTCGTAGTCGTAGCTGATGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
>gene_c
CGATGCTATTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCAAAAAAAAGCTGTAGCTAGCATCTATCTACTATCGTATCACTGACGGCGGCGCGGCGCGCGGGCCGGAGAGCGGACGTCATACGTCAGTACTGTCAGTCAGTACTGCAGTACGTCGTGTCGTAGTCATGACTGTACGTGTAGCACGTCGTACGTACGTCGAGCGAGCTGATTCGATCTGATCTGATTCGTAGCGGCGTGCGCGATCGTAGCTATCGATGCTGATCGAGCTGAACTGGTACTCTGGTAGCTAGCTAGCTTGACGTCATGTACGTCATGTCAGTCAGTACTGTCAGTACGTACTGACTGACGGTCAGTACG

#gene_a のみを切り出してシングルファスタファイルにしたかったらこう
$ more multi.fasta | sed -e '/>/s/$/@/g' | tr -d "\n" | tr ">" "\n" | grep gene_a | awk -F "@" '{print ">"$1"\n"$2}' > gene_a.fasta

#gene_a.fasta をみてみる
$ more gene_a.fasta
>gene_a
CTAGCTAGTCTAGCTATCTATCGTATCGTAGCTGATGCTAGTCGTAGCTGATCGTAGCTGATGCTGTAGCGATGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTAGTCGATGTGTGTGTGGGGGCGCGCCGGCGCGGCGC

#これを全部の遺伝子について書こうとするとこう
$ more multi.fasta | grep ">" | tr -d ">" | awk '{print "more multi.fasta | sed -e %/>/s/$/@/g% | tr -d \"\\n\" | tr \"\>\" \"\\n\" | grep "$1" | awk -F \"@\" %{print \"\>\"$1\"\\n\"$2}% > "$1".fasta"}' | tr "%" "'" | more
more multi.fasta | sed -e '/>/s/$/@/g' | tr -d "\n" | tr ">" "\n" | grep gene_a | awk -F "@" '{print ">"$1"\n"$2}' > gene_a.fasta
more multi.fasta | sed -e '/>/s/$/@/g' | tr -d "\n" | tr ">" "\n" | grep gene_b | awk -F "@" '{print ">"$1"\n"$2}' > gene_b.fasta
more multi.fasta | sed -e '/>/s/$/@/g' | tr -d "\n" | tr ">" "\n" | grep gene_c | awk -F "@" '{print ">"$1"\n"$2}' > gene_c.fasta

#それをm2s.sh という名前で保存します
$ more multi.fasta | grep ">" | tr -d ">" | awk '{print "more multi.fasta | sed -e %/>/s/$/@/g% | tr -d \"\\n\" | tr \"\>\" \"\\n\" | grep "$1" | awk -F \"@\" %{print \"\>\"$1\"\\n\"$2}% > "$1".fasta"}' | tr "%" "'" > m2s.sh

#m2s.sh を実行可能ファイルにします
$ chmod 777 m2s.sh

#ファイルが揃いました
$ ls
gene_a.fastam2s.shmulti.fasta

#さっきつくったgene_a.fasta いらない
$ rm gene_a.fasta

#すっきり
$ ls
m2s.shmulti.fasta

#m2s.sh の実行
$ sh m2s.sh

#実行結果。遺伝子毎にバラバラになってる
$ ls
gene_a.fastagene_b.fastagene_c.fastam2s.shmulti.fasta

#中身を確認
$ head gene_*
==> gene_a.fasta <==
>gene_a
CTAGCTAGTCTAGCTATCTATCGTATCGTAGCTGATGCTAGTCGTAGCTGATCGTAGCTGATGCTGTAGCGATGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTAGTCGATGTGTGTGTGGGGGCGCGCCGGCGCGGCGC

==> gene_b.fasta <==
>gene_b
TGACTGATGCTATACGACGATGCTAGCTACGTAGCTAGCTAGCTAGCTGATGAATATATTTATTTCTGCTAGCTATAGCTGATGCTAGTCGTAGTGGGGGCTAGCTAGTCGTAGCTAGTCGTAGTCGTAGCTGATGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

==> gene_c.fasta <==
>gene_c
CGATGCTATTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCAAAAAAAAGCTGTAGCTAGCATCTATCTACTATCGTATCACTGACGGCGGCGCGGCGCGCGGGCCGGAGAGCGGACGTCATACGTCAGTACTGTCAGTCAGTACTGCAGTACGTCGTGTCGTAGTCATGACTGTACGTGTAGCACGTCGTACGTACGTCGAGCGAGCTGATTCGATCTGATCTGATTCGTAGCGGCGTGCGCGATCGTAGCTATCGATGCTGATCGAGCTGAACTGGTACTCTGGTAGCTAGCTAGCTTGACGTCATGTACGTCATGTCAGTCAGTACTGTCAGTACGTACTGACTGACGGTCAGTACG


実際のファスタファイルにはスペースが入っていたりすることが多いので、事前に殺すか-F "\t" で回避。
元のマルチファスタと同じところにシングルファスタがもりもり入るとうざいから事前に
mkdir single
として箱つくって
> "$1".fasta" のところを > single/"$1".fasta" にするとすっきり



人気ブログランキングに参加しています。
応援よろしくお願いします。
FC2 Blog Ranking

[PR]

by koretoki | 2013-12-20 14:24
2013年 12月 18日

エクセルなメタボロームデータを化合物IDで揃えたテーブルにしたい

#同じ値に複数の化合物が含まれているデータのモデルデータ
$ more test2
C00001 ID_1 ID_2 ID_3 2.1E-04 10.1E-04 0.4E-04 N.D. N.D. 1.2E-04
C011010 ID_1 ID_2 ID_3 1.11 6.4E-04 6.3E-04 N.D. N.D. N.D.
C001983,C03231,C04213 ID_1 ID_2 ID_3 N.D. 2.1E-05 1.1E-02 0.9E-02 1.1E-02 1.1E-02
C034143 ID_1 ID_2 ID_3 6.8E-10 8.1E-99 N.D. N.D. N.D. N.D.
C044555,C022233 ID_1 ID_2 ID_3 4.8 2.1E-04 4.8E-04 1.3E-04 N.D. N.D.


#ひとまずいらないIDを消して化合物毎にわけます。
$ more test2 | awk '{print $1"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}' | tr "," "\n" | awk '{print NR"\t"$1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7}' | sort -nr
8 C022233 4.8 2.1E-04 4.8E-04 1.3E-04 N.D. N.D.
7 C044555
6 C034143 6.8E-10 8.1E-99 N.D. N.D. N.D. N.D.
5 C04213 N.D. 2.1E-05 1.1E-02 0.9E-02 1.1E-02 1.1E-02
4 C03231
3 C001983
2 C011010 1.11 6.4E-04 6.3E-04 N.D. N.D. N.D.
1 C00001 2.1E-04 10.1E-04 0.4E-04 N.D. N.D. 1.2E-04


#データが空になっちゃった化合物のところにデータを戻します。
$ more test2 | awk '{print $1"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}' | tr "," "\n" | awk '{print NR"\t"$1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7}' | sort -nr | awk '(NR==1){a=$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8; print; next}{if($3<0){print $1"\t"$2"\t"a} if($3>0){print; a=$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8} }' | more
8 C022233 4.8 2.1E-04 4.8E-04 1.3E-04 N.D. N.D.
7 C044555 4.8 2.1E-04 4.8E-04 1.3E-04 N.D. N.D.
6 C034143 6.8E-10 8.1E-99 N.D. N.D. N.D. N.D.
5 C04213 N.D. 2.1E-05 1.1E-02 0.9E-02 1.1E-02 1.1E-02
4 C03231 N.D. 2.1E-05 1.1E-02 0.9E-02 1.1E-02 1.1E-02
3 C001983 N.D. 2.1E-05 1.1E-02 0.9E-02 1.1E-02 1.1E-02
2 C011010 1.11 6.4E-04 6.3E-04 N.D. N.D. N.D.
1 C00001 2.1E-04 10.1E-04 0.4E-04 N.D. N.D. 1.2E-04


#前の行を参照する方法
awk '
(NR==1){a=$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8; print; next} #1行目は絶対入ってるのでaに値を代入してそのまま吐いて次の行にいく
{
if($3<0){print $1"\t"$2"\t"a} #行のデータが空だったらaを吐く
if($3>0){print; a=$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8} #行のデータが空じゃなかったらそのまま吐いて自分の行の値をa に代入する
}' | more


セル内改行は手と目で殺した
[PR]

by koretoki | 2013-12-18 13:31
2013年 12月 12日

gen. pre.

genome.fasta->geneprediction.gff->codingseq.fasta

export AUGUSTUS_CONFIG_PATH=~/tools/augustus/augustus-2.4/config/

~/tools/augustus/augustus-2.4/src/augustus --cds=on --codingseq=on --species=homo_n0rr homo_n0rr.fasta > homo_n0rr.gff

more homo_n0rr.gff | grep "#" | tr -d "#" | tr -d " " | sed -e 's/startgene/startgene@/g' -e 's/endgene/endgene@/g' | tr -d "\n" | tr "@" "\n" | grep endgene | tr "]" "[" | sed -e 's/coding/[coding/g' | tr atcg ATCG | awk -F "[" '{print ">augustus_hn_"$1"\n"$3}' > codingseq_homo_n0rr.fasta

#$5でpepseq

人気ブログランキングに参加しています。
応援よろしくお願いします。
FC2 Blog Ranking

[PR]

by koretoki | 2013-12-12 01:13
2013年 12月 06日

NLTKで論文を読んでみた

今日の材料はこれ
Transcriptome Responses of Insect Fat Body Cells to Tissue Culture Environment

映画「The Social Network」の脚本をNLTKで解析して遊んでみたや、
「魔法少女まどか☆マギカ」の台詞をNLTK(Natural Language Toolkit)で解析するに習ってやってみました。
途中グラフ描くのがやばかったので、Matpotlibのインストールで詰んだあたなへを参考にEPDを使って楽しました。

論文のテキストをコピーペーストし、'pc論文.txt' として保存しました。
iPythonを起動し、

#20131206追記
import nltk
from nltk.book import *

#データの読み込み
raw=o p e n ('pc論文.txt').read()
#単語認識
tokens=nltk.word_tokenize(raw)
text=nltk.Text(tokens)

#単語数
len(tokens)
Out[12]: 5791

#語彙数
tokens_1 = [w.lower() for w in tokens] #20131206修正('s/I/1/g')
len(set(tokens_1))
Out[14]: 1935

訳2000単語知ってれば論文が書ける!!

#何回'cell' って言ってるか
tokens_1.count('cell')
Out[15]: 36

tokens_1.count('cells')
Out[16]: 27

tokens_1.count('fat')
Out[17]: 81

tokens_1.count('body')
Out[18]: 70

bodies がいる予感...

#culture とprimary が含まれる文
text.concordance("culture" and "primary",lines=5)
Displaying 5 of 8 matches:
es. Cell lines are established from primary culture of tissue when a population
roliferating cells derived from the primary tissue explant undergo immortalizat
undergo immortalization [ 5 ] . In primary cultures of insect cells , it usual
ught that during the early stage of primary culture , isolated explants activat
ant are different from cells in the primary tissue. Champy ( 1913 ) proposed th

#単語の分布をみる
fdist=nltk.FreqDist(w.lower() for w in text)
fdist.plot(50,cumulative=True)

e0160319_6553397.png


よく出てくる単語の出てくるタイミングがみたくなりますね。

#単語を指定
terms=['cells','cell','transcriptome','genes','culture','intact']

#指定した単語の分布を見る
text.dispersion_plot(terms)

transcriptomeについて話してからgenesに移行するのが多いらしい
e0160319_6564779.png


人気ブログランキングに参加しています。
応援よろしくお願いします。
FC2 Blog Ranking

[PR]

by koretoki | 2013-12-06 10:20