1 RNA Sequence Data

2014 年 5 月 23 日第 7 回
1
統計科学同演習第 7 回演習問題
(鳥インフルエンザの RNA sequence）
学籍番号:
氏名：
1
RNA Sequence Data
鳥インフルエンザ（Avian influenza, Avian flu, bird flu）とは，A 型インフルエンザウイルスが鳥類に感染
して起きる鳥類の感染症である（ウィキペディアより）．Avian influenza, sometimes avian flu, and commonly
bird flu, refers to “influenza caused by viruses adapted to birds.” ...... While its most highly pathogenic
strain (H5N1) had been spreading throughout Asia since 2003, Avian Influenza reached Europe in 2005,
and the Middle East as well as Africa the following year. （Wikipedia より）
演習のページに置いてある H5N1.fa は，鳥インフルエンザの RNA sequence （変異）のファイルであり，
次のようにして National Center for Biotechnology Information (NCBI) のサイトから入手したものであ
る．まず，
http://www.ncbi.nlm.nih.gov/genomes/FLU/
のサイトを開き，青字で書かれた Database をクリックする．開いたサイトで Nucleotide sequence を選び,
Type: A, Host: Avian, Country/Region: any, Segment: HA, Subtype H:5, Subtype N:1, Sequence length
Max:1742, Min:1742, Collection Date: 2009 年 1 月 1 日から 2014 年 5 月 20 日を指定し，Show Results
をクリックして，該当する sequence の一覧を得る．そのうえで，Nucleotide (FASTA) を指定し，ダウン
ロードした．
演習では，既に用意されている H5N1.fa を用いること（NCBI のデータは日々更新されているため）．
問題 1. ファイル H5N1.fa をテキストエディタでながめ，どんなファイル形式となっているか説明しなさい．
2
R への取り込み
パッケージ seqinr をインストールし，library(seqinr) で使用できるようにしなさい．このライブラ
リーにはゲノム配列を解析するのに必要なさまざまな関数が定義されており，関数 read.fasta は，FASTA
形式のファイルを R に取り込む関数であるが，これは各ゲノム配列をリストの枝として取り込むので，少々
不便である．そこで演習では，ファイル e7.R に定義された read.fasta2 の方を用いることにする．他にも
必要な諸関数定義がこのファイルにあるので，このファイルの内容を source で取り込みなさい．その上で,
>H5N1=read.fasta2("H5N1.fa")
によってファイル H5N1.fa をオブジェクト H5N1 として取り込みなさい.
2014 年 5 月 23 日第 7 回
2
問題 2. オブジェクト H5N1 は行列であるが，行数，列数と各行は何を表しているか答えなさい．
関数 Dist(H5N1) で変異間の一つの距離を計算できる．結果をオブジェクト dd に付値しなさい．
（この計
算には多少時間がかかるかもしれない．
）
問題 3. どのような距離を求めたのか，関数 Dist の定義や dd の値を見て説明しなさい．なお、dd の値を
部分的に見るには，as.matrix(dd) のように行列に直し，添字演算で取り出す必要がある．
さらに関数 hclust にオブジェクト dd を与えれば，階層型クラスタリングした結果が得られる．これを
オブジェクト hh に付値しなさい．すると，plot(hh) ( plclust(hh) でもよい）で階層型クラスタリング
の結果を図示できる．
この図は葉の数が多すぎてやや眺めにくいので，不要なデータを除くことを考える．
問題 4. オブジェクト hh を用いると，距離 0，つまり完全一致を条件にした変異のグループ分けができる：
> ct = cutree(hh, h=0.1)
> lab = dimnames(H5N1)[[1]]
# lab = rownames(H5N1) でもよい
> split(lab, ct)
なぜ，これで目的のグループ分けができるのか説明するとともに，結果的にどのような特徴をもったグルー
プ分けになっているか記しなさい．
（ヒント：plot(hh) の後に rect.hclust(hh, h=0.1) を実行してみよ）
2014 年 5 月 23 日第 7 回
3
e7.R で定義されている関数 unique.fasta を使って
> H5N1.u = unique.fasta(H5N1)
としなさい．これによって「塩基配列と名前が完全に一致するような 2 つの変異」は同一視される．
（行数
は 96 となるはずである．
）
問題 5. H5N1.u に対して階層型クラスタリングを実行した結果を hh.u とおき，その結果を plot(hh.u)
で表示しなさい．ただし，各葉に適切な大きさのラベルをつけて図示するには
> plot(hh.u, labels=rownames(H5N1.u), cex=0.5)
などとするとよい．この結果を印刷し，このクラスタ木の概要を説明しなさい．
今日の演習の感想:
この用紙と印刷出力を提出せよ．

Download Report