User Tools

Site Tools


variousperlscriptsjapanese

This is an old revision of the document!


はじめに

このチュートリアルはPerlプログラミングを用いて生物系統情報や複数のゲノム情報を扱う方法、系統樹の作成方法、統計的手法をつかった解析を行うためのサンプルスクリプトを示すものです。ここでは、'g-language-1.8.10/share/genomes/'に含まれるGenBankファイル(*.gbk)を利用します。使用する関数の一覧はここを参照して下さい。http://kobesearch.cpan.org/htdocs/Bio-Glite/Bio/Glite.html

Example 1 - 生物系統情報の取得

以下のPerlスクリプトは大腸菌のGenBankファイルを読み込み、生物系統情報を出力します。

  use G;
  my $gb = load("g-language-1.8.10/share/genomes/ecoli.gbk");
  foreach(sort keys %{$gb->{TAXONOMY}}) {
    print "$_: ", $gb->{TAXONOMY}->{$_}, "\n";
  }

このスクリプトを実行すると以下のような出力結果が得られます。

     __/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/
              
                 G-language  Genome Analysis Environment v.1.8.10


                           http://www.g-language.org/

            Please cite: 
               Arakawa K. et al. (2003) Bioinformatics.
               Arakawa K. et al. (2006) Journal of Pestice Science.
               Arakawa K. et al. (2008) Genes, Genomes and Genomics.

	      License: GNU General Public License
      Copyright (C) 2001-2009 G-language Project
      Institute for Advanced Biosciences, Keio University, JAPAN 

     __/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/
	   



Accession Number: NC_000913

  Length of Sequence :   4639675
           A Content :   1142228 (24.62%)
           T Content :   1140970 (24.59%)
           G Content :   1176923 (25.37%)
           C Content :   1179554 (25.42%)
              Others :         0 (0.00%)
          AT Content :    49.21%
          GC Content :    50.79%

0: Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacteriales;Enterobacteriaceae;Escherichia.
1: Bacteria
2: Proteobacteria
3: Gammaproteobacteria
4: Enterobacteriales
5: Enterobacteriaceae
6: Escherichia.
all: Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacteriales;Enterobacteriaceae;Escherichia.
class: Gammaproteobacteria
domain: Bacteria
family: Enterobacteriaceae
genus: Escherichia
order: Enterobacteriales
phylum: Proteobacteria
species: coli

このように、$gb→{TAXONOMY}→{0}には系統の項目すべてが含まれ、個別の項目にアクセスするためにはランクを0から1、2、3、4、5に変更することで、ドメイン、門、網、目、科の名称を取得することができます。

この系統情報はシェルこのように入力することで簡単に確認することができます。

  $ grep -A 2 ORGANISM g-language-1.8.10/share/genomes/ecoli.gbk
    ORGANISM  Escherichia coli str. K-12 substr. MG1655
              Bacteria; Proteobacteria; Gammaproteobacteria; Enterobacteriales;
              Enterobacteriaceae; Escherichia.

Example 2 - 複数のGenbankファイルを用いた解析

以下のスクリプトは任意のディレクトリ(ここでは'g-language-1.8.10/share/genomes/')に含まれるすべてのGenBankファイル(*.gbk)を読み込み、生物種の定義(DEFINITION)を出力します。

  use G;
  my $dir = 'g-language-1.8.10/share/genomes/';
  opendir(DIR, $dir) || die "directory open error:$!";
  foreach my $file (sort readdir(DIR)) {
    next if ($file !~ /\.gbk$/);
    my $gb = load("$dir/$file", "no msg");
    print $gb->{DEFINITION}, "\n";
  }
  closedir(DIR);

このスクリプトを実行すると以下のような出力結果が得られます。

Borrelia burgdorferi B31, complete genome.
Bacillus subtilis subsp. subtilis str. 168, complete genome.
Buchnera aphidicola str. APS (Acyrthosiphon pisum), complete genome.
Synechococcus sp. WH 8102, complete genome.
Escherichia coli str. K-12 substr. MG1655, complete genome.
Mycoplasma genitalium G37, complete genome.
Plasmid F, complete sequence.
Pyrococcus furiosus DSM 3638, complete genome.

