カテゴリ:未分類( 99 )


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
2013年 06月 25日

ホルマリン固定パラフィン包埋なサンプルからmRNAの発現が見れるかどうかは配列による

生検用の肝臓のサンプルはホルマリン固定してパラフィン包埋するそうなのですが、それからRNAを取って発現を見ようというお話があるそうです。
試薬も売ってますね→ホルマリン固定パラフィン包埋組織からのRNA抽出
http://www.roche-biochem.jp/biochemica/2006/10/newrna.html

ホルマリンによる固定は塩基にモノメチロールを付加させることで変質させてしまうことが知られています。
著者らは、ホルマリン固定パラフィン包埋サンプルから得られる発現データについて検証するために、培養細胞を用いて実験を行いました。

J Clin Pathol. 2013 Jun 11.
Development of a robust protocol for gene expression analysis using formalin-fixed, paraffin-embedded liver transplant biopsy specimens.
Thompson E, Burt AD, Barker CE, Kirby JA, Brain JG.


実験には胆管上皮細胞(BECs)を用いました。培養にH2O2を添加することでストレスを与え、実験区と比較区のそれぞれについて、p21, TGFβ1, TGFβ2 の3つの遺伝子の発現量を測定しました。
次に、同様に培養した胆管上皮細胞をホルマリン固定後パラフィン包埋したのちにRNAを抽出して同様に発現量を調べました。

結果、p21 とTGFβ1 については、固定や包埋を行わなかった区と、固定後包埋した区で再現性が取れましたが、TGFβ2では再現性が取れませんでした。
固定後包埋したサンプルのmRNA発現を見たい時は、事前にin vitro の系でいいからそのmRNA配列が固定包埋に耐えるか調べた方がいいよ、と著者らは述べています。


どのくらいの期間持つんだろうね。


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

[PR]

by koretoki | 2013-06-25 23:26
2013年 06月 24日

awkとsedで塩基配列をアミノ酸配列にする

手元のマルチファスタな塩基配列をアミノ酸配列に変えたい。
でも、色んなパッケージを取ってくるのはよくわかんないし依存関係とか考えてるより培養したい...

awk とかsed とか誰のmac にでも入ってるものだけでなんとかしてみました。
どんな環境でも動きますけど死ぬほど遅いです。
ひとつの配列を処理するのに1秒ぐらいかかります。

使う時は、#########に囲われてる部分をコピーして"nucl2prot.sh" という名前で保存。
手元のマルチファスタな塩基配列を"transcripts.fasta" という名前にして

sh nucl2prot.sh transcripts.fasta proteins.fasta

と打てば、"proteins.fasta" という名前のマルチファスタなアミノ酸配列が手に入ります。

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

#nucl2prot.sh

#塩基配列をアミノ酸に。コドンテーブルを目的に合わせて書き直せ
#数字の3桁"," 区切りを利用しているので合わない時は考えて
#すごく遅い
#$1=[input.file] $2=[output.file]

