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

タグ:バイオインフォマティクス ( 8 ) タグの人気記事


2016年 12月 14日

R作図 軸を省略する

軸を途中波線で省略したい時の作図
データは下記の通り。data.txtとする。

pbentropyhyst
09.9582011
010.391451
010.550411
0.2510.177021
0.2510.36021
0.2510.434361
18.1718641
17.5698931
18.0024861
2.57.6817491
2.57.3344651
2.58.0597491
12.58.0045861
12.57.5191571
12.57.9314761
010.128374
09.7516554
010.264994
0.259.3499724
0.2510.220824
0.259.1107974

data<-read.table("data.txt",header=T)
attach(data)
options(repos="http://cran.md.tsukuba.ac.jp")
library(plotrix)
plot(log(pb+0.09),entropy,xaxt="n",xlab="mM (Drug concentration)",ylab="Hij (transcriptome diversity)",pch=ifelse(hyst<2,1,3))
axis(side=1,at=c(log(0.09),log(0.25+0.09),log(1+0.09),log(2.5+0.09),log(12.5+0.09)),labels=c("0","0.25","1.0","2.5","12.5")) #軸に好きな点を数字を書き足してる
axis.break(1,log(0.15)) #線に切り込みを入れてる

e0160319_14401994.png

http://dx.doi.org/10.1371/journal.pone.0144822.g001

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

[PR]

by koretoki | 2016-12-14 14:41
2016年 06月 24日

培養改善を目指した開発を行う前に確認しておきたいトランスクリプトーム

数十mlのフラスコでの培養の結果と、数Lのジャーでの培養の結果が一致しないことはよくあります。というか、ほとんど一致しないと思います。
その原因をアレイで調べ、培養方法を変えてフラスコとジャーの結果を一致させた研究です。

PMID: 25271019
Hypoxia influences protein transport and epigenetic repression of CHO cell cultures in shake flasks

著者らの先行研究で、フラスコ培養において銅の添加が細胞増殖や組み替えタンパク質生産量が改善できるというものがありましたが、ジャーで再現することができませんでした。
フラスコとジャーが一致しなかった原因を調べるため、それぞれで培養したCHO細胞をマイクロアレイで分析しました。
分析の結果、
タンパク質の細胞内輸送に関わる遺伝子の発現量がフラスコ培養<ジャー培養、
エピジェネティック制御に関わる遺伝子の発現量がフラスコ培養>ジャー培養、
となっていた他、低酸素状態で発現する遺伝子の発現がフラスコ培養>ジャー培養となっていました。

そこで、先の銅添加による培養改善効果をジャーで再現するため、ジャーでの培養を元々のDO60%からDO15%に下げて実験を行いました。
DOを15%まで下げた結果、金属イオン(銅とか鉄とか、と論文中にはフワッとした記述)による培養改善効果がジャーで再現したとのことです。

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


培養改善のための開発を行う前に、実生産タンク培養中のCHO細胞のデータと開発時の評価系に用いる培養中のCHO細胞のデータを比べておいて、ここならこの系で開発しても実生産タンクに持っていけるなっていうのを確認しようと思いました。
[PR]

by koretoki | 2016-06-24 17:15
2016年 06月 18日

医薬品製造用CHOにも使えるゲノム編集技術:RNA-seqから見つけた遺伝子が実際に良い遺伝子であった例

「RNA-seqから見つけた遺伝子が実際に良い遺伝子であった例」をよく求められるのですが、これはRNA-seqで見つけた6遺伝子のうち2個でknockoutすると良いことがあった例です。
あと、CHOの論文で5番目の著者がCHOさん。

PMID: 26854539
Targeted gene deletion using DNA-free RNA-guided Cas9 nuclease accelerates
adaptation of CHO cells to suspension culture
工業的なタンパク質生産に用いるCHO細胞は無血清培地に馴化させており、また、足場のない浮遊細胞を行います。著者らは、ゲノム編集技術を用いてCHO-K1細胞の無血清・浮遊培養への馴化促進を試みました。

まず開発の対象を得るため、CHO-K1細胞を無血清培地で浮遊培養して馴化させました。継代は10回行い、7回目ほどで馴化が完了したようです。馴化中の細胞及び馴化後の細胞についてRNA-seqを行い、DESeqを用いて発言変動遺伝子(DEG)を抽出しました。得られたDEGより、馴化に伴って発現の減少するDEGとして減少幅の大きい6遺伝子を選抜しました。
次に、それらの6遺伝子についてDNA-free RNA-guided Cas 9の系を用いてknockout したCHO-K1細胞を作成し、無血清・浮遊培養への馴化を行いました。結果、6遺伝子中2個の遺伝子については、それぞれのknochoutによって馴化が促進されました。また、アンプリコンシークエンスにより、馴化前後で細胞集団中のknochoutが成功している細胞の割合を調べたところ、馴化促進効果のあった2遺伝子については、knochoutの成功している細胞が集団中の割合を増加させていました。

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