生物種の定義はGenBankファイルに記述されているので、以下のようなシェルコマンドでシンプルに確認することもできます。

$ grep DEFINITION g-language-1.8.10/share/genomes/*.gbk

ジヌクレオチドの相対的な存在量(ゲノミックシグニチャー)は存在することが予測されるジヌクレオチド量との比として表されます。以下のように求めることができます。

  &signature($gb);

このスクリプトを実行すると、すべてのゲノムに対して開始位置、終結位置、G+C組成、16種類のジヌクレオチドの相対的存在量を行列として生成します(http://www.pnas.org/content/96/16/9184/F1.large.jpg)。

相対的同義語コドン使用頻度(Relative synonymous codon usage; RSCU)は予測されるコドン使用頻度との比として表されるため、コドンの使用頻度に偏りが生じないときは1となります。RSCUは以下のように求めることができます。

スクリプトを実行すると、以下のように各遺伝子のすべてのコドン(開始コドンと終止コドンは除く)に対してRSCU値の行列を出力します。

  &codon_compiler($gb, -output=>'stdout', -data=>'RSCU');

ライム病ボレリア(Borrelia burgdorferi)のゲノムに対するRSCU値は論文の表1に示しています(http://nar.oxfordjournals.org/cgi/content/full/27/7/1642)。

Example 3 - rRNA配列から系統樹の作成

$gb→rRNA()はrRNAのすべてのfeature IDを配列として返し、'16S'、'23S'、'5S'、'SSU'、'LSU'のように、特定のrRNAを指定することもできます。 例えば、$gb→rRNA('16S')とすることで塩基配列順にソートされた16S rRNAのfeature IDを配列として返します。 以下のスクリプトを実行すると、各生物種の中からもっとも長い塩基配列の16S rRNAをfasta形式のファイルとして出力します。

  use G;
 
  my $rRNA = '16S';
 
  open(OUT,">gene.fasta");
 
  my $dir = 'g-language-1.8.10/share/genomes/';
  opendir(DIR, $dir) || die "directory open error:$!";
  foreach my $file (sort readdir(DIR)){
    next if($file !~ /\.gbk$/);
    my $gb = load("$dir/$file","no msg");
    next unless($gb->rRNA($rRNA));
    print OUT ">$gb->{ORGANISM}\n", $gb->get_geneseq($gb->rRNA($rRNA)),"\n";
  }
  closedir(DIR);
 
  close(OUT);

fastaファイルにいくつの配列が含まれるのか確認するためには、シェルで以下のように入力してください。

$ grep '>' gene.fasta | wc -l
7

系統樹を作成するためには2行だけ以下に示すスクリプトを書き加えてください。1行目でマルチプルアラインメントを実行し、2行目で近隣結合法(neighbor-joining method; NJ法)を計算します。

system(qq!
       clustalw -align -infile=gene.fasta -outfile=gene.aln -outorder=input
       clustalw -tree -infile=gene.aln
       !);

このスクリプトを実行することで、3つのファイルを出力します。gene.aln(マルチアライメントファイル)とgene.dnd(デンドログラムファイル)とgene.ph(推定された系統樹ファイル)です。以下の系統樹はgene.phをTreeView(http://taxonomy.zoology.gla.ac.uk/rod/treeview.html)で描画したものです。

Example 4 - ゲノム配列の情報と相関解析

以下のスクリプトはさまざまなゲノム情報を解析します。rRNAやtRNA遺伝子、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)といった解析です。

  use G;
 
  my @key = qw(rRNA tRNA cds size GC gcsi);
  my %stat;
 
  open(OUT,">stat.txt");
  print OUT join("\t", @key), "\n";
 
  my $dir = 'g-language-1.8.10/share/genomes/';
  opendir(DIR, $dir) || die "directory open error:$!";
  foreach my $file (sort readdir(DIR)){
    next if($file !~ /\.gbk$/);
    my $gb = load("$dir/$file","no msg");
 
    $stat{rRNA} = scalar $gb->rRNA(); #	rRNA gene number
    $stat{tRNA} = scalar $gb->tRNA(); #	tRNA gene number
    $stat{cds} = scalar $gb->cds(); # protein-coding sequence (cds) number
    $stat{size} = length($gb->{SEQ}); #	genome size
 
    my $acgt = $gb->{SEQ} =~ tr/acgt/acgt/; # A+C+G+T
    my $gc = $gb->{SEQ} =~ tr/gcGC/gcGC/; # G+C
    $stat{GC} = sprintf "%.1f", 100 * $gc / $acgt; # GC	content [100x(G+C)/(A+C+G+T)]
 
    $stat{gcsi} = (&gcsi($gb))[0]; # GC	skew index (gcsi)
 
    my @tmp;
    foreach (@key){ push(@tmp, $stat{$_}); }
    print OUT join("\t", @tmp), "\n";
  }
  closedir(DIR);
 
  close(OUT);

このスクリプトを実行すると、以下のようなタブ区切りの結果がファイル('stat.txt')に保存されます。

rRNA    tRNA    cds     size    GC      gcsi
5       18      851     910724  28.6    0.44353725960019
30      86      4176    4215606 43.5    0.214855185905019
3       32      564     640681  26.3    0.126305531560014
6       44      2519    2434428 59.4    0.0837418081678612
22      89      4321    4639675 50.8    0.0966615833014818
3       36      476     580076  31.7    0.127587425631069
0       0       108     99159   48.2    0.20973363827319
4       46      2125    1908256 40.8    0.0203799266306591

得られた結果の散布図や相関係数は以下のようなRコマンド(http://www.bioconductor.org/mogr/chapter-code/TwoColorPre.R)を書き加えることで計算できます。

  my $rcmd = new Rcmd;
  my @result = $rcmd->exec(
                         qq!
                         dat = read.delim("stat.txt")
                         postscript("Rplot.ps", pointsize=10, width=10, height=50)
                         upPan <- function(...) points(..., col="darkblue")
                         lowPan <- function(x, y, ...) text(mean(par("usr")[1:2]), mean(par("usr")[3:4]),signif(cor(x, y),2),cex=2)
                         pairs(dat, pch=20, lower.panel=lowPan, upper.panel=upPan)
                         dev.off()
                         !
                         );

実行するとこのようなグラフが生成されます。

先行研究で指摘されるように(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- 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).

use G;
 
my $gb = load("g-language-1.8.10/share/genomes/mgen.gbk");
 
foreach my $cds ($gb->cds()){
    &bui($gb, -output=>'n', -id=>$cds, -position=>3); # base usage indices (gcc3, gcs3) at 3rd codon posotions
}
&cai($gb, -output=>'n', -w_output=>'n');
&Ew($gb, -output=>'n');
 
my @key = qw(start gcc3 gcs3 cai Ew);
open(OUT,">stat.txt");
print OUT join("\t", @key), "\n";
foreach my $cds ($gb->cds()){
    my @tmp;
    foreach (@key){ push(@tmp, $gb->{$cds}->{$_}); }
    print OUT join("\t", @tmp), "\n";
}
close(OUT);

Executing this script generates the following tab-delimited output file ('stat.txt').

$ head stat.txt
start   gcc3    gcs3    cai     Ew
686     0.1583  -0.1333 0.7351  0.6295
1828    0.1780  -0.1636 0.7443  0.6675
2845    0.2296  -0.1409 0.7257  0.7239
4812    0.1868  -0.0769 0.7164  0.7150
7294    0.1827  +0.0526 0.7566  0.6604
8551    0.1818  +0.0526 0.7391  0.6084
9156    0.1304  -0.3939 0.7644  0.5625
9923    0.1519  +0.0746 0.7635  0.6090
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). Executing this script generates the following 'Rplot.ps' file.

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). In Borrelia burgdorferi ("g-language-1.8.10/share/genomes/bbur.gbk"), GC skew (gcs3) is correlated with codon usage bias.

variousperlscriptsjapanese.1285938443.txt.gz · Last modified: 2014/01/18 07:44 (external edit)