more $1 | grep -v ">" | tr "ATCGN" "12345" | sed -e :a -e 's/\(.*[0-9]\)\([0-9]\{3\}\)/\1,\2/;ta' | tr "12345" "ATCGN" | awk '{print $1","}' | sed -e 's/TCA,/S/g' -e 's/TCC,/S/g' -e 's/TCG,/S/g' -e 's/TCT,/S/g' -e 's/TTC,/F/g' -e 's/TTT,/F/g' -e 's/TTA,/L/g' -e 's/TTG,/L/g' -e 's/TAC,/Y/g' -e 's/TAT,/Y/g' -e 's/TAA,/*/g' -e 's/TAG,/*/g' -e 's/TGC,/C/g' -e 's/TGT,/C/g' -e 's/TGA,/*/g' -e 's/TGG,/W/g' -e 's/CTA,/L/g' -e 's/CTC,/L/g' -e 's/CTG,/L/g' -e 's/CTT,/L/g' -e 's/CCA,/P/g' -e 's/CCC,/P/g' -e 's/CCG,/P/g' -e 's/CCT,/P/g' -e 's/CAC,/H/g' -e 's/CAT,/H/g' -e 's/CAA,/Q/g' -e 's/CAG,/Q/g' -e 's/CGA,/R/g' -e 's/CGC,/R/g' -e 's/CGG,/R/g' -e 's/CGT,/R/g' -e 's/ATA,/I/g' -e 's/ATC,/I/g' -e 's/ATT,/I/g' -e 's/ATG,/M/g' -e 's/ACA,/T/g' -e 's/ACC,/T/g' -e 's/ACG,/T/g' -e 's/ACT,/T/g' -e 's/AAC,/N/g' -e 's/AAT,/N/g' -e 's/AAA,/K/g' -e 's/AAG,/K/g' -e 's/AGC,/S/g' -e 's/AGT,/S/g' -e 's/AGA,/R/g' -e 's/AGG,/R/g' -e 's/GTA,/V/g' -e 's/GTC,/V/g' -e 's/GTG,/V/g' -e 's/GTT,/V/g' -e 's/GCA,/A/g' -e 's/GCC,/A/g' -e 's/GCG,/A/g' -e 's/GCT,/A/g' -e 's/GAC,/D/g' -e 's/GAT,/D/g' -e 's/GAA,/E/g' -e 's/GAG,/E/g' -e 's/GGA,/G/g' -e 's/GGC,/G/g' -e 's/GGG,/G/g' -e 's/GGT,/G/g' -e 's/---,/-/g' -e 's/NNN,/N/g' -e 's/N..,/N/g' -e 's/.N.,/N/g' -e 's/..N,/N/g' > temp_nucl2prot_seq.txt

more $1 | grep ">" > temp_nucl2prot_head.txt

paste temp_nucl2prot_head.txt temp_nucl2prot_seq.txt | awk '{print $1"\n"$2}' > $2

rm temp_nucl2prot_head.txt
rm temp_nucl2prot_seq.txt

echo "(^q^)/"


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

以上です。
ほんとに遅いから心配になったら

wc -l temp_nucl2prot_seq.txt

とかやって進行状況を見ても安心して下さい。

#20140109追記
塩基配列の文字数が3の倍数になってないとダメです。
その場合は生成されたアミノ酸配列に「,」が入ってきます。
数字の桁区切りの「,」を利用して3文字を区切っているため、例えば4文字であれば
ATCG→1234→1,234→A,S
と残念な結果になります。
5文字だと、
AATCG→11234→11,234→AA,S

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

[PR]

by koretoki | 2013-06-24 21:31
2013年 06月 23日

ゲノムアノテーションファイルがエクセルだった

いただいたゲノムアノテーションファイルがエクセルでした。
genome.faとそれでtranscripts.faなファイルを作りたいと思います。
以下作業記録

# カラムとして管理していく。IDは付けない。
# まずエクセルから欲しそうなところを取り出してhenkan.txt にした。

# henkan.txt はこういうファイルでタブ区切り。エクセルから取ってきた。
# $1は遺伝子の名前、$2はアミノ酸配列、$3はどっちの鎖に付けられるか、$4ポジションのカラムでスタートとエンドが..で繋がってる。ポジションで複数のエキソンから成る場合は半角スペースで区切られてる。$5は全体のスタート位置、$6は全体のエンド位置を示す。$7はエキソンの数、$7は$4と同じ。$8はアノテーションが付いてる
Gene#Size (aa)StrandPositionsstartendexon#PositionsDescription
scaffold00001.g1.t1306-297..121429712141297..1214-
scaffold00001.g10.t1427+28319..295992831929599128319..29599-
scaffold00001.g100.t1364-224230..2253212242302253211224230..225321-
scaffold00001.g1000.t1337-2042669..20436792042669204367912042669..2043679-
scaffold00001.g1001.t11174+2043981..20475022043981204750212043981..2047502hypothetical protein hogehoge
scaffold00001.g1002.t11106+2047818..2048726 2048852..2049681 2049836..20514142047818205141432047818..2048726 2048852..2049681 2049836..2051414genehogehoge
scaffold00001.g173.t1914-371907..373954 374592..375274 375331..3753413719073753413371907..373954 374592..375274 375331..375341Putative protein of unknown function [hoge hoge]
scaffold00001.g18.t1645-43949..44175 44371..44416 44513..461744394946174343949..44175 44371..44416 44513..46174hogehoge

