User Tools

Site Tools


variousperlscriptsjapanese

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revision Previous revision
Next revision
Previous revision
variousperlscriptsjapanese [2010/10/01 12:57]
ike
variousperlscriptsjapanese [2014/01/18 07:44] (current)
Line 2: Line 2:
  
 このチュートリアルはPerlプログラミングを用いて生物系統情報や複数のゲノム情報を扱う方法、系統樹の作成方法、統計的手法をつかった解析を行うためのサンプルスクリプトを示すものです。ここでは、'​g-language-1.8.10/​share/​genomes/'​に含まれるGenBankファイル(*.gbk)を利用します。使用する関数の一覧はここを参照して下さい。http://​kobesearch.cpan.org/​htdocs/​Bio-Glite/​Bio/​Glite.html このチュートリアルはPerlプログラミングを用いて生物系統情報や複数のゲノム情報を扱う方法、系統樹の作成方法、統計的手法をつかった解析を行うためのサンプルスクリプトを示すものです。ここでは、'​g-language-1.8.10/​share/​genomes/'​に含まれるGenBankファイル(*.gbk)を利用します。使用する関数の一覧はここを参照して下さい。http://​kobesearch.cpan.org/​htdocs/​Bio-Glite/​Bio/​Glite.html
 +====== Example 1 - 生物系統情報の取得 ​ ======
  
- +以下のPerlスクリプトは大腸菌のGenBankファイルを読み込み、生物系統情報を出力します。
-====== Example 1- Retrieving taxonomic information ​ ====== +
- +
-The following ​Perl script loads the GenBank-formatted file of E.coli genome, and prints taxonomic information.+
  
 <code perl> <code perl>
-use G; +  ​use G; 
-my $gb = load("​g-language-1.8.10/​share/​genomes/​ecoli.gbk"​);​ +  my $gb = load("​g-language-1.8.10/​share/​genomes/​ecoli.gbk"​);​ 
-foreach(sort keys %{$gb->​{TAXONOMY}}){ print "$_: ",​$gb->​{TAXONOMY}->​{$_},"​\n";​ }+  foreach(sort keys %{$gb->​{TAXONOMY}}) { 
 +    ​print "$_: ", $gb->​{TAXONOMY}->​{$_},​ "​\n";​ 
 +  ​}
 </​code>​ </​code>​
  
-Executing this script prints the following output.+このスクリプトを実行すると以下のような出力結果が得られます。
  
       ​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​       ​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​__/​
Line 64: Line 64:
   species: coli   species: coli
  
-Thus, $gb->​{TAXONOMY}->​{0} ​contains '​all'​ taxonomic information,​ and taxonomic ranks 1234, and correspond with domain, phylum, class, order, and family, respectively.+このように、$gb->​{TAXONOMY}->​{0}には系統の項目すべてが含まれ、個別の項目にアクセスするためにはランクを0から12345に変更することで、ドメイン、門、網、目、科の名称を取得することができます。
  
-This taxonomic information is simply checked by typing:+この系統情報はシェルこのように入力することで簡単に確認することができます。
  
 +<code sh>
   $ grep -A 2 ORGANISM g-language-1.8.10/​share/​genomes/​ecoli.gbk   $ grep -A 2 ORGANISM g-language-1.8.10/​share/​genomes/​ecoli.gbk
     ORGANISM ​ Escherichia coli str. K-12 substr. MG1655     ORGANISM ​ Escherichia coli str. K-12 substr. MG1655
               Bacteria; Proteobacteria;​ Gammaproteobacteria;​ Enterobacteriales;​               Bacteria; Proteobacteria;​ Gammaproteobacteria;​ Enterobacteriales;​
               Enterobacteriaceae;​ Escherichia.               Enterobacteriaceae;​ Escherichia.
 +</​code>​
 +====== Example 2 - 複数のGenbankファイルを用いた解析 ​ ======
  
- +以下のスクリプトは任意のディレクトリ(ここでは'​g-language-1.8.10/​share/​genomes/'​)に含まれるすべてのGenBankファイル(*.gbk)を読み込み、生物種の定義(DEFINITION)を出力します。
-====== Example 2- Multiple genome analysis ​ ====== +
- +
-The following script loads all GenBank files (*.gbk) in the directory ​($dir = '​g-language-1.8.10/​share/​genomes/'​), and prints ​DEFINITION ​of these genomes.+
  
 <code perl> <code perl>
