1. 解析の準備
C: においてnucleosomeというフォルダを作成し、そこにダウンロードしたすべてのファイルを入れてください。
Rを立ち上げ、getwd( )と打ち込みworking directoryを確認してください。
working directoryをC:/nucleosomeに変更してください（ディレクトリの変更はsetwd("C:/nucleosome")とすることもできます）。
2. ヌクレオソームを形成していた網羅的位置情報
生物統計解析基礎の9. 行列とデータフレームにおいて、Mixia osmundaeのコンティグ0009（BABT02000007）を使用しましたので、このコンティグにマップされたヌクレオソームを抽出します。
nuc9.txtというファイルをつけていますので、こちらを使用してください。
これをダウンロードして、nucleosomeのフォルダに入れてください。
Rを起動して、ディレクトリの変更を行い、nucleosomeにしてください。
dir( )として、先ほどのファイルが入っていることを確認してください。
nuc9 <- read.table("nuc9.txt")
ファイルの様子を見て見ましょう。
nuc9[1:10, ] #最初の10行を見てみます。
このファイルは1列目にヒストン８量体に結合していたDNA断片のコンティグ配列上での最初の位置、2列目に最後の位置を示しています。
コンティグ0009の長さは286369でしたので、nuc9における2列目の最大値がそれに近いことを確認します。
max(nuc9[ ,2])
により、286351となり、問題ないことを確認します。
コンティグ0009にヌクレオソームマップされた数は
nrow(nuc9)
により、656815であることより、塩基当りのマップ数が2を超えていることがわかります。
コンティグ0009のヌクレオソームDNAの長さの分布を見てみます。
hist(nuc9[,2]-nuc9[,1]+1)
boxplot(nuc9[,2]-nuc9[,1]+1)
mean(nuc9[,2]-nuc9[,1]+1)
median(nuc9[,2]-nuc9[,1]+1)
次に、個々のヌクレオソームDNA断片に位置情報を得るために、列1と列2からマップされた中央位置を算出します。
pos9 <- (nuc9[,1] + nuc9[,2])/2
length(pos9)が656815であることを確認してください。
ヌクレオソームの分布を見るための方法はいくつかありますが、データが大きいため、Rで全体を見ることはあきらめてください。
ただ、部分を見ることはできます。例えば、1000塩基位置から2000塩基位置のヌクレオソームの頻度を見たい場合、
x <- 1000
y <- 2000
a1 <- nuc9[nuc9[,1] >= x & nuc9[,2] <= y, ]
データを整理します。
a2 <- data.frame(a1[,1]-x+1, a1[,2]-x+1)
a3 <- y-x+1
b <- numeric(nrow(a2)*a3)
for(i in 1:nrow(a2)){
b[seq(a3*i+a2[i,1]-a3, a3*i+a2[i,2]-a3)] <- 1
}
b1 <- matrix(b, nrow=nrow(a2), byrow=T)
b2 <- apply(b1, 2, sum)
プロファイルを確認します。
plot(b2, type="h")
として、その様子を見ることはできます。
3. 各遺伝子領域におけるヌクレオソーム数のカウント
生物統計解析基礎の9. 行列とデータフレームにおいて作成したgene.txtを使って、それぞれの遺伝子領域にマップされたヌクレオソームDNA断片の数をカウントします。
gene <- read.table("gene.txt")
nuc <- c( )
for(i in 1:nrow(gene)) {
nuc[i] <- length(pos9[pos9 > min(gene[i, 1], gene[i, 2]) & pos9 < max(gene[i, 1], gene[i,2])])
}
ここでgeneの位置における最小値と最大値を持ち出してる理由を考えてください。
このgeneの位置は、遺伝子が発現する際の方向の情報も入っているということです。
4. GC含量とヌクレオソーム密度の関係
生物統計解析基礎の9. 行列とデータフレームにおいて計算した各遺伝子のGC含量と各遺伝子におけるヌクレオソーム密度を比較します。
gc <- scan("gc.txt")
size <- c( )
for(i in 1:nrow(gene)) size[i] <- length(gc[gene[i,1]:gene[i,2]])
GC <- c( )
for(i in 1:nrow(gene)) GC[i] <- sum(gc[gene[i,1]:gene[i,2]])
plot(GC/size, nuc/size) #この様子を見て気づくことはありますか？
cor(GC/size, nuc/size) #相関係数、この結果から何か言えるか考えてください
この関係についてまとめたものが下記の論文です。
Nishida H et al. (2013) Base composition and nucleosome density in exonic and intronic regions in genes of the filamentous ascomycetes Aspergillus nidulans and Aspergillus oryzae. Gene 525, 5-10
また、ヌクレオソームの位置と遺伝子発現の情報を組み合わせることによって、下記のような結果を得ることできます。
ここで強調したいことは、「生物統計解析基礎」をはじめ、microbioinformaticsで示したものは、もちろん、この「ヌクレオソーム解析」を含め、すべてRの基本パッケージであるbaseしか使っていません（それにこだわって書いてきました）。
Rには、様々なパッケージ・ライブラリーが用意されており、それが魅力的なところの１つですが、それを使わなくてもこの程度はできるということを知っていただき、さらにそれ以上を目指す人は、是非そのように進めてほしいと思います。
------------------------------------------------------------------------------------------------
【学生実験予備】
0. はじめに
真核細胞では、DNAはヒストン８量体に結合してヌクレオソーム構造をとっています。
ヌクレオソームの形成位置が厳密に決まっているわけではないのですが、密なところと疎のところがあります。
ヌクレオソーム位置やその位置の決定方法などについては、インターネットなどを使って調べてください。
ここでは、詳細に知る必要はありませんが、ヒストンとDNAの「タンパク質・DNA相互作用」の一つです。
ヒストンがゲノムにおけるどの位置に結合しているか、というテーマなので、ゲノム上のどこからどこまでに結合していたという情報が必要です。
ヒストン８量体（H2A, H2B, H3, H4が２セット）にDNAは約２周巻き付いています。
ここでは、ヌクレオソーム形成が遺伝子領域毎に違っているかどうかを調べ、GC含量との関連についても調べることを目的とします。
Rを使用します。
Rを立ち上げてください。
「ファイル」をクリックして、「ディレクトリの変更」をクリックし、デスクトップを選択し、OKとしてください。
次に、nuc9.txtのファイルをダウンロードしてください。
ダウンロードが完了後、そのファイルをデスクトップに置いてください。
Rにおいて、
nuc9 <- read.table("nuc9.txt")
を実行してください。
次に、メールで添付したファイルを開けて、すべてを選択し、Rにおいて、ペーストしてください。
これによって、
nuc9
gene
size
GC
の４つのファイルがRにおいて作られました。
1. nuc9の確認
このファイルは1列目にヒストン８量体に結合していたDNA断片のコンティグ配列上での最初の位置、2列目に最後の位置を示しています。
コンティグ0009の長さは286369でしたので、nuc9における2列目の最大値がそれに近いことを確認します。
max(nuc9[ ,2])
により、286351となり、問題ないことを確認します。
コンティグ0009にヌクレオソームマップされた数は
nrow(nuc9)
により、656815であることより、塩基当りのマップ数が2を超えていることがわかります。
コンティグ0009のヌクレオソームDNAの長さの分布を見てみます。
hist(nuc9[,2]-nuc9[,1]+1)
boxplot(nuc9[,2]-nuc9[,1]+1)
mean(nuc9[,2]-nuc9[,1]+1)
median(nuc9[,2]-nuc9[,1]+1)
次に、個々のヌクレオソームDNA断片に位置情報を得るために、列1と列2からマップされた中央位置を算出します。
pos9 <- (nuc9[,1] + nuc9[,2])/2
length(pos9)
が656815であることを確認してください。
ヌクレオソームの分布を見るための方法はいくつかありますが、データが大きいため、Rで全体を見ることはあきらめてください。ただ、部分を見ることはできます。例えば、
plot(sort(pos9)[1:1000], rep(1, 1000))
2. 各遺伝子領域におけるヌクレオソーム数のカウント
geneは２列からなっており、遺伝子の両端の位置を示しています。
nuc <- c( )
for(i in 1:nrow(gene)) {
nuc[i] <- length(pos9[pos9 > min(gene[i, 1], gene[i, 2]) & pos9 < max(gene[i, 1], gene[i,2])])
}
ここでgeneの位置における最小値と最大値を持ち出してる理由は、geneにおける列１は開始点、列２は終結点を示しているからです（遺伝子発現の向きの情報も入っているということです）。
3. GC含量とヌクレオソーム密度の関係
GCは各遺伝子領域におけるGとCの数を示しています。
よって、各遺伝子領域のGC含量は、GC/sizeで計算できます。
plot(GC/size, nuc/size) #GC含量とヌクレオソーム密度の関係
cor(GC/size, nuc/size) #相関係数
これらと
plot(GC/size, size) #GC含量と遺伝子サイズの関係
plot(size, nuc/size) #遺伝子サイズとヌクレオソーム密度の関係
を比較してください。
課題： plot(GC/size, nuc/size)、plot(GC/size, size)、plot(size, nuc/size)の図を送ってください。また、cor(GC/size, nuc/size)、cor(GC/size, size)、cor(size, nuc/size)の数値を示し、どのようなことが言えるかについて述べてください。