# 変換元のデータから場所を示すデータを切り出す。
more henkan.txt | grep scaffold | awk -F "\t" '{print $4}' | tr " " "." | awk -F "." '{print $1".."$3"."$4+0".."$6+0"."$7+0".."$9+0"."$10+0".."$12+0"."$13+0".."$15+0"."$16+0".."$18+0"."$19+0".."$21+0}' > kiridashi.txt


# kiridashi.txt の中身。
# 最大で7つのエキソンがあるので、startとendを7セット指定している(s1)..(e1).(s2)..(e2).(s3)..(e3).(s4)..(e4).(s5)..(e5).(s6)..(e6).(s7)..(e7)
head -n 2 kiridashi.txt
297..1214.0..0.0..0.0..0.0..0.0..0.0..0
28319..29599.0..0.0..0.0..0.0..0.0..0.0..0

# 何番染色体かを示す
head -n 2 scaffold_number.txt
1
1

# 上記の2ファイルを結合させて、stard end の指定と染色体番号を示すファイルを作る。
paste kiridashi.txt scaffold_number.txt | tr "\t" "." > temp002.txt

head -n 2 temp002.txt
297..1214.0..0.0..0.0..0.0..0.0..0.0..0.1
28319..29599.0..0.0..0.0..0.0..0.0..0.0..0.1

# ゲノムファイルを変換して一行のファイルにしてある。もちろん疑似配列な
more temp003.txt | awk '{print substr($1,1,50)}'
scaffold00001@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00002@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00003@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00004@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00005@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00006@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00007@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00008@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00010@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00011@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00012@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00013@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00014@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...


# 上のファイルを使って、こういうコマンドラインをつくりたい。
more temp003.txt | grep 1 | awk -F "@" '{print substr($2,297,918)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)}' >> transcriptSEQ.txt
more temp003.txt | grep 1 | awk -F "@" '{print substr($2,28319,1281)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)}' >> transcriptSEQ.txt
more temp003.txt | grep 1 | awk -F "@" '{print substr($2,224230,1092)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)}' >> transcriptSEQ.txt

# 上記のコマンドラインを作るためのコマンドライン
# 植木算的な意味で+1してる。エキソンが7未満の奴はawk のところだと($2,0,1)になってしまうので、sedでごまかした。
more temp002.txt | awk -F "." '{print "more temp003.txt | grep "$22" | awk -F \"@\" _{print substr($2,"$1","$3-$1+1")substr($2,"$4","$6-$4+1")substr($2,"$7","$9-$7+1")substr($2,"$10","$12-$10+1")substr($2,"$13","$15-$13+1")substr($2,"$16","$18-$16+1")substr($2,"$19","$21-$19+1")}_ >> transcriptSEQ.txt"}' | tr "_" "'" | sed -e 's/($2,0,1)/($2,0,0)/g' > make_transcriptsSEQ.sh

#-鎖用のファイルを作るよ
tr "ATCG" "TAGC" < transcriptSEQ.txt | rev > minus_transcriptSEQ.txt

#プラス鎖かマイナス鎖かを指定するファイル
more henkan.txt | grep scaffold | awk -F "\t" '{print $3}' > plus_or_minus.txt

#中身
more henkan.txt | grep scaffold | awk -F "\t" '{print $3}' | head -n 3
-
+
-

# プラスとマイナスを見てどっちを使うのか指定する。
paste plus_or_minus.txt transcriptSEQ.txt minus_transcriptSEQ.txt | awk '{if ($1=="+"){print $2} else if ($1=="-"){print $3}}' > transcriptSEQUENCE.txt

#遺伝子の名前取ってくる
more henkan.txt | grep scaffold | awk -F "\t" '{print $1}' | tr "." "_" > gene_names.txt

#合体させて欲しかったtranscripts.fasta をつくる
paste gene_names.txt transcriptSEQUENCE.txt | awk '{print ">"$1"\n"$2}' > transcripts.fasta


これで完成です。
ちゃんときれいにやる方法を勉強しないとなーと思い続けて2年目...



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

[PR]

by koretoki | 2013-06-23 18:20
2013年 06月 15日

内部寄生蜂のvenom をLC-MSとRNA-seqで調べたよ