元々CHO細胞の無血清浮遊化にかかる期間は1ヶ月程度ですので、ゲノム編集の手間を考えると実用性には疑問がありますが、RNA-seqで見つけた良さそうな遺伝子の答え合わせの例題としてクリアであると思います。通常のCRISPR/Cas9の系ではプラスミドやレンチウィルスの配列がホスト細胞ゲノムに導入されてしまうことがあり抗体医薬製造用の細胞開発では不適切であるため、CHO細胞のゲノム編集にはDNA-free な編集系の適用が必要であると著者らは述べており、ここが本研究の重心のようです。なので、DEGを得るためにどのサンプル間で比較したのかがわからなかった点、シークエンスした細胞のうち1つの実験区でのみ培養初期(培地に細胞を播種してすぐ)にサンプリングを行ってしまった点(なのでひとつだけPCAで外れてる)も致し方ないのかなと思いました。
[PR]

by koretoki | 2016-06-18 11:49
2014年 05月 23日

genome.fa とtranscripts.gtf でtranscripts.fa をつくる


genome.fa とtranscripts.gtf でtranscripts.fa をつくる

編集 | 削除
tophat-cufflinks を使ってtranscripts.gtf は手に入れたもののtranscripts.fa がないから遺伝子見るの大変、とか、何でか知らないけどtranscripts.gtf しか貰えなかったパターンのときにtranscripts.fa を作らないといけない状況になります。

#transcripts.gtf とgenome.fa からtranscripts.fa を作る。
#transcripts.gtf を加工1
more transcripts.gtf | awk -F "\t" '{print $1"\t"$3"\t"$4"\t"$5"\t"$9}' | grep exon | awk '{print $1"\t"$3"\t"$4"\t"$6"\t"$8"\t"$10}' | tr -d "\"" | tr -d ";" | awk '{print "transcript_id_"$5"\t"$1"\t"$2"\t"$3"\t"$6}' | awk '{print $1"\t"$2"\tsubstr($2,"$3","$4-$3+1")\t"$5/100"0"}' | awk '{print $1"\t"$2"\t"$3"\tn"substr($4,3,2)}' | sed -e '/n01/s/subs/@subs/g' | awk '{print $3}' | tr -d "\n" | tr "@" "\n" | grep substr > temp1.txt
#transcripts.gtf の加工2
more transcripts.gtf | awk -F "\t" '{print $1"\t"$3"\t"$4"\t"$5"\t"$9}' | grep exon | awk '{print $1"\t"$3"\t"$4"\t"$6"\t"$8"\t"$10}' | tr -d "\"" | tr -d ";" | awk '{print "transcript_id_"$5"\t"$1"\t"$2"\t"$3"\t"$6}' | awk '{print $1"\t"$2}' | uniq > temp2.txt

#genome.fa を1行にしておく
more genome.fa | sed -e '/>/s/$/@/g' | tr -d "\n" | tr ">" "\n" > line_genome.txt

#transcripts.gtf の情報を使って
paste temp2.txt temp1.txt | head -n 10
transcript_id_CUFF.1.1NODE_1000_length_89123_cov_46.018040substr($2,2,2497)
transcript_id_CUFF.2.1NODE_1000_length_89123_cov_46.018040substr($2,3291,1003)substr($2,4409,646)

#こういうコマンドをつくりたい
more line_genome.txt | grep NODE_1000_length_89123_cov_46.018040 | awk -F "@" '{print ">transcript_id_CUFF.1.1@"substr($2,2,2497)}' | tr "@" "\n" >> transcripts.fa
>transcript_id_CUFF.1.1
CGTAGTAGTATCGATCGTATCTACTATCGTACGTAGCTAGTCGATGCTAGTCGTAGCTAACTCGATCGATCGTAGTCTGTTAGCTAGCTGATCGTAGCTAGCTGATCGTAGCTAGCTGATCGTAGCTAGCTGATCGTAGCTGATCGTAGCTGATGCTGTAGCTGATCGTAGCTGATCGTAGCTAGTCTAGCTGATGCTAGCTGATCGTAGCTAGCTGATCGTAGCTAGCTGATGCTAGCTGATCGTAGCTGATCGTAGCTAGCTGATCgTAGCTAGCTTATCGTGTAGTCGTAGCTAGCTGATCGTAGCTAGATAAAAACGGGGATCGATCGTAGCTAAAAAGCTTTTTTTGAAAAATCGTAGCTGATGCTAGCTGATGCTGTAGTGTGTGATTTATATTACTAGCTAGCCCCATGCTAGTCGATCACACAC

#transcripts.gtf の情報を使ってコマンドを作る
paste temp2.txt temp1.txt | awk '{print "more line_genome.txt | grep "$2" | awk -F \"@\" %{print \">"$1"@\""$3"}% | tr \"@\" \"\\n\" >> transcripts.fa"}' | tr "%" "'" > make_transcripts.sh

#準備ができたら実行
sh make_transcripts.sh

#できる
ls

genome.fa
transcripts.gtf
temp1.txt
temp2.txt
line_genome.txt
make_transcripts.sh
transcripts.fa <- NEW!!!


終わってから気づいたけどこれでいいじゃん

http://cufflinks.cbcb.umd.edu/gff.html

Extracting transcript sequences
The gffread utility can be used to generate a FASTA file with the DNA sequences for all transcripts in a GFF file. For this operation a fasta file with the genomic sequences have to be provided as well. For example, one might want to extract the sequence of all transfrags assembled from a Cufflinks assembly session. This can be accomplished with a command line like this:

gffread -w transcripts.fa -g /path/to/genome.fa transcripts.gtf

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

[PR]

by koretoki | 2014-05-23 15:00
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年 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月 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