-use G; +  ​use G; 
-my $dir = '​g-language-1.8.10/​share/​genomes/';​ +  my $dir = '​g-language-1.8.10/​share/​genomes/';​ 
-opendir(DIR,​ $dir) || die "​directory open error:​$!";​ +  opendir(DIR,​ $dir) || die "​directory open error:​$!";​ 
-foreach my $file (sort readdir(DIR)){ +  foreach my $file (sort readdir(DIR)) { 
-    next if($file !~ /​\.gbk$/​);​ +    next if ($file !~ /​\.gbk$/​);​ 
-    my $gb = load("​$dir/​$file","​no msg"​);​ +    my $gb = load("​$dir/​$file",​ "no msg"​);​ 
-    print "$gb->​{DEFINITION}\n";​ +    print $gb->​{DEFINITION}, "\n"; 
-+  
-closedir(DIR);​+  closedir(DIR);​
 </​code>​ </​code>​
  
-Executing this script prints the following output.+このスクリプトを実行すると以下のような出力結果が得られます。
  
   Borrelia burgdorferi B31, complete genome.   Borrelia burgdorferi B31, complete genome.
Line 101: Line 101:
   Pyrococcus furiosus DSM 3638, complete genome.   Pyrococcus furiosus DSM 3638, complete genome.
  
-The DEFINITION is quickly checked by typing:+生物種の定義はGenBankファイルに記述されているので、以下のようなシェルコマンドでシンプルに確認することもできます。
  
   $ grep DEFINITION g-language-1.8.10/​share/​genomes/​*.gbk   $ grep DEFINITION g-language-1.8.10/​share/​genomes/​*.gbk
  
  
 +ジヌクレオチドの相対的な存在量(ゲノミックシグニチャー)は存在することが予測されるジヌクレオチド量との比として表されます。以下のように求めることができます。
  
-Dinucleotide relative abundance (genomic signature) is a ratio of observed to expected dinucleotide frequencies. It can be calculated by: 
 <code perl> <code perl>
-    ​&​signature($gb);​+  ​&​signature($gb);​
 </​code>​ </​code>​
  
-Executing the script should generate a vector of start and end positions, ​G+C content, and the relative abundance values of 16 dinucleotides for a whole genome region, as shown in http://​www.pnas.org/​content/​96/​16/​9184/​F1.large.jpg.+このスクリプトを実行すると、すべてのゲノムに対して開始位置、終結位置、G+C組成、16種類のジヌクレオチドの相対的存在量を行列として生成します(http://​www.pnas.org/​content/​96/​16/​9184/​F1.large.jpg)。
  
 +相対的同義語コドン使用頻度(Relative synonymous codon usage; RSCU)は予測されるコドン使用頻度との比として表されるため、コドンの使用頻度に偏りが生じないときは1となります。RSCUは以下のように求めることができます。
  
 +スクリプトを実行すると、以下のように各遺伝子のすべてのコドン(開始コドンと終止コドンは除く)に対してRSCU値の行列を出力します。
  
-Relative synonymous codon usage (RSCU) is a ratio of observed to expected codon frequencies,​ and is equal to 1 in the absence of any codon usage bias. It can be calculated by: 
 <code perl> <code perl>
-    ​&​codon_compiler($gb,​ -output=>'​stdout',​ -data=>'​RSCU'​);​+  ​&​codon_compiler($gb,​ -output=>'​stdout',​ -data=>'​RSCU'​);​
 </​code>​ </​code>​
  
-Executing the script should generate a vector of the RSCU values of all codons ​(omitting start and stop codons) for a collection of all protein-coding genes in each genome. +ライム病ボレリア(Borrelia burgdorferi)のゲノムに対するRSCU値は論文の表1に示しています(http://​nar.oxfordjournals.org/​cgi/​content/​full/​27/​7/​1642)。 
-The RSCU values for Borrelia burgdorferi ​genome are shown in Table at http://​nar.oxfordjournals.org/​cgi/​content/​full/​27/​7/​1642+====== Example 3 - rRNA配列から系統樹の作成 ======
  
- +$gb->​rRNA()はrRNAのすべてのfeature ​IDを配列として返し、'​16S'​'​23S'​'​5S'​'​SSU'​'​LSU'​のように、特定のrRNAを指定することもできます。 
-====== Example 3- Retrieving rRNA sequences and phylogenetic tree reconstruction ​ ====== +例えば、$gb->​rRNA('​16S'​)とすることで塩基配列順にソートされた16S rRNAのfeature ​IDを配列として返します。 
- +以下のスクリプトを実行すると、各生物種の中からもっとも長い塩基配列の16S rRNAをfasta形式のファイルとして出力します。
-$gb->​rRNA() ​returns the array of all feature ​IDs of rRNAs, and can specify rRNA species (e.g. '​16S'​'​23S'​'​5S'​'​SSU'​'​LSU'​). For exmaple, ​$gb->​rRNA('​16S'​) ​returns a list of feature ​IDs of 16S rRNAs, sorted by the nucleotide sequence lengths. +
- +
-The following script outputs a fasta-formatted file of nucleotide sequences of the longest ​16S rRNA genes for multiple bacterial strains.+
  
 <code perl> <code perl>
-use G;+  ​use G;
  
-my $rRNA = '​16S';​+  ​my $rRNA = '​16S';​
  
-open(OUT,">​gene.fasta"​);​+  ​open(OUT,">​gene.fasta"​);​
  
-my $dir = '​g-language-1.8.10/​share/​genomes/';​ +  ​my $dir = '​g-language-1.8.10/​share/​genomes/';​ 
-opendir(DIR,​ $dir) || die "​directory open error:​$!";​ +  opendir(DIR,​ $dir) || die "​directory open error:​$!";​ 
-foreach my $file (sort readdir(DIR)){+  foreach my $file (sort readdir(DIR)){
     next if($file !~ /\.gbk$/);     next if($file !~ /\.gbk$/);
     my $gb = load("​$dir/​$file","​no msg");     my $gb = load("​$dir/​$file","​no msg");
     next unless($gb->​rRNA($rRNA));​     next unless($gb->​rRNA($rRNA));​
     print OUT ">​$gb->​{ORGANISM}\n",​ $gb->​get_geneseq($gb->​rRNA($rRNA)),"​\n";​     print OUT ">​$gb->​{ORGANISM}\n",​ $gb->​get_geneseq($gb->​rRNA($rRNA)),"​\n";​
-+  ​
-closedir(DIR);​+  closedir(DIR);​
  
-close(OUT);+  ​close(OUT);
  
 </​code>​ </​code>​
  
-To check the number of sequences in the fasta-formatted output file (gene.fasta),​ type:+fastaファイルにいくつの配列が含まれるのか確認するためには、シェルで以下のように入力してください。
  
   $ grep '>'​ gene.fasta | wc -l   $ grep '>'​ gene.fasta | wc -l
   7   7
  
-Phylogenetic tree can be reconstructed by adding the following two lines: the first line invokes multiple alignment, and the second line calculates Neighbor joining ​(NJ) tree.+系統樹を作成するためには2行だけ以下に示すスクリプトを書き加えてください。1行目でマルチプルアラインメントを実行し、2行目で近隣結合法(neighbor-joining method; ​NJ)を計算します。
  
 <code perl> <code perl>
Line 166: Line 164:
 </​code>​ </​code>​
  
-Executing this script generates the following three output files: the gene.aln ​file (multiple alignment)gene.dnd ​file (guide tree), and gene.ph ​file (NJ tree)+このスクリプトを実行することで、3つのファイルを出力します。gene.aln(マルチアライメントファイル)gene.dnd(デンドログラムファイル)gene.ph(推定された系統樹ファイル)です。以下の系統樹はgene.phTreeView(http://​taxonomy.zoology.gla.ac.uk/​rod/​treeview.html)で描画したものです。
-The following phylogram was drawn using the gene.ph ​file as an input file of TreeView (http://​taxonomy.zoology.gla.ac.uk/​rod/​treeview.html).+
  
 {{hs_nj_rrs.png?​300|}} {{hs_nj_rrs.png?​300|}}
 +====== Example 4 - ゲノム配列の情報と相関解析 ======
  
- +以下のスクリプトはさまざまなゲノム情報を解析します。rRNAtRNA遺伝子、CDSの数、ゲノムサイズ、塩基組成[100x(G+C)/​(A+C+G+T)]GC skew [(G-C)/​(G+C)] ​の度合いを示すgcsi ​(GC skew index; ​http://​www.ncbi.nlm.nih.gov/​pubmed/​19461976)といった解析です。
-====== Example 4- Calculating genome statistics and their correlation analyses ​ ====== +
- +
-The following script calculates various genome statistics such as rRNA gene number, ​tRNA gene number, protein-coding sequence (cds) number, genome size, GC content ​[100x(G+C)/​(A+C+G+T)], and GC skew index (gcsi), which is a measure of the degree of GC skew [(G-C)/​(G+C)] (http://​www.ncbi.nlm.nih.gov/​pubmed/​19461976).+
  
 <code perl> <code perl>
-use G;+  ​use G;
  
-my @key = qw(rRNA tRNA cds size GC gcsi); +  ​my @key = qw(rRNA tRNA cds size GC gcsi); 
-my %stat;+  my %stat;
  
-open(OUT,">​stat.txt"​);​ +  ​open(OUT,">​stat.txt"​);​ 
-print OUT join("​\t",​ @key), "​\n";​+  print OUT join("​\t",​ @key), "​\n";​
  
-my $dir = '​g-language-1.8.10/​share/​genomes/';​ +  ​my $dir = '​g-language-1.8.10/​share/​genomes/';​ 
-opendir(DIR,​ $dir) || die "​directory open error:​$!";​ +  opendir(DIR,​ $dir) || die "​directory open error:​$!";​ 
-foreach my $file (sort readdir(DIR)){+  foreach my $file (sort readdir(DIR)){
     next if($file !~ /\.gbk$/);     next if($file !~ /\.gbk$/);
     my $gb = load("​$dir/​$file","​no msg");     my $gb = load("​$dir/​$file","​no msg");
Line 205: Line 200:
     foreach (@key){ push(@tmp, $stat{$_}); }     foreach (@key){ push(@tmp, $stat{$_}); }
     print OUT join("​\t",​ @tmp), "​\n";​     print OUT join("​\t",​ @tmp), "​\n";​
-+  ​
-closedir(DIR);​+  closedir(DIR);​
  
-close(OUT);+  ​close(OUT);
 </​code>​ </​code>​
  
-Executing this script generates the following tab-delimited output file ('​stat.txt'​).+このスクリプトを実行すると、以下のようなタブ区切りの結果がファイル('​stat.txt'​)に保存されます。
  
   rRNA    tRNA    cds     ​size ​   GC      gcsi   rRNA    tRNA    cds     ​size ​   GC      gcsi
Line 223: Line 218:
   4       ​46 ​     2125    1908256 40.8    0.0203799266306591   4       ​46 ​     2125    1908256 40.8    0.0203799266306591
  
-Scatter plots and correlation coefficients between all pairs of these values can be computed by adding the following ​commands ​(http://​www.bioconductor.org/​mogr/​chapter-code/​TwoColorPre.R).+得られた結果の散布図や相関係数は以下のようなRコマンド(http://​www.bioconductor.org/​mogr/​chapter-code/​TwoColorPre.R)を書き加えることで計算できます。
  
 <code perl> <code perl>
-my $rcmd = new Rcmd; +  ​my $rcmd = new Rcmd; 
-my @result = $rcmd->​exec(+  my @result = $rcmd->​exec(
                          qq!                          qq!
                          dat = read.delim("​stat.txt"​)                          dat = read.delim("​stat.txt"​)
Line 239: Line 234:
 </​code>​ </​code>​
  
-Executing this script generates the following '​Rplot.ps'​ file.+実行するとこのようなグラフが生成されます。
  
 {{genomestatrplot.png?​600|}} {{genomestatrplot.png?​600|}}
  
-There are a strong positive correlation between rRNA gene number and tRNA gene number (see rRNA-tRNA plot, r=0.9), as reported previouly ​(http://​nar.oxfordjournals.org/​cgi/​content/​full/​33/​4/​1141/​FIG4), and between the number of protein-coding sequences and genome size (see cds-size plot, r=1), as reported previouly ​(http://​www.pnas.org/​content/​101/​9/​3160/​F2.expansion.html).+先行研究で指摘されるように(http://​nar.oxfordjournals.org/​cgi/​content/​full/​33/​4/​1141/​FIG4)、rRNA遺伝子数とtRNA遺伝子数に強い正の相関があることがわかります 
 +(rRNA-tRNAのプロットが相関係数r=0.9を示している)。同様に、別の先行研究で指摘されていますように(http://​www.pnas.org/​content/​101/​9/​3160/​F2.expansion.html)、CDS数とゲノムサイズにも強い正の相関が見られます(cds-sizeのプロットが相関係数r=1を示している)。 
 +====== Example 5 - 遺伝子情報の統計解析 ======
  
- +以下のスクリプトは、マイコプラズマ菌(Mycoplasma genitalium)の遺伝子についてさまざまな統計的な解析をおこなったものです。この解析には開始位置、コドンの3番目の塩基におけるGC組成、コドンの3番目の塩基におけるGC skew、遺伝子発現量を推定するcodon adaptation index (caihttp://​www.ncbi.nlm.nih.gov/​pubmed/​3547335)、同義コドンの均一性を示す重み付けされた相対エントロピーの和(http://​www.ncbi.nlm.nih.gov/​pubmed/​15194186)が含まれます。
-====== Example 5- Calculating gene statistics ​ ====== +
- +
-The following script prints various statistics for protein-coding sequences ​(cds) in Mycoplasma genitalium ​genome such as start position (start)GC content at 3rd codon positions (gcc3), ​GC skew at 3rd codon positions (gcs3), the codon adaptation index (cai) (http://​www.ncbi.nlm.nih.gov/​pubmed/​3547335), and the weighted sum of relative entropy (Ew) as a measure of the degree of synonymous codon evenness ​(http://​www.ncbi.nlm.nih.gov/​pubmed/​15194186).+
  
 <code perl> <code perl>
-use G;+  ​use G;
  
-my $gb = load("​g-language-1.8.10/​share/​genomes/​mgen.gbk"​);​+  ​my $gb = load("​g-language-1.8.10/​share/​genomes/​mgen.gbk"​);​
  
-foreach my $cds ($gb->​cds()){+  ​foreach my $cds ($gb->​cds()){
     &​bui($gb,​ -output=>'​n',​ -id=>​$cds,​ -position=>​3);​ # base usage indices (gcc3, gcs3) at 3rd codon posotions     &​bui($gb,​ -output=>'​n',​ -id=>​$cds,​ -position=>​3);​ # base usage indices (gcc3, gcs3) at 3rd codon posotions
-+  ​
-&​cai($gb,​ -output=>'​n',​ -w_output=>'​n'​);​ +  &​cai($gb,​ -output=>'​n',​ -w_output=>'​n'​);​ 
-&​Ew($gb,​ -output=>'​n'​);​+  &​Ew($gb,​ -output=>'​n'​);​
  
-my @key = qw(start gcc3 gcs3 cai Ew); +  ​my @key = qw(start gcc3 gcs3 cai Ew); 
-open(OUT,">​stat.txt"​);​ +  open(OUT,">​stat.txt"​);​ 
-print OUT join("​\t",​ @key), "​\n";​ +  print OUT join("​\t",​ @key), "​\n";​ 
-foreach my $cds ($gb->​cds()){+  foreach my $cds ($gb->​cds()){
     my @tmp;     my @tmp;
     foreach (@key){ push(@tmp, $gb->​{$cds}->​{$_});​ }     foreach (@key){ push(@tmp, $gb->​{$cds}->​{$_});​ }
     print OUT join("​\t",​ @tmp), "​\n";​     print OUT join("​\t",​ @tmp), "​\n";​
-+  ​
-close(OUT);+  close(OUT);
 </​code>​ </​code>​
  
-Executing this script generates the following tab-delimited output file ('​stat.txt'​).+このスクリプトを実行すると以下のようなタブ区切りの結果が'​stat.txt'​として保存されます。
  
   $ head stat.txt   $ head stat.txt
Line 286: Line 280:
   11251   ​0.1533 ​ -0.1500 0.7432 ​ 0.6235   11251   ​0.1533 ​ -0.1500 0.7432 ​ 0.6235
  
-Scatter plots and correlation coefficients between all pairs of these values can be computed by adding the R commands described above (Example 4). +得られた結果の散布図と相関係数は、Example 4で示したようにRコマンドを書き加えることで描画されます。このスクリプトを実行すると'​Rplot.ps'​というファイルが出力されます。
-Executing this script generates the following ​'​Rplot.ps' ​file.+
  
 {{cdsstatrplot_mgen.png?​600|}} {{cdsstatrplot_mgen.png?​600|}}
  
-M. genitalium has a very distinctive variation in GC content within the genome (see start-gcc3 ​plot), and GC content is correlated with codon usage bias (see gcc3-cai ​and gcc3-Ew ​plots), as reported previously ​(Figure 2 in http://​www.horizonpress.com/​cimb/​v/​v3/​v3n403.pdf)+図のstart-gcc3プロットより、マイコプラズマ菌のGC組成は特徴的な傾向を示していることが分かります。そして、gcc3-caigcc3-Ewプロっトより、GC組成とコドン使用頻度の偏りに相関があることが分かります。これらは先行研究の知見とも一致します(http://​www.horizonpress.com/​cimb/​v/​v3/​v3n403.pdf・図2)。ライム病ボレリア(Borrelia burgdorferi"​g-language-1.8.10/​share/​genomes/​bbur.gbk"​)ではGC skewとコドン使用頻度の偏りに相関が見られます。 
-In Borrelia burgdorferi ​("​g-language-1.8.10/​share/​genomes/​bbur.gbk"​)GC skew (gcs3) is correlated with codon usage bias.+
variousperlscriptsjapanese.1285937836.txt.gz · Last modified: 2014/01/18 07:44 (external edit)