内部寄生蜂というのは、宿主の体内に住んで育つ蜂のグループのことです。
当然宿主の体内に入ってれば、宿主の免疫系から「おまえだれ?」って聞かれたりするわけで、侵入してるのがバレたらボコボコにされてしまいます。

内部寄生蜂の親は子どもが宿主の体内でボコられないようにvenom (なんか"毒"とは直訳しない雰囲気ある)を宿主の体内に打ち込んで、宿主の血球によう包囲を抑制しようとします。

キイロショウジョウバエを宿主として寄生する蜂がいて、キイロショウジョウバエはモデル生物で情報が多いので、キイロショウジョウバエ側の免疫の研究は結構されてました。でも、寄生蜂の方はほとんど情報がなかったので、今回著者らはvenom がどんなタンパク質なのかを調べました。

PLoS One. 2013 May 23;8(5):e64125. doi: 10.1371/journal.pone.0064125. Print 2013.
Integrative approach reveals composition of endoparasitoid wasp venoms.
Goecks J, Mortimer NT, Mobley JA, Bowersock GJ, Taylor J, Schlenke TA.


実験では2種類の寄生蜂、L. b. とL. h. を使いました。
RNA-seq で寄生蜂腹部で発現してる遺伝子を特定し、LC-MS で特定したvenomu 線の中身のペプチドをそこにはりつけることで、venom のタンパク質を特定しようとしています。

RNA-seq から作ったトランスクリプトデータに、mass で得たアミノ酸配列を貼付けたところ、L.b. で129のvenom 遺伝子候補、L.h. で176のvenom 遺伝子候補がそれぞれとれました。

それらの中から、雌蜂特異的に発現するもの(雄vs雌のRT-PCR比較から)、
venom 腺特異的に発現するもの(雌体vsvenom 腺のRT-PCR比較から)を絞り込むことにしました。
※他の蜂のvenomとして、Calreticulin や、Heat-shock protein 70 が知られているので、上記の方法ではそれらの遺伝子のような、他の組織で他の仕事をしているvenom 遺伝子を取りこぼすリスクがあるそうです。

ランダムに10個のvenom 遺伝子候補を選んでRT-PCRしたところ、
L.b. でもL.h. でも10個中の8個は雌特異的に発現をしていました。

L.b. ではそれら8個の遺伝子はvenom 腺以外のところでも発現していました。L.b. では、ショウジョウバエ寄生蜂Asobara japonica と同様に、venom 遺伝子が卵巣など他の組織でも発現しているのでしょう。
一方で、L.h. では8個の遺伝子はvenom 腺特異的な発現をしていました。


次に、既知のvenom 遺伝子と相同性のある遺伝子を集めて、分泌シグナルついてるやつだけを集めました。
GO解析をしたところ、L.b. のサンプルでは、antioxidant activity がエンリッチしていました。
ショウジョウバエの免疫応答では、酸素ラジカルを発生させていると考えられているので、L.b. は抗酸化タンパクをvenom として用いて、免疫に対抗していると考えられます。



RNA-seq とMASS を組み合わせた方法は、ゲノムが分かっていない生物の分泌タンパクを特定するのにすごい効くと著者らは述べています。




寄生蜂の話をするといっつもvenom っていう謎物質が出てきて苦しめられてきたのですが、蜂によって全然違うんですね。
今回の研究で見つかった200個くらいのvenom 候補をひとつづつ試してもらって、venom の成分完全同定されるようになると期待しています。
あと、いままでGO階席でぱりっと意味をとる論文に出会えてなかったので、すごい絞り込んだ遺伝子群でG0解析→エンリッチみっけ→エンリッチさせている遺伝子これな→知見と会わせてばっちり考察してるのが楽しかったです。





以下メモ

L.b., L. h. について
幼虫か蛹に産卵
L. b. は卵を免疫から隠して、血球が卵の表面にふれるのを防ぐ。宿主の免疫は寄生されていても強い。
L. h. は宿主の免疫を破壊する。

最初はcDNAをクローニングして進めた。4つとれた。でも、100くらいはありそうなので、ほんの表面だけしか見れていない。

P. puparum のmass もやったけどトランスクリプトがないから56個の予測venom のうち12個しか取れなかった。

Chelonus inanitus のESTsやった。29個みつかって、いくつかは完全にユニークな配列だった。

Microctonus hyperodae でRNA-seqやった。8個の全長venom 遺伝子見つかった。

Nasonia vitripennis でmass やってゲノムのにはっつけたら79venom 遺伝子見つかった。

→核酸+アミノ酸のアプローチがいい
→NGSでせめるのがいい



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

[PR]

by koretoki | 2013-06-15 10:30
2013年 06月 13日

ショウジョウバエの脂肪酸合成酵素群の配列をまるっと欲しい場合

特定の機能に関連した遺伝子の配列がまるっと欲しい気持ちになることがしばしばあります。
そういうときはKEGGを使うのがいいよって教えていただいたのでときどき使っています。
例えば、ショウジョウバエの脂肪酸合成酵素群をまるっと欲しいときは、

KEGG pathway maps のメニュー(以下リンク)を開く
http://www.genome.jp/kegg-bin/get_htext?query=00061&htext=br08901.keg&option=-a

脂肪酸合成パスウェイが見たいので"fatty acid biosynthesis"をクリックして開く(以下リンク)
http://www.genome.jp/kegg-bin/show_pathway?map00061

"Reference patheway" のところをクリックしてショウジョウバエの選んでG0をクリック(以下リンク)
http://www.genome.jp/kegg-bin/show_pathway?org_name=dme&mapno=00061&mapscale=&show_description=hide

脂肪酸合成パスウェイのうち、ショウジョウバエで見つかっている物は緑色に塗られています。
緑に塗られているところをクリックするとそれらの配列を見ることができるので、それをポチポチコピペしていけば脂肪酸合成パスウェイ上の配列を集めることができます。



でもやりたくないです。
30個くらいやるのに1時間かかったことあるし。
twitter で聞いたら@nonshu さんが教えて下さったのでまとめました。


マップの番号が00061であることは分かっています。(マップの左下に書いてある)
キイロショウジョウバエの略称はdmeなので(マップ上のどれか適当に選んでクリックして配列のファスタのヘッダーに書いてあるのがそれ。例えば、">dme:Dmel_CG11198 ACC; Acetyl-CoA carboxylase (EC:6.3.4.14 6.4.1.2); K11262 acetyl-CoA carboxylase / biotin carboxylase [EC:6.4.1.2 6.3.4.14] (A) "。)

基本URL1[http://www.genome.jp/dbget-bin/get_linkdb?-t+genes+path:]のおしりに[dme00061]をくっつけてURLを書きます。

例)http://www.genome.jp/dbget-bin/get_linkdb?-t+genes+path:dme00061

すると、IDとDefinition がでてくるのでIDだけコピペして+で繋げます。
例)dme:Dmel_CG11198+dme:Dmel_CG12170+dme:Dmel_CG17374+dme:Dmel_CG3523+dme:Dmel_CG7842

それを、基本URL2[http://www.genome.jp/dbget-bin/www_bget?-f+-n+a+]のおしりにくっつける。

例)http://www.genome.jp/dbget-bin/www_bget?-f+-n+a+dme:Dmel_CG11198+dme:Dmel_CG12170+dme:Dmel_CG17374+dme:Dmel_CG3523+dme:Dmel_CG7842

リンク先に欲しい配列がファスタフォーマットで待ってます。

もし塩基で欲しい時には、aをnに変えます。
基本URL3[http://www.genome.jp/dbget-bin/www_bget?-f+-n+n+]のおしりにIDをくっつける。

例)http://www.genome.jp/dbget-bin/www_bget?-f+-n+n+dme:Dmel_CG11198+dme:Dmel_CG12170+dme:Dmel_CG17374+dme:Dmel_CG3523+dme:Dmel_CG7842

リンク先に欲しい配列がファスタ形式で待ってます。


はやく迎えに行ってあげて





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

[PR]

by koretoki | 2013-06-13 19:52
2013年 02月 14日

経費な本リスト

1 2013/12/06 数量化革命 アルフレッドWクロスビー amazon
2 2013/12/12 ベイズ統計データ解析 姜興起 amazon
[PR]

by koretoki | 2013-02-14 17:46