> x <- c(1:5)
> y <- c(1:3)
> y[4:5] <- NA
> y
[1] 1 2 3 NA NA
> z <- cbind(x, y)
> z
x y
[1,] 1 1
[2,] 2 2
[3,] 3 3
[4,] 4 NA
[5,] 5 NA
> boxplot(x,y)
医療系の仕事をしています。生命の尊さ、美しさがどのようなメカニズムで生じるのかに興味があります。科学の方法論を用いて、このような問いに応えたい、私はこう思って医学生物学の基礎研究のトレーニングを受けてきました。生命を科学的手法を用いて理解を試みる上で、genomeを始めとした種々の大量データの処理が必要不可欠であることを痛感しました。また、生命科学が物理学、数学、統計学、有機化学などの種々の学問と深い関わりを持つことを実感しました。そのため、このブログは広範囲の学問領域に関しての記事を載せています。日々の学習内容を文書に書き残し、それを読み返すことによって、体系化された知識を身に付けることを目標としています。どうぞよろしくお願いします。
美しい散布図の描き方
source("http://aoki2.si.gunma-u.ac.jp/R/src/dot_plot.R", encoding="euc-jp")
a <- rnorm(10, mean = 3)
b <- rnorm(10, mean = 9)
area <- c(a,b)
fac <- c(rep("control", 10), rep("treat", 10))
factor <- as.factor(fac)
png("area.png")
dot.plot(factor, area, accu=0.1, stp=0.5, xlab="Group", ylab="Area", main="Staining Area")
dev.off()
a <- rnorm(10, mean = 3)
b <- rnorm(10, mean = 9)
area <- c(a,b)
fac <- c(rep("control", 10), rep("treat", 10))
factor <- as.factor(fac)
png("area.png")
dot.plot(factor, area, accu=0.1, stp=0.5, xlab="Group", ylab="Area", main="Staining Area")
dev.off()
shellでもテキストからでもHello,worldをRで!
(1)シェルで
$ R
> print("Hello, world!")
$ R
> print("Hello, world!")
[1] "Hello, world!"
> quit()
(2)テキストで(バッチ処理)
$ vim hello.R
a
print("Hello, world!")
Esc
:wq
$ R CMD BATCH hello.R output.txt
1
2 R version 2.10.1 (2009-12-14)
3 Copyright (C) 2009 The R Foundation for Statistical Computing
4 ISBN 3-900051-07-0
5
6 Rは、自由なソフトウェアであり、「完全に無保証」です。
7 一定の条件に従えば、自由にこれを再配布することができます。
8 配布条件の詳細に関しては、'license()'あるいは'licence()'と入力してください。
9
10 Rは多くの貢献者による共同プロジェクトです。
11 詳しくは'contributors()'と入力してください。
12 また、RやRのパッケージを出版物で引用する際の形式については
13 'citation()'と入力してください。
14
15 'demo()'と入力すればデモをみることができます。
16 'help()'とすればオンラインヘルプが出ます。
17 'help.start()'でHTMLブラウザによるヘルプがみられます。
18 'q()'と入力すればRを終了します。
19
20 > print("Hello, world!")
21 [1] "Hello, world!"
22 >
23 >
24 > proc.time()
25 ユーザ システム 経過
26 0.456 0.012 0.553
:q
shellでもテキストからでもHello,worldをrubyで!
(1)テキストで
$ vim hello.rb
a
1 print "Hello, world\n!";
Esc
:wq
$ ruby hello.rb
Hello, worl!
(2)シェルで
$ ruby <<END
> print "hello\n";
> END
hello
$ vim hello.rb
a
1 print "Hello, world\n!";
Esc
:wq
$ ruby hello.rb
Hello, worl!
(2)シェルで
$ ruby <<END
> print "hello\n";
> END
hello
shellでもテキストからでもHello,worldをなperl
1)まずはテキストで編集して実行
$ vim hello.pl
a
1 use strict;
2 use warnings;
3
4 print "Hello, world!\n";
5
6 exit;
$ vim hello.pl
a
1 use strict;
2 use warnings;
3
4 print "Hello, world!\n";
5
6 exit;
Esc
:wq
$ perl hello.pl
Hello, world!
1)次はシェルの中で打ち込んで。
$ perl -e "print 'Hello, world'"
Hello, world
いい感じ♪
シェルでhtmlを作成した後、firefoxで内容を表示してみる♪
HTMLファイルを作ります。
$ vim test.html
a
test
Esc
:wq
$ firefox test.html
これでブラウザが開いて、test.htmの内容lがウェブブラウザの表示形式で確認できます( ^ ^ )
ububutuにgtk+を突っ込む。
$ sudo aptitude install libgtk2.0-dev
でいいけます。knoppix for bioにはこのコノコマンドでは入りませんでした。依存関係が複雑過ぎて、やる気が失せました。
ちなみにgtk+を使ったコンパイルは、たとえば「test.c」というファイルに対して行うのであれば、
$ (pkg-config --cflags --libs gtk+-2.0) test.c -o test
と行えばよい。
javaでHello!
javaに手を出してみました.....
$ vim hello.java
a
$ vim hello.java
a
public class hello {
public static void main (String[] args) {
System.out.println("Hello!");
}
}
Esc
:wq
$ javac hello.java
$ ls
hello.class hello.java
$ java hello
Hello!
♪( ̄∠  ̄ )ノ
gccでmath.hのヘッダーファイルがリンクしない(涙)
どうやらmath.hは巨大であり、かつ一部のユーザーしか使わないために、デフォルトではリンクしないようなのである。
以下のように -lm というオプションをつけることでリンクされる。
例)
$ vim math.c
a
#include <stdio.h>
#include <math.h>
int main(void)
{
double suti, hieo;
printf("実数を入力してください");
scanf("%lf", &suti);
printf("%lfの平方根は%lfです\n", suti, heiho);
return o;
}
Esc
$ gcc -lf -o math math.c
$ ./math
実数を入力してください9
9.000000の平方根は3.00000です。
以下のように -lm というオプションをつけることでリンクされる。
例)
$ vim math.c
a
#include <stdio.h>
#include <math.h>
int main(void)
{
double suti, hieo;
printf("実数を入力してください");
scanf("%lf", &suti);
printf("%lfの平方根は%lfです\n", suti, heiho);
return o;
}
Esc
$ gcc -lf -o math math.c
$ ./math
実数を入力してください9
9.000000の平方根は3.00000です。
マイホームページが完成したぞいヘ(*゚∇゚)ノ
Bioinformaticsのお勉強(ブログと名前は一緒♪)
http://kappabucha.web.fc2.com/
よろしくお願いしますm(_ _)m
http://kappabucha.web.fc2.com/
よろしくお願いしますm(_ _)m
FC2プロバイダーを用いたホームページの作成
1.テキストエディタ等のテキスト編集ソフトでhtmlファイルを作成
2.FC2(http://fc2.com/)でユーザー登録
3.FFFTPをダウンロード&インストール(http://www2.biglobe.ne.jp/~sota/ffftp.html)
4.FFFTPにより、htmlファイルをアップロード
2.FC2(http://fc2.com/)でユーザー登録
3.FFFTPをダウンロード&インストール(http://www2.biglobe.ne.jp/~sota/ffftp.html)
4.FFFTPにより、htmlファイルをアップロード
knoppixで、usbのマウント、アンマウントのコマンド
usbをマウントするときは、
$ su
# cd /mnt
# mkdir usbmem
# mount -t vfat /dev/sda1 /mnt/usbmem
# mount これでマウント状況を確認できる。
usbをアマウント(これをしないと書き換えが正常に反映されない)
# umount /mnt/usbmem
最後に、/etc/fstab に設定を書き入れてオプションせずにマウントできるようにする。
# vim /etc/fstab
で /etc/fstabを開いて、一番最後の行に
/dev/sda1 /mnt/usbfm auto noauto,user 0 0
と書き入れる。
これで以下のようにしてusbのマウントができる。
# mount /mnt/usbmem
これで確認せよ。
# mount
$ su
# cd /mnt
# mkdir usbmem
# mount -t vfat /dev/sda1 /mnt/usbmem
# mount これでマウント状況を確認できる。
usbをアマウント(これをしないと書き換えが正常に反映されない)
# umount /mnt/usbmem
最後に、/etc/fstab に設定を書き入れてオプションせずにマウントできるようにする。
# vim /etc/fstab
で /etc/fstabを開いて、一番最後の行に
/dev/sda1 /mnt/usbfm auto noauto,user 0 0
と書き入れる。
これで以下のようにしてusbのマウントができる。
# mount /mnt/usbmem
これで確認せよ。
# mount
Linuxのシェルでコマンドを保存する方法(^-^)
$ R | tee -a test.log
> q()
Save workspace image? [y/n/c]: y
これでOKです♪
カレントディレクトリに、test.logというファイルが生成している。
> q()
Save workspace image? [y/n/c]: y
これでOKです♪
カレントディレクトリに、test.logというファイルが生成している。
KNOBがUSBを認識しない(涙)
ディバイスのマウントという作業が必要名ようである。
$ cd /mnt
$ mkdir usbmem
$ mount -t vfat /dev/sda1 /mnt/usbmem
Windowsや最近のLinuxディストリビューションはマウントなんたる作業は必要ない。
なので、KNOB(knoppix)なんか使うと、想定外のことが発生して楽しいな♪
$ cd /mnt
$ mkdir usbmem
$ mount -t vfat /dev/sda1 /mnt/usbmem
Windowsや最近のLinuxディストリビューションはマウントなんたる作業は必要ない。
なので、KNOB(knoppix)なんか使うと、想定外のことが発生して楽しいな♪
KNOBがネットにつながらない(涙)
結局、解決しました。
ペンギンのアイコン→NetWork/Internet→ネットワークカードの設定
でやりますと・・・・
「su がエラーを返しました」
と悲しいお知らせ。
そこで、
$ vim /home/(ユーザー名)/.kde/share/config/kdesurc
とやって、vimで設定ファイルを開いて、
super-user-command=sudo
の一文を
super-user-command=su
と書き換える。
これで再度
ペンギンのアイコン→NetWork/Internet→ネットワークカードの設定
とやるとうまくいく。
あとは指示にしたがって、アドレスを取得すれば完了♪♪
ヾ(〃⌒ ー———⌒)ノ~~
ペンギンのアイコン→NetWork/Internet→ネットワークカードの設定
でやりますと・・・・
「su がエラーを返しました」
と悲しいお知らせ。
そこで、
$ vim /home/(ユーザー名)/.kde/share/config/kdesurc
とやって、vimで設定ファイルを開いて、
super-user-command=sudo
の一文を
super-user-command=su
と書き換える。
これで再度
ペンギンのアイコン→NetWork/Internet→ネットワークカードの設定
とやるとうまくいく。
あとは指示にしたがって、アドレスを取得すれば完了♪♪
ヾ(〃⌒ ー———⌒)ノ~~
knoppix for bioをHDDにインストールする方法を見つけた!!!!
ルートになってから実行します。
$ su
# knx2hd
ちなみにKNOB自体は
「オープンソースで学ぶバイオインフォマティクス」
の付属のDVDより。
$ su
# knx2hd
ちなみにKNOB自体は
「オープンソースで学ぶバイオインフォマティクス」
の付属のDVDより。
分散・平均・正規分布の検討
$ R
cilia <- rnorm(50)
range(cilia)
max(cilia)
min(cilia)
meancilia1 <- sum(cilia)/length(cilia)
meancilia2 <- mean(cilia)
median(cilia)
png("110317.figure1.png")
plot(cilia)
lines(c(1,20),c(min(cilia),min(cilia)),lty = 2)
lines(c(1,20),c(max(cilia),max(cilia)),lty = 2)
lines(c(10,10), c(min(cilia),max(cilia)))
text(10,max(cilia)-1,"Range = min(cilia) to max(cilia)")
dev.off()
png("110317.figure2.png")
plot(cilia)
abline(h= mean(cilia))
for(i in 1:length(cilia)){lines(c(i,i),c(mean(cilia),cilia[i]))}
dev.off()
SUM <- sum((cilia-mean(cilia))^2)
LENGTH <- length(cilia)
VALIANCE1 <- SUM/(LENGTH-1)
VALIANCE2 <- var(cilia)
VALIANCE3 <- (sum(cilia^2)-sum(cilia)^2/length(cilia))/(length(cilia)-1)
VALIANCE1
VALIANCE2
VALIANCE3
png("110317_figure3.png")
par(mfrow = c(1,2))
plot(rnorm(1000),rnorm(1000), ylab = "Label for y axis",xlab = "Label for x axis")
plot(rnorm(1000),rnorm(1000), ylab = "Label for y axis",xlab = "Label for x axis",las = 1 , cex.lab = 1.5, main = "rnorm plot")
dev.off()
x <- seq(-4,4, len =101)
plot(x, sin(x))
plot(x, sin(x), type ="l", xaxt = "n", xlab = expression(pasete("Phase Angle",phi)), ylab = expression("sin"phi))
kansu <- function(x){
maximum <- max(x)
minimum <- min(x)
mean <- mean(x)
median <- median(x)
length <- length(x)
cat("Maximum", maximum, "\n")
cat("Minimum", minimum, "\n")
cat("Mean", mean, "\n")
cat("Median", median, "\n")
cat("Length",length, "\n")
}
x <- rnorm(100)
kansu(x)
plot(c(0,32), c(0,15), type="n", xlab = "Sample size", ylab = "Variance")
for(df in seq(3,31,2)){
for(i in 1:30){
x <- rnorm(df, mean =10, sd =2)
points(df, var(x))
}
}
png("110317_figure4.png")
plot(c(0,32), c(0,15), type="n", xlab = "Sample size", ylab = "Variance")
for(df in seq(3,31,2)){
for(i in 1:300){
x <- rnorm(df, mean =10, sd =2)
points(df, var(x))
}
}
dev.off()
png("110317_figure6.png")
par(mfrow = c(1,5))
for(i in 1:5){
plot(runif(10^i, min = 0, max = 100))
}
dev.off()
png("110317_figure7.png")
means <- numeric(10000)
for (i in 1:10000){
means[i] <- mean(runif(5, min = 0, max = 10))
}
hist(means)
dev.off()
means <- numeric(10000)
for (i in 1:10000){
means[i] <- mean(runif(5, min=0, max=10))
}
hist(means)
yv <- dnorm(seq(0,10,0.1), mean = mean(means), sd = sd(means))*5000
line(seq(0,10,0.1), yv)
cilia <- rnorm(50)
range(cilia)
max(cilia)
min(cilia)
meancilia1 <- sum(cilia)/length(cilia)
meancilia2 <- mean(cilia)
median(cilia)
png("110317.figure1.png")
plot(cilia)
lines(c(1,20),c(min(cilia),min(cilia)),lty = 2)
lines(c(1,20),c(max(cilia),max(cilia)),lty = 2)
lines(c(10,10), c(min(cilia),max(cilia)))
text(10,max(cilia)-1,"Range = min(cilia) to max(cilia)")
dev.off()
png("110317.figure2.png")
plot(cilia)
abline(h= mean(cilia))
for(i in 1:length(cilia)){lines(c(i,i),c(mean(cilia),cilia[i]))}
dev.off()
SUM <- sum((cilia-mean(cilia))^2)
LENGTH <- length(cilia)
VALIANCE1 <- SUM/(LENGTH-1)
VALIANCE2 <- var(cilia)
VALIANCE3 <- (sum(cilia^2)-sum(cilia)^2/length(cilia))/(length(cilia)-1)
VALIANCE1
VALIANCE2
VALIANCE3
png("110317_figure3.png")
par(mfrow = c(1,2))
plot(rnorm(1000),rnorm(1000), ylab = "Label for y axis",xlab = "Label for x axis")
plot(rnorm(1000),rnorm(1000), ylab = "Label for y axis",xlab = "Label for x axis",las = 1 , cex.lab = 1.5, main = "rnorm plot")
dev.off()
x <- seq(-4,4, len =101)
plot(x, sin(x))
plot(x, sin(x), type ="l", xaxt = "n", xlab = expression(pasete("Phase Angle",phi)), ylab = expression("sin"phi))
kansu <- function(x){
maximum <- max(x)
minimum <- min(x)
mean <- mean(x)
median <- median(x)
length <- length(x)
cat("Maximum", maximum, "\n")
cat("Minimum", minimum, "\n")
cat("Mean", mean, "\n")
cat("Median", median, "\n")
cat("Length",length, "\n")
}
x <- rnorm(100)
kansu(x)
plot(c(0,32), c(0,15), type="n", xlab = "Sample size", ylab = "Variance")
for(df in seq(3,31,2)){
for(i in 1:30){
x <- rnorm(df, mean =10, sd =2)
points(df, var(x))
}
}
png("110317_figure4.png")
plot(c(0,32), c(0,15), type="n", xlab = "Sample size", ylab = "Variance")
for(df in seq(3,31,2)){
for(i in 1:300){
x <- rnorm(df, mean =10, sd =2)
points(df, var(x))
}
}
dev.off()
png("110317_figure6.png")
par(mfrow = c(1,5))
for(i in 1:5){
plot(runif(10^i, min = 0, max = 100))
}
dev.off()
png("110317_figure7.png")
means <- numeric(10000)
for (i in 1:10000){
means[i] <- mean(runif(5, min = 0, max = 10))
}
hist(means)
dev.off()
means <- numeric(10000)
for (i in 1:10000){
means[i] <- mean(runif(5, min=0, max=10))
}
hist(means)
yv <- dnorm(seq(0,10,0.1), mean = mean(means), sd = sd(means))*5000
line(seq(0,10,0.1), yv)
Wilcoxon検定による二群比較
#データの読み込み
> f.test.data <- read.table("/home/kappa/デスクトップ/kappa/2011/1103/統計学:Rを用いた入門書/f.test.data.txt", header = T)
> f.test.data
gardenB gardenC
1 5 3
2 5 3
3 6 2
4 7 1
5 4 10
6 4 4
7 3 3
8 5 11
9 6 3
10 5 10
> attach(f.test.data)
#標本データの概観
> png("110318.figure3.png")
> par(mfrow = c(1,2))
> hist(gardenB)
> hist(gardenc)
> hist(gardenC)
> dev.off()
> png("110318.figure4.png")
> boxplot(gardenB, main = "gardenB")
> boxplot(gardenC, main = "gardenC")
> dev.off()
null device
1
#正規性の検定
> png("110318.figure1.png")
> par(mfrow = c(1,2))
> qqnorm(gardenB, main = "gardenB")
> qqline(gardenB)
> qqnorm(gardenC, main = "gardenC")
> qqline(gardenC)
> dev.off()
> shapiro.test(gardenB)
Shapiro-Wilk normality test
data: gardenB
W = 0.9529, p-value = 0.7026
> sahpiro.test(gardenC)
Shapiro-Wilk normality test
data: gardenC
W = 0.78, p-value = 0.008284
#Wilcoxonの順位輪検定(gardenCの標本は正規性を示さないので)
> ozone <- c(gardenB, gardenC)
> ozone
[1] 5 5 6 7 4 4 3 5 6 5 3 3 2 1 10 4 3 11 3 10
> label <- c(rep("B", length(gardenB)), rep("C", length(gardenC)))
> label
[1] "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "C" "C" "C" "C" "C" "C" "C" "C" "C"
[20] "C"
> combined.ranks <- rank(ozone)
> combined.ranks
[1] 12.5 12.5 15.5 17.0 9.0 9.0 5.0 12.5 15.5 12.5 5.0 5.0 2.0 1.0 18.5
[16] 9.0 5.0 20.0 5.0 18.5
> tapply(combined.ranks, label, sum)
B C
121 89
#えられた2つの値の最小値89をウィルコクソン順位和の数表にある数値と比較する。標本吸うが10,10の場合は5%値は78である。したがって、二つの標本平均には有意差が認められない。
> wilcox.test(gardenB, gardenC)
Wilcoxon rank sum test with continuity correction
data: gardenB and gardenC
W = 66, p-value = 0.2349
alternative hypothesis: true location shift is not equal to 0
警告メッセージ:
In wilcox.test.default(gardenB, gardenC) :
タイがあるため、正確な p 値を計算することができません
> f.test.data <- read.table("/home/kappa/デスクトップ/kappa/2011/1103/統計学:Rを用いた入門書/f.test.data.txt", header = T)
> f.test.data
gardenB gardenC
1 5 3
2 5 3
3 6 2
4 7 1
5 4 10
6 4 4
7 3 3
8 5 11
9 6 3
10 5 10
> attach(f.test.data)
#標本データの概観
> png("110318.figure3.png")
> par(mfrow = c(1,2))
> hist(gardenB)
> hist(gardenc)
> hist(gardenC)
> dev.off()
> png("110318.figure4.png")
> boxplot(gardenB, main = "gardenB")
> boxplot(gardenC, main = "gardenC")
> dev.off()
null device
1
#正規性の検定
> png("110318.figure1.png")
> par(mfrow = c(1,2))
> qqnorm(gardenB, main = "gardenB")
> qqline(gardenB)
> qqnorm(gardenC, main = "gardenC")
> qqline(gardenC)
> dev.off()
> shapiro.test(gardenB)
Shapiro-Wilk normality test
data: gardenB
W = 0.9529, p-value = 0.7026
> sahpiro.test(gardenC)
Shapiro-Wilk normality test
data: gardenC
W = 0.78, p-value = 0.008284
#Wilcoxonの順位輪検定(gardenCの標本は正規性を示さないので)
> ozone <- c(gardenB, gardenC)
> ozone
[1] 5 5 6 7 4 4 3 5 6 5 3 3 2 1 10 4 3 11 3 10
> label <- c(rep("B", length(gardenB)), rep("C", length(gardenC)))
> label
[1] "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "C" "C" "C" "C" "C" "C" "C" "C" "C"
[20] "C"
> combined.ranks <- rank(ozone)
> combined.ranks
[1] 12.5 12.5 15.5 17.0 9.0 9.0 5.0 12.5 15.5 12.5 5.0 5.0 2.0 1.0 18.5
[16] 9.0 5.0 20.0 5.0 18.5
> tapply(combined.ranks, label, sum)
B C
121 89
#えられた2つの値の最小値89をウィルコクソン順位和の数表にある数値と比較する。標本吸うが10,10の場合は5%値は78である。したがって、二つの標本平均には有意差が認められない。
> wilcox.test(gardenB, gardenC)
Wilcoxon rank sum test with continuity correction
data: gardenB and gardenC
W = 66, p-value = 0.2349
alternative hypothesis: true location shift is not equal to 0
警告メッセージ:
In wilcox.test.default(gardenB, gardenC) :
タイがあるため、正確な p 値を計算することができません
Student's t-test による二群間比較
#データの読み込み
> t.test.data <- read.table("/home/kappa/デスクトップ/kappa/2011/1103/統計学:Rを用いた入門書/t.test.data.txt", header = T)
#データの概観
> t.test.data
gardenA gardenB
1 3 5
2 4 5
3 4 6
4 3 7
5 2 4
6 3 4
7 1 3
8 3 5
9 5 6
10 2 5
par(mfrow =c(1,2))
plot(gardenA, main = "gardenA")
abline( h = mean(gardenA))
for ( i in 1:length(gardenA)){
lines(c(i,i), c(mean(gardenA), gardenA[i]))
}
plot(gardenB, main = "gardenB")
abline( h = mean(gardenB))
for ( i in 1:length(gardenB)){
lines(c(i,i), c(mean(gardenB), gardenB[i]))
}
> par(mfrow = c(1,3))
> attach(t.test.data)
> hist(gardenA)
> hist(gardenB)
> boxplot(t.test.data)
#正規性の検定
> qqnorm(gardenA)
> qqline(gardenA)
> qqnorm(gardenB)
> qqline(gardenB)
> shapiro.test(gardenA)
Shapiro-Wilk normality test
data: gardenA
W = 0.9529, p-value = 0.7026
> shapiro.test(gardenB)
Shapiro-Wilk normality test
data: gardenB
#等分散の検定
> 2 * (1 - pf(var(gardenA)/var(gardenB), length(gardenA)-1, length(gardenB)-1)) > 0.05
[1] TRUE
> var.test(gardenA, gardenB)
F test to compare two variances
data: gardenA and gardenB
F = 1, num df = 9, denom df = 9, p-value = 1
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.2483859 4.0259942
sample estimates:
ratio of variances
1
#Student t検定(gardenA, gardenBは正規分布。さらに、二つの分散は等しいため)
> t.test(gardenA, gardenB, var.equal = T)
Two Sample t-test
data: gardenA and gardenB
t = -3.873, df = 18, p-value = 0.001115
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3.0849115 -0.9150885
sample estimates:
mean of x mean of y
3 5
#ちなみに分散が等しくないときはWelchのt検定
> t.test(gardenA, gardenB, var.equal = F)
Welch Two Sample t-test
data: gardenA and gardenB
t = -3.873, df = 18, p-value = 0.001115
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3.0849115 -0.9150885
sample estimates:
mean of x mean of y
3 5
#ちなみにこのサンプルにWilicoxonを用いると
> wilcox.test(gardenA, gardenB)
Wilcoxon rank sum test with continuity correction
data: gardenA and gardenB
W = 11, p-value = 0.002988
alternative hypothesis: true location shift is not equal to 0
警告メッセージ:
In wilcox.test.default(gardenA, gardenB) :
タイがあるため、正確な p 値を計算することができません
#Wilcoxon検定のp値の方がt検定のp値よりも大きい。つまり、Wilcoxonの方が保守的であるといえる。"誤差が正規的でないとき"、ノンパラメトリック検定はt検定よりも優れて適切な検定である。
> t.test.data <- read.table("/home/kappa/デスクトップ/kappa/2011/1103/統計学:Rを用いた入門書/t.test.data.txt", header = T)
#データの概観
> t.test.data
gardenA gardenB
1 3 5
2 4 5
3 4 6
4 3 7
5 2 4
6 3 4
7 1 3
8 3 5
9 5 6
10 2 5
par(mfrow =c(1,2))
plot(gardenA, main = "gardenA")
abline( h = mean(gardenA))
for ( i in 1:length(gardenA)){
lines(c(i,i), c(mean(gardenA), gardenA[i]))
}
plot(gardenB, main = "gardenB")
abline( h = mean(gardenB))
for ( i in 1:length(gardenB)){
lines(c(i,i), c(mean(gardenB), gardenB[i]))
}
> par(mfrow = c(1,3))
> attach(t.test.data)
> hist(gardenA)
> hist(gardenB)
> boxplot(t.test.data)
#正規性の検定
> qqnorm(gardenA)
> qqline(gardenA)
> qqnorm(gardenB)
> qqline(gardenB)
> shapiro.test(gardenA)
Shapiro-Wilk normality test
data: gardenA
W = 0.9529, p-value = 0.7026
> shapiro.test(gardenB)
Shapiro-Wilk normality test
data: gardenB
#等分散の検定
> 2 * (1 - pf(var(gardenA)/var(gardenB), length(gardenA)-1, length(gardenB)-1)) > 0.05
[1] TRUE
> var.test(gardenA, gardenB)
F test to compare two variances
data: gardenA and gardenB
F = 1, num df = 9, denom df = 9, p-value = 1
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.2483859 4.0259942
sample estimates:
ratio of variances
1
#Student t検定(gardenA, gardenBは正規分布。さらに、二つの分散は等しいため)
> t.test(gardenA, gardenB, var.equal = T)
Two Sample t-test
data: gardenA and gardenB
t = -3.873, df = 18, p-value = 0.001115
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3.0849115 -0.9150885
sample estimates:
mean of x mean of y
3 5
#ちなみに分散が等しくないときはWelchのt検定
> t.test(gardenA, gardenB, var.equal = F)
Welch Two Sample t-test
data: gardenA and gardenB
t = -3.873, df = 18, p-value = 0.001115
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3.0849115 -0.9150885
sample estimates:
mean of x mean of y
3 5
#ちなみにこのサンプルにWilicoxonを用いると
> wilcox.test(gardenA, gardenB)
Wilcoxon rank sum test with continuity correction
data: gardenA and gardenB
W = 11, p-value = 0.002988
alternative hypothesis: true location shift is not equal to 0
警告メッセージ:
In wilcox.test.default(gardenA, gardenB) :
タイがあるため、正確な p 値を計算することができません
#Wilcoxon検定のp値の方がt検定のp値よりも大きい。つまり、Wilcoxonの方が保守的であるといえる。"誤差が正規的でないとき"、ノンパラメトリック検定はt検定よりも優れて適切な検定である。
機械的に発生させた正規乱数に対してSAMを用いてみたシミュレーション
cnt1 <- rnorm(40000, mean = 200, sd = 1)
cnt2 <- rnorm(40000, mean = 255, sd = 1)
cnt3 <- rnorm(40000, mean = 230, sd = 1)
cnt4 <- rnorm(40000, mean = 344, sd = 1)
treat1 <- rnorm(40000, mean = 250, sd = 1)
treat2 <- rnorm(40000, mean = 210, sd = 1)
treat3 <- rnorm(40000, mean = 240, sd = 1)
treat4 <- rnorm(40000, mean = 210, sd = 1)
data <- cbind(cnt1, cnt2, cnt3, cnt4, treat1,treat2, treat3, treat4)
data[1,]
cnt1 cnt2 cnt3 cnt4 treat1 treat2 treat3 treat4
199.8450 255.8026 230.7610 343.8670 248.4101 210.0879 239.5559 211.8849
boxplot(data)
par(mfrow = c(2,4))
hist(data[,1], main = "control 1",xlab = "gene expression")
hist(data[,2], main = "control 2",xlab = "gene expression")
hist(data[,3], main = "control 3",xlab = "gene expression")
hist(data[,4], main = "control 4",xlab = "gene expression")
hist(data[,5], main = "treatment 1",xlab = "gene expression")
hist(data[,6], main = "treatment 2",xlab = "gene expression")
hist(data[,7], main = "treatment 3",xlab = "gene expression")
hist(data[,8], main = "treatment 4",xlab = "gene expression")
dev.off()
for (i in 1:8){
data[,i] <- data[,i]/median(data[,i])
}
boxplot(data)
for (i in 1: length(data[,1])){
data[i,] <- data[i,]/median(data[i,])
}
for (i in 1: length(data[,1])){
data[i,] <- log(data[i,], 2)
}
norm <- log(data, 2)
library(siggenes)
n1 <- length(data[1,])/2
n2 <- length(data[1,])/2
cl <- rep(c(0,1), c(n1,n2))
sam.out <- sam(norm, cl, rand = 123)
summary(sam.out)
print(sam.out, seq(1,3,0.1)
plot(sam.out, 1)
plot(sam.out, 0,1)
plot(sam.out, 0.01)
print(sam.out, seq(0.01, 0.1, 0.01))
SAM Analysis for the Two-Class Unpaired Case Assuming Unequal Variances
Delta p0 False Called FDR
1 0.01 0.995 94.714 106 0.889
2 0.02 0.995 2.629 5 0.523
3 0.03 0.995 2.629 5 0.523
4 0.04 0.995 0 0 0
5 0.05 0.995 0 0 0
6 0.06 0.995 0 0 0
7 0.07 0.995 0 0 0
8 0.08 0.995 0 0 0
9 0.09 0.995 0 0 0
10 0.10 0.995 0 0 0
sum.sam.out <- summary(sam.out, 0.03)
sum.sam.out
slotNames(sum.sam.out)
[1] "row.sig.genes" "mat.fdr" "mat.sig" "list.args"
significant.genes <- sum.sam.out@mat.sig
> significant.genes[1:5,]
Row d.value stdev rawp q.value R.fold
1 766 2.067911 1.888463 1.285714e-05 0.5142857 23.17078675
2 29588 2.036899 1.519079 2.571429e-05 0.5142857 13.12129065
3 34635 -2.017873 1.497832 4.285714e-05 0.5714286 0.08042155
4 7750 -1.992394 1.562408 9.428571e-05 0.7836735 0.07593873
5 17922 -1.986813 1.423697 1.028571e-04 0.7836735 0.09258940
cnt2 <- rnorm(40000, mean = 255, sd = 1)
cnt3 <- rnorm(40000, mean = 230, sd = 1)
cnt4 <- rnorm(40000, mean = 344, sd = 1)
treat1 <- rnorm(40000, mean = 250, sd = 1)
treat2 <- rnorm(40000, mean = 210, sd = 1)
treat3 <- rnorm(40000, mean = 240, sd = 1)
treat4 <- rnorm(40000, mean = 210, sd = 1)
data <- cbind(cnt1, cnt2, cnt3, cnt4, treat1,treat2, treat3, treat4)
data[1,]
cnt1 cnt2 cnt3 cnt4 treat1 treat2 treat3 treat4
199.8450 255.8026 230.7610 343.8670 248.4101 210.0879 239.5559 211.8849
boxplot(data)
par(mfrow = c(2,4))
hist(data[,1], main = "control 1",xlab = "gene expression")
hist(data[,2], main = "control 2",xlab = "gene expression")
hist(data[,3], main = "control 3",xlab = "gene expression")
hist(data[,4], main = "control 4",xlab = "gene expression")
hist(data[,5], main = "treatment 1",xlab = "gene expression")
hist(data[,6], main = "treatment 2",xlab = "gene expression")
hist(data[,7], main = "treatment 3",xlab = "gene expression")
hist(data[,8], main = "treatment 4",xlab = "gene expression")
dev.off()
for (i in 1:8){
data[,i] <- data[,i]/median(data[,i])
}
boxplot(data)
for (i in 1: length(data[,1])){
data[i,] <- data[i,]/median(data[i,])
}
for (i in 1: length(data[,1])){
data[i,] <- log(data[i,], 2)
}
norm <- log(data, 2)
library(siggenes)
n1 <- length(data[1,])/2
n2 <- length(data[1,])/2
cl <- rep(c(0,1), c(n1,n2))
sam.out <- sam(norm, cl, rand = 123)
summary(sam.out)
print(sam.out, seq(1,3,0.1)
plot(sam.out, 1)
plot(sam.out, 0,1)
plot(sam.out, 0.01)
print(sam.out, seq(0.01, 0.1, 0.01))
SAM Analysis for the Two-Class Unpaired Case Assuming Unequal Variances
Delta p0 False Called FDR
1 0.01 0.995 94.714 106 0.889
2 0.02 0.995 2.629 5 0.523
3 0.03 0.995 2.629 5 0.523
4 0.04 0.995 0 0 0
5 0.05 0.995 0 0 0
6 0.06 0.995 0 0 0
7 0.07 0.995 0 0 0
8 0.08 0.995 0 0 0
9 0.09 0.995 0 0 0
10 0.10 0.995 0 0 0
sum.sam.out <- summary(sam.out, 0.03)
sum.sam.out
slotNames(sum.sam.out)
[1] "row.sig.genes" "mat.fdr" "mat.sig" "list.args"
significant.genes <- sum.sam.out@mat.sig
> significant.genes[1:5,]
Row d.value stdev rawp q.value R.fold
1 766 2.067911 1.888463 1.285714e-05 0.5142857 23.17078675
2 29588 2.036899 1.519079 2.571429e-05 0.5142857 13.12129065
3 34635 -2.017873 1.497832 4.285714e-05 0.5714286 0.08042155
4 7750 -1.992394 1.562408 9.428571e-05 0.7836735 0.07593873
5 17922 -1.986813 1.423697 1.028571e-04 0.7836735 0.09258940
C言語:入力した単語を文字列として認識させる。
$ vim 110310.c
a
1 #include <stdio.h>
2
3 int main(void)
4 {
5 char namae[10];
6 int nagasa;
7 int i;
8 int tuzuki;
9
10 tuzuki = 1;
11 printf("果物の名前入力>>");
12 for(nagasa = 0; tuzuki == 1; nagasa++)
13 {
14 scanf("%c", &namae[nagasa]);
15 if(namae[nagasa] == '\n')
16 {
17 tuzuki = 0;
18 }
19 }
20
21 nagasa--;
22
23 for(i = 0; i < nagasa; i++)
25 printf("%c", namae[i]);
26 }
27 printf("が大好き\n");
28 return 0;
29 }
Esc
$ gcc -o 110310 110310.c
$ ./110310
果物の名前入力>>ばなな
ばななが大好き
a
1 #include <stdio.h>
2
3 int main(void)
4 {
5 char namae[10];
6 int nagasa;
7 int i;
8 int tuzuki;
9
10 tuzuki = 1;
11 printf("果物の名前入力>>");
12 for(nagasa = 0; tuzuki == 1; nagasa++)
13 {
14 scanf("%c", &namae[nagasa]);
15 if(namae[nagasa] == '\n')
16 {
17 tuzuki = 0;
18 }
19 }
20
21 nagasa--;
22
23 for(i = 0; i < nagasa; i++)
25 printf("%c", namae[i]);
26 }
27 printf("が大好き\n");
28 return 0;
29 }
Esc
$ gcc -o 110310 110310.c
$ ./110310
果物の名前入力>>ばなな
ばななが大好き
入力した数の合計を計算させる☆
$ vim 110309.c
a
1 #include <stdio.h>
2
3 int main(void)
4 {
5 int ten[5];
6
7 int i;
8 int gokei;
9
10 for(i = 0; i < 5; i++)
11 {
12 printf("%d番目の点数の入力>>", i );
13 scanf("%d", &ten[i]);
14 }
15
16 gokei = 0;
17 for(i = 0; i < 5; i++)
18 {
19 gokei = gokei + ten[i];
20 }
21
22 printf("合計は%dです\n", gokei);
23 return 0;
24 }
Esc
:wq
$ gcc -o 110309 110309.c
$ ./110309
0番目の点数の入力>>2
1番目の点数の入力>>3
2番目の点数の入力>>5
3番目の点数の入力>>2
4番目の点数の入力>>5
合計は17です
$ vim 110309.2.c
a
1 #include <stdio.h>
2
3 int main(void)
4 {
5 int ten[5];
6
7 int i;
8 int gokei;
9
10 for(i = 0; i < 5; i++)
11 {
12 printf("%d番目の点数の入力>>", i );
13 scanf("%d", &ten[i]);
14 }
15
16 gokei = 0;
17 for(i = 0; i < 5; i++)
18 {
19 gokei = gokei + ten[i];
20 }
21
22 printf("合計は%dです\n", gokei);
23 return 0;
24 }
Esc
:wq
$ gcc -o 110309 110309.c
$ ./110309
0番目の点数の入力>>2
1番目の点数の入力>>3
2番目の点数の入力>>5
3番目の点数の入力>>2
4番目の点数の入力>>5
合計は17です
$ vim 110309.2.c
Cで配列の各要素にキーボードから入力
$ vim 110309.c
a
1 #include<stdio.h>
2
3 int main(void)
4 {
5 int hairetu[5];
6 int i;
7 for(i = 0; i < 5; i++)
8 {
9 printf("要素%dの入力<<", i);
10 scanf("%d", &hairetu[i]);
11 }
12 for(i = 0; i < 5; i++)
13 {
14 printf("要素%dは%dです\n", i, hairetu[i]);
15 }
16 return 0;
17 }
Esc
:wq
$ gcc -o 110309 110309.c
$ ./110309
要素0の入力<<5
要素1の入力<<3
要素2の入力<<8
要素3の入力<<6
要素4の入力<<4
要素0は5です
要素1は3です
要素2は8です
要素3は6です
要素4は4です
a
1 #include<stdio.h>
2
3 int main(void)
4 {
5 int hairetu[5];
6 int i;
7 for(i = 0; i < 5; i++)
8 {
9 printf("要素%dの入力<<", i);
10 scanf("%d", &hairetu[i]);
11 }
12 for(i = 0; i < 5; i++)
13 {
14 printf("要素%dは%dです\n", i, hairetu[i]);
15 }
16 return 0;
17 }
Esc
:wq
$ gcc -o 110309 110309.c
$ ./110309
要素0の入力<<5
要素1の入力<<3
要素2の入力<<8
要素3の入力<<6
要素4の入力<<4
要素0は5です
要素1は3です
要素2は8です
要素3は6です
要素4は4です
繰り返し文を作ったり、入力した数値の積を返したり、配列を導入してみたり....
$ vim 110307.1.c
a
1 #include<stdio.h>
2
3 int main(void)
4 {
5 int i;
6
7 for(i = 1; i <= 5; i++)
8 {
9 printf("%d回目の繰り返しです\n", i);
10 }
11 return 0;
12 }
13
Esc
:wq
gcc -o 110307.1 110307.1.c
$ ./110307.1
1回目の繰り返しです
2回目の繰り返しです
3回目の繰り返しです
4回目の繰り返しです
5回目の繰り返しです
$ vim 110307.2.c
a
1 #include <stdio.h>
2
3 int main(void)
4 {
5 int suti1, suti2;
6 int kekka;
7 printf("整数を2つ入力してください>>");
8 scanf("%d %d", &suti1, &suti2);
9
10 kekka = suti1 * suti2;
11
12 printf("%dと%dの積はは%dです\n",suti1,suti2,kekka);
13 return 0;
14 }
Esc
:wq
$ gcc -o 110307.2 110307.2.c
$ ./110307.2
整数を2つ入力してください>>3 9
3と9の積はは27です
$ vim 110307.3.c
a
1 #include <stdio.h>
2
3 int main(void)
4 {
5 char mozi1;
6 char mozi2;
7
8 printf("一文字目を入力してください");
9 scanf("%c", &mozi1);
10
11 printf("二文字目を入力してください");
12 scanf("%c", &mozi2);
13
14
15 printf("入力された値は%cと%cです\n", mozi1, mozi2);
16
17 return 0;
18 }
Esc
:wq
$ ./110307.3
一文字目を入力してくださいa
二文字目を入力してください入力された値はaと
です
$ ./110307.3
一文字目を入力してくださいab
二文字目を入力してください入力された値はaとbです
$ ./110307.3
一文字目を入力してくださいa b
二文字目を入力してください入力された値はaと です
$ vim 110307.4.c
a
1 #include<stdio.h>
2
3 int main(void)
4 {
5 int suti1,suti2,suti3,suti4,suti5;
6
7 suti1=1;
8 suti2=2;
9 suti3=3;
10 suti4=4;
11 suti5=5;
12
13 printf("数値1は%dです\n", suti1);
14 printf("数値2は%dです\n", suti2);
15 printf("数値3は%dです\n", suti3);
16 printf("数値4は%dです\n", suti4);
17 printf("数値5は%dです\n", suti5);
18
19 return 0;
20 }
Esc
:wq
$ gcc -o 110307.4 110307.4.c
$ ./110307.4
数値1は1です
数値2は2です
数値3は3です
数値4は4です
数値5は5です
$ vim 110307.5.c
a
1 #include <stdio.h>
2
3 int main(void)
4 {
5 int hairetu[5];
6 int i;
7
8 hairetu[0] = 1;
9 hairetu[1] = 2;
10 hairetu[2] = 3;
11 hairetu[3] = 4;
12 hairetu[4] = 5;
13
14 for(i = 0; i <5; i++)
15 {
16 printf("要素%dは%dです\n", i, hairetu[i]);
17 }
18 return 0;
19 }
Esc
:wq
$ gcc -o 110307.5 110307.5.c
$ ./110307.5
要素0は1です
要素1は2です
要素2は3です
要素3は4です
要素4は5です
a
1 #include<stdio.h>
2
3 int main(void)
4 {
5 int i;
6
7 for(i = 1; i <= 5; i++)
8 {
9 printf("%d回目の繰り返しです\n", i);
10 }
11 return 0;
12 }
13
Esc
:wq
gcc -o 110307.1 110307.1.c
$ ./110307.1
1回目の繰り返しです
2回目の繰り返しです
3回目の繰り返しです
4回目の繰り返しです
5回目の繰り返しです
$ vim 110307.2.c
a
1 #include <stdio.h>
2
3 int main(void)
4 {
5 int suti1, suti2;
6 int kekka;
7 printf("整数を2つ入力してください>>");
8 scanf("%d %d", &suti1, &suti2);
9
10 kekka = suti1 * suti2;
11
12 printf("%dと%dの積はは%dです\n",suti1,suti2,kekka);
13 return 0;
14 }
Esc
:wq
$ gcc -o 110307.2 110307.2.c
$ ./110307.2
整数を2つ入力してください>>3 9
3と9の積はは27です
$ vim 110307.3.c
a
1 #include <stdio.h>
2
3 int main(void)
4 {
5 char mozi1;
6 char mozi2;
7
8 printf("一文字目を入力してください");
9 scanf("%c", &mozi1);
10
11 printf("二文字目を入力してください");
12 scanf("%c", &mozi2);
13
14
15 printf("入力された値は%cと%cです\n", mozi1, mozi2);
16
17 return 0;
18 }
Esc
:wq
$ ./110307.3
一文字目を入力してくださいa
二文字目を入力してください入力された値はaと
です
$ ./110307.3
一文字目を入力してくださいab
二文字目を入力してください入力された値はaとbです
$ ./110307.3
一文字目を入力してくださいa b
二文字目を入力してください入力された値はaと です
$ vim 110307.4.c
a
1 #include<stdio.h>
2
3 int main(void)
4 {
5 int suti1,suti2,suti3,suti4,suti5;
6
7 suti1=1;
8 suti2=2;
9 suti3=3;
10 suti4=4;
11 suti5=5;
12
13 printf("数値1は%dです\n", suti1);
14 printf("数値2は%dです\n", suti2);
15 printf("数値3は%dです\n", suti3);
16 printf("数値4は%dです\n", suti4);
17 printf("数値5は%dです\n", suti5);
18
19 return 0;
20 }
Esc
:wq
$ gcc -o 110307.4 110307.4.c
$ ./110307.4
数値1は1です
数値2は2です
数値3は3です
数値4は4です
数値5は5です
$ vim 110307.5.c
a
1 #include <stdio.h>
2
3 int main(void)
4 {
5 int hairetu[5];
6 int i;
7
8 hairetu[0] = 1;
9 hairetu[1] = 2;
10 hairetu[2] = 3;
11 hairetu[3] = 4;
12 hairetu[4] = 5;
13
14 for(i = 0; i <5; i++)
15 {
16 printf("要素%dは%dです\n", i, hairetu[i]);
17 }
18 return 0;
19 }
Esc
:wq
$ gcc -o 110307.5 110307.5.c
$ ./110307.5
要素0は1です
要素1は2です
要素2は3です
要素3は4です
要素4は5です
annotationファイルに基づいて発現変動遺伝子のプローブにアノテーションをつける方法のモデルの考察(2)
$ R
genename <- c("gene1","geen2","gene3","gene4","gene5")
genename
[1] "gene1" "geen2" "gene3" "gene4" "gene5"
expres <- c(1.1,1.2,1.9,1.7,1.8)
probename <- c("A","B","C","D","E")
probename
[1] "A" "B" "C" "D" "E"
exprelist <-cbind(genename, expres)
exprelist
genename expres
[1,] "gene1" "1.1"
[2,] "geen2" "1.2"
[3,] "gene3" "1.9"
[4,] "gene4" "1.7"
[5,] "gene5" "1.8"
rownames(genename) <- probename
rownames(exprelist) <- probename
exprelist #これがアノテーションのためのマトリックス
genename expres
A "gene1" "1.1"
B "geen2" "1.2"
C "gene3" "1.9"
D "gene4" "1.7"
E "gene5" "1.8"
difgenes <- c("A","C","E") #実験でこのような発現変動プローブが出たとする(A,C,Eのプローブが発現変動を示したとする)
difgenes
[1] "A" "C" "E"
exprelist[difgenes,] #プローブ名を指定してやって抜き出してやる!!
genename expres
A "gene1" "1.1"
C "gene3" "1.9"
E "gene5" "1.8"
よし!!!!!発現変動遺伝子リストに遺伝子名のアノテーションをつけられたぞ!!!
genename <- c("gene1","geen2","gene3","gene4","gene5")
genename
[1] "gene1" "geen2" "gene3" "gene4" "gene5"
expres <- c(1.1,1.2,1.9,1.7,1.8)
probename <- c("A","B","C","D","E")
probename
[1] "A" "B" "C" "D" "E"
exprelist <-cbind(genename, expres)
exprelist
genename expres
[1,] "gene1" "1.1"
[2,] "geen2" "1.2"
[3,] "gene3" "1.9"
[4,] "gene4" "1.7"
[5,] "gene5" "1.8"
rownames(genename) <- probename
rownames(exprelist) <- probename
exprelist #これがアノテーションのためのマトリックス
genename expres
A "gene1" "1.1"
B "geen2" "1.2"
C "gene3" "1.9"
D "gene4" "1.7"
E "gene5" "1.8"
difgenes <- c("A","C","E") #実験でこのような発現変動プローブが出たとする(A,C,Eのプローブが発現変動を示したとする)
difgenes
[1] "A" "C" "E"
exprelist[difgenes,] #プローブ名を指定してやって抜き出してやる!!
genename expres
A "gene1" "1.1"
C "gene3" "1.9"
E "gene5" "1.8"
よし!!!!!発現変動遺伝子リストに遺伝子名のアノテーションをつけられたぞ!!!
annotationファイルに基づいて発現変動遺伝子のプローブにアノテーションをつける方法のモデルの考察(1)
$ R
probename <- c("A","B","C","D","E")
probename
[1] "A" "B" "C" "D" "E"
genename <- c("gene1","geen2","gene3","gene4","gene5")
genename
[1] "gene1" "geen2" "gene3" "gene4" "gene5"
annotationlist <- cbind(probename, genename)
annotationlist
probename genename
[1,] "A" "gene1"
[2,] "B" "geen2"
[3,] "C" "gene3"
[4,] "D" "gene4"
[5,] "E" "gene5"
class(annotationlist)
[1] "matrix"
annot <- as.data.frame(annotationlist)
annot
probename genename
1 A gene1
2 B geen2
3 C gene3
4 D gene4
5 E gene5
class(annot)
[1] "data.frame"
annot$probename
[1] A B C D E
Levels: A B C D E
difgenes <- c("A", "E")
difgenes
[1] "A" "E"
annot[annot$probename == "A",]
probename genename
1 A gene1
probename <- c("A","B","C","D","E")
probename
[1] "A" "B" "C" "D" "E"
genename <- c("gene1","geen2","gene3","gene4","gene5")
genename
[1] "gene1" "geen2" "gene3" "gene4" "gene5"
annotationlist <- cbind(probename, genename)
annotationlist
probename genename
[1,] "A" "gene1"
[2,] "B" "geen2"
[3,] "C" "gene3"
[4,] "D" "gene4"
[5,] "E" "gene5"
class(annotationlist)
[1] "matrix"
annot <- as.data.frame(annotationlist)
annot
probename genename
1 A gene1
2 B geen2
3 C gene3
4 D gene4
5 E gene5
class(annot)
[1] "data.frame"
annot$probename
[1] A B C D E
Levels: A B C D E
difgenes <- c("A", "E")
difgenes
[1] "A" "E"
annot[annot$probename == "A",]
probename genename
1 A gene1
Webページをpdf化して保存するには・・・・・
Save as PDF
https://addons.mozilla.org/ja/firefox/addon/save-as-pdf/
があります(^-^)
https://addons.mozilla.org/ja/firefox/addon/save-as-pdf/
があります(^-^)
ソースコードとバイナリコードを味わう(。-_-。)ノ☆
ソースコード:
プログラム言語でかかれたプログラムそのもの。
人間にはわかるが、コンピュータにはわからない為コンピュータにわかる言語体系に変換(トランスレート)する必要がある。
トランスレータにはアセンブラ、コンパイラ、インタプリタといった種類がある。
バイナリコード:
一般にコンピュータが直接理解できる言語(=機械語)で記述されたプログラムをさす。
しかし、現在では言葉が多様化しておりバイナリコード=機械語プログラムとは限らない。
$ vim Hello.c
a
1 #include <stdio.h>
2 int main(void)
3 {
4 printf("こんにちは!\n");
5 return 0;
6 }
Esc
:wq
$ gcc -o Hello Hello.c
$ vim Hello
^?ELF^A^A^A^@^@^@^@^@^@^@^@^@^B^@^C^@^A^@^@^@0<83>^D^H4^@^@^@(^Q^@^@^@^@^@^@4^@ ^@^H^@(^@^^^@^[^@^F^@^@^@4^@^@^@4<80>^D^H4<80>^D^H^@^A^@^@^@^A^@^@^E^@^@^@^D^@
・・・以下省略・・・
$ ./Hello.c
bash: ./Hello.c: Permission denied
$ ./Hello
こんにちは!
$ /home/kappa/2011/C/Hello
こんにちは!
プログラム言語でかかれたプログラムそのもの。
人間にはわかるが、コンピュータにはわからない為コンピュータにわかる言語体系に変換(トランスレート)する必要がある。
トランスレータにはアセンブラ、コンパイラ、インタプリタといった種類がある。
バイナリコード:
一般にコンピュータが直接理解できる言語(=機械語)で記述されたプログラムをさす。
しかし、現在では言葉が多様化しておりバイナリコード=機械語プログラムとは限らない。
$ vim Hello.c
a
1 #include <stdio.h>
2 int main(void)
3 {
4 printf("こんにちは!\n");
5 return 0;
6 }
Esc
:wq
$ gcc -o Hello Hello.c
$ vim Hello
^?ELF^A^A^A^@^@^@^@^@^@^@^@^@^B^@^C^@^A^@^@^@0<83>^D^H4^@^@^@(^Q^@^@^@^@^@^@4^@ ^@^H^@(^@^^^@^[^@^F^@^@^@4^@^@^@4<80>^D^H4<80>^D^H^@^A^@^@^@^A^@^@^E^@^@^@^D^@
・・・以下省略・・・
$ ./Hello.c
bash: ./Hello.c: Permission denied
$ ./Hello
こんにちは!
$ /home/kappa/2011/C/Hello
こんにちは!
Cでswitch文
$ vim 110305.switch.c
#include <stdio.h>
2
3 int main(void)
4 {
5 int suti;
6
7 printf("整数を入力してください>>");
8 scanf("%d", &suti);
9
10 switch(suti)
11 {
12 case 1987:
13 printf("その数字は私の生まれた年号です\n");
14 break;
15 case 5:
16 printf("その数字は私の生まれた月です\n");
17 break;
18 case 22:
19 printf("その数字は私の生まれた日です\n");
20 break;
21 default:
22 printf("残念ですが、その数字について何も話すネタはありません\n");
23 break;
24 }
25 return 0;
26 }
Esc
:wq
$ ./110305.switch
整数を入力してください>>1987
その数字は私の生まれた年号です
$ ./110305.switch
整数を入力してください>>5
その数字は私の生まれた月です
$ ./110305.switch
整数を入力してください>>22
その数字は私の生まれた日です
$ ./110305.switch
整数を入力してください>>3
残念ですが、その数字について何も話すネタはありません
#include <stdio.h>
2
3 int main(void)
4 {
5 int suti;
6
7 printf("整数を入力してください>>");
8 scanf("%d", &suti);
9
10 switch(suti)
11 {
12 case 1987:
13 printf("その数字は私の生まれた年号です\n");
14 break;
15 case 5:
16 printf("その数字は私の生まれた月です\n");
17 break;
18 case 22:
19 printf("その数字は私の生まれた日です\n");
20 break;
21 default:
22 printf("残念ですが、その数字について何も話すネタはありません\n");
23 break;
24 }
25 return 0;
26 }
Esc
:wq
$ ./110305.switch
整数を入力してください>>1987
その数字は私の生まれた年号です
$ ./110305.switch
整数を入力してください>>5
その数字は私の生まれた月です
$ ./110305.switch
整数を入力してください>>22
その数字は私の生まれた日です
$ ./110305.switch
整数を入力してください>>3
残念ですが、その数字について何も話すネタはありません
Cで分岐条件
$ vim 110305.分岐処理.c
a
1 #include <stdio.h>
2 int main(void)
3 {
4 int suti;
5
6 printf("整数を入力してください>>");
7 scanf("%d", &suti);
8
9 if (suti >= 0)
10 {
11 printf("sutiは0以上の整数です\n");
12 }
13 else
14 {
15 printf("sutiは0未満の整数です\n");
16 }
17 return 0;
18 }
Esc
$ gcc -o 110305.分岐処理 110305.分岐処理.c
$ ./110305.分岐処理
$ ./110305.分岐処理
整数を入力してください>>8
sutiは0以上の整数です
$ ./110305.分岐処理
整数を入力してください>>8
sutiは0以上の整数です
$ ./110305.分岐処理
整数を入力してください>>0
sutiは0以上の整数です
$ ./110305.分岐処理
整数を入力してください>>-1
sutiは0未満の整数です
a
1 #include <stdio.h>
2 int main(void)
3 {
4 int suti;
5
6 printf("整数を入力してください>>");
7 scanf("%d", &suti);
8
9 if (suti >= 0)
10 {
11 printf("sutiは0以上の整数です\n");
12 }
13 else
14 {
15 printf("sutiは0未満の整数です\n");
16 }
17 return 0;
18 }
Esc
$ gcc -o 110305.分岐処理 110305.分岐処理.c
$ ./110305.分岐処理
$ ./110305.分岐処理
整数を入力してください>>8
sutiは0以上の整数です
$ ./110305.分岐処理
整数を入力してください>>8
sutiは0以上の整数です
$ ./110305.分岐処理
整数を入力してください>>0
sutiは0以上の整数です
$ ./110305.分岐処理
整数を入力してください>>-1
sutiは0未満の整数です
Cでprintfの書式指定だよ☆
$ vim printf.c
a
1 #include <stdio.h>
2
3 int main(void)
4 {
5 char mozi;
6 int seisu;
7 double zissu;
8
9
10 mozi = 'A';
11 printf("char mozi: %c\n", mozi);
12
13 seisu = 10000;
14 printf("int seisu: %d\n", seisu);
15
16 zissu = 12345.678;
17 printf("double zissu: %lf\n", zissu);
18
19 return 0;
20 }
Esc
:wq
$ gcc -o printf printf.c
$ ./printf
a
1 #include <stdio.h>
2
3 int main(void)
4 {
5 char mozi;
6 int seisu;
7 double zissu;
8
9
10 mozi = 'A';
11 printf("char mozi: %c\n", mozi);
12
13 seisu = 10000;
14 printf("int seisu: %d\n", seisu);
15
16 zissu = 12345.678;
17 printf("double zissu: %lf\n", zissu);
18
19 return 0;
20 }
Esc
:wq
$ gcc -o printf printf.c
$ ./printf
perlでpeptideを記入したファイルを読み込んでシェルに表示させるプログラムを書いた♪
$ vim 110305_test.pep
a
1 DDKLEGLFLKCGGIDEMQSS
Esc
;wq
$ vim 110305_test.perl
a
1 use strict;
2 use warnings;
3
4 #Reading protein sequence data from a file
5
6 #The filename of the file containing the protein sequence data
7 my $proteinfilename = '110305_test.pep';
8
9 #First we have to #open" the file, and associate
10 #a "filehandle" with it. We choose the filehandle
11 #PROTEINFILE for readability.
12 open(PROTEINFILE, $proteinfilename);
13
14 #Now we do the actual reading of the protein sequence data fromj the file,
15 #by using the angle brackets <and> to get teh input from the
16 #filehandle. We store the data into our variable $protein.
17 my $protein = <PROTEINFILE>;
18
19 #Now that we've got our data, we can close the file.
20 close PROTEINFILE;
21
22 #Print the protein onto the screen
23 print "Here is the protein:\n\n";
24
25 print $protein;
Esc
:wq
$ perl 110305_test perl
Here is the protein:
DDKLEGLFLKCGGIDEMQSS
うん。美しい♪♪
Pymolを用いて、PDBのデータを美しく表示してみせる♪♪
pymolaをapt-getを用いてインストール
$ sudo apt-get install pymol
さあ、PDBjへ行こう♪♪
GUIにて
http://www.pdbj.org/
を入力
COX-2 (PDB ID = 5cox)をPDBjよりダウンロードしてみる。
次は、pymolを用いてタンパク質の立体構造を描写してみます(^^)
pdb5cox.ent
をダウンロードして、カレントディレクトリにおいて
$ pymol pdb5cox.ent
H - Hide - everythin (一度すべての表示を消すっていう意味)
とした後
S - Show - cartoon
とすると下のような表示(cartoon)になる。すごくきれい(^^)
さらにcartoonの図にlinesを重ねることもできる。
このままの状態で
S - Show - lines
さらには
L - Label - resudues
とした後、もう一つのwindowの方で、
Display - zoom
ってすれば拡大した絵ができる♪♪
美しすぎるタンパク質\(○^ω^○)/
$ sudo apt-get install pymol
さあ、PDBjへ行こう♪♪
GUIにて
http://www.pdbj.org/
を入力
説明を追加 |
COX-2 (PDB ID = 5cox)をPDBjよりダウンロードしてみる。
次は、pymolを用いてタンパク質の立体構造を描写してみます(^^)
pdb5cox.ent
をダウンロードして、カレントディレクトリにおいて
$ pymol pdb5cox.ent
H - Hide - everythin (一度すべての表示を消すっていう意味)
とした後
S - Show - cartoon
とすると下のような表示(cartoon)になる。すごくきれい(^^)
さらにcartoonの図にlinesを重ねることもできる。
このままの状態で
S - Show - lines
さらには
L - Label - resudues
とした後、もう一つのwindowの方で、
Display - zoom
ってすれば拡大した絵ができる♪♪
美しすぎるタンパク質\(○^ω^○)/
IntelのLinux用Fortranをインストールを試みた♪
IntelのLinux用Fortranをインストールを試みた♪
http://software.intel.com/en-us/articles/non-commercial-software-development/に行きます。
IntelのLinux用Fortranをインストールを試みた♪
Intelからメールがくるはず!!
メールには、サンプルFortranのシリアルナンバーとアプリケーションのインストーラーのダウンロード先が書いてある。それにしたがって以下のライセンスとパッケージをダウンロードする。
NCOM___NR2M-P257VJBJ.lic (ライセンスは適当なディレクトリにおいておけばいい)
l_fcompxe_2011.2.137.tgz
これを適当なディレクトリに展開するとこんな感じになる。
kappa@kappa-desktop:~/2011/ダウンロード関連/l_fcompxe_2011.2.137$ ls
install.sh pset support.txt cd_eject.sh license rpms
ここで、
$ ./install.sh
でインストーラーをを起動させる。
あとは、指示にしたがってEnterキーを押し続ければよい。
途中でシリアルナンバーを入力しなければならないので、注意が必要である。
インストールが完了したら、パスを指定してあげる必要がある。
一番簡単なのは、
$ source /opt/intel/bin/ifortvars.sh ia32
と入力してあげれば、Fortranは使用可能になる。
例えば、
$ ifort -help
・・・・激しく文字が表示される・・・・
Intelからメールがくるはず!!
メールには、サンプルFortranのシリアルナンバーとアプリケーションのインストーラーのダウンロード先が書いてある。それにしたがって以下のライセンスとパッケージをダウンロードする。
NCOM___NR2M-P257VJBJ.lic (ライセンスは適当なディレクトリにおいておけばいい)
l_fcompxe_2011.2.137.tgz
これを適当なディレクトリに展開するとこんな感じになる。
kappa@kappa-desktop:~/2011/ダウンロード関連/l_fcompxe_2011.2.137$ ls
install.sh pset support.txt cd_eject.sh license rpms
ここで、
$ ./install.sh
でインストーラーをを起動させる。
あとは、指示にしたがってEnterキーを押し続ければよい。
途中でシリアルナンバーを入力しなければならないので、注意が必要である。
インストールが完了したら、パスを指定してあげる必要がある。
一番簡単なのは、
$ source /opt/intel/bin/ifortvars.sh ia32
と入力してあげれば、Fortranは使用可能になる。
例えば、
$ ifort -help
・・・・激しく文字が表示される・・・・
終了(^^)
perlで与えられたDNAをひっくり返すプログラムの作成♪♪
$ vim 110302.perl
a
1 use strict;
2 use warnings;
3
4 #Calculating the reverse complement of a strand od DNA
5
6 # The DNA
7 my $DNA = 'ACGGGAGGACGGGGAAAAATTACTACGGCATTAGC';
8
9 #Print the DNA onto the screen
10 print "Here is the starting DNA:\n\n";
11
12 print "$DNA\n\n";
13
14 #Calculate the reverse complement
15 #Warning : this attempt will fail!
16 #
17 #First, copy the DNA into new variable $revcom
18
19 my $revcom = reverse $DNA;
20
21 #Nexxt substitute all bases by their complements,
22 #A->T, T->A, G->C, C->G
23
24 $revcom =~ s/A/T/g;
25 $revcom =~ s/T/A/g;
26 $revcom =~ s/G/C/g;
27 $revcom =~ s/C/G/g;
28
29
30 #Print the reverse complement DNA onto the screen
31 print "Here is the reverse complement DNA:\n\n";
32
33 print "$revcom\n";
34
35 print "\nThat was a bad algorithm, and the reverse complement was wrong!\n";
36 print "Try again...\n\n";
37
38 #Make a new copy of the DNA
39 my $revcom = reverse $DNA;
40 $revcom =~ tr/ACGTacgt/TGCAtgca/;
41
42 #Print the reverse complement DNA onto the screen
43 print "Here is the reverse complement DNA:\n\n";
44
45 print "$revcom\n";
46 print"\nThis time ist warked!!\n\n";
Esc
$ perl 110302.perl
Here is the starting DNA:
ACGGGAGGACGGGGAAAAATTACTACGGCATTAGC
Here is the reverse complement DNA:
GGAAAAGGGGAAGAAAAAAAAGGGGGAGGAGGGGA
That was a bad algorithm, and the reverse complement was wrong!
Try again...
Here is the reverse complement DNA:
GCTAATGCCGTAGTAATTTTTCCCCGTCCTCCCGT
This time ist warked!!
素敵過ぎる♪♪♪
a
1 use strict;
2 use warnings;
3
4 #Calculating the reverse complement of a strand od DNA
5
6 # The DNA
7 my $DNA = 'ACGGGAGGACGGGGAAAAATTACTACGGCATTAGC';
8
9 #Print the DNA onto the screen
10 print "Here is the starting DNA:\n\n";
11
12 print "$DNA\n\n";
13
14 #Calculate the reverse complement
15 #Warning : this attempt will fail!
16 #
17 #First, copy the DNA into new variable $revcom
18
19 my $revcom = reverse $DNA;
20
21 #Nexxt substitute all bases by their complements,
22 #A->T, T->A, G->C, C->G
23
24 $revcom =~ s/A/T/g;
25 $revcom =~ s/T/A/g;
26 $revcom =~ s/G/C/g;
27 $revcom =~ s/C/G/g;
28
29
30 #Print the reverse complement DNA onto the screen
31 print "Here is the reverse complement DNA:\n\n";
32
33 print "$revcom\n";
34
35 print "\nThat was a bad algorithm, and the reverse complement was wrong!\n";
36 print "Try again...\n\n";
37
38 #Make a new copy of the DNA
39 my $revcom = reverse $DNA;
40 $revcom =~ tr/ACGTacgt/TGCAtgca/;
41
42 #Print the reverse complement DNA onto the screen
43 print "Here is the reverse complement DNA:\n\n";
44
45 print "$revcom\n";
46 print"\nThis time ist warked!!\n\n";
Esc
$ perl 110302.perl
Here is the starting DNA:
ACGGGAGGACGGGGAAAAATTACTACGGCATTAGC
Here is the reverse complement DNA:
GGAAAAGGGGAAGAAAAAAAAGGGGGAGGAGGGGA
That was a bad algorithm, and the reverse complement was wrong!
Try again...
Here is the reverse complement DNA:
GCTAATGCCGTAGTAATTTTTCCCCGTCCTCCCGT
This time ist warked!!
素敵過ぎる♪♪♪
Perlで与えたDNAの塩基配列をRNAの塩基配列に変更する。
$ vim 110301.pl
a
1 use strict;
2 use warnings;
3
4 #The DNA
5 my $DNA = 'ACGGGAGGACGGGAAAATTACTACGGCATTAGC';
6
7 # Print the DNA onto the screen
8 print "Here is the starting DNA : \n\n";
9
10 print"$DNA\n\n";
11
12 #Transcribe the DNA to RNA by substituting all T's with U's.
13 my $RNA = $DNA;
14
15 $RNA =~ s/T/U/g;
16
17 #Print the RNA onto the screen
18 print "Here is the result of transcribing the DNA to RNA : \n\n";
19
20 print"$RNA\n";
21
Esc
$ perl 110301.pl
Here is the starting DNA :
ACGGGAGGACGGGAAAATTACTACGGCATTAGC
Here is the result of transcribing the DNA to RNA :
ACGGGAGGACGGGAAAAUUACUACGGCAUUAGC
a
1 use strict;
2 use warnings;
3
4 #The DNA
5 my $DNA = 'ACGGGAGGACGGGAAAATTACTACGGCATTAGC';
6
7 # Print the DNA onto the screen
8 print "Here is the starting DNA : \n\n";
9
10 print"$DNA\n\n";
11
12 #Transcribe the DNA to RNA by substituting all T's with U's.
13 my $RNA = $DNA;
14
15 $RNA =~ s/T/U/g;
16
17 #Print the RNA onto the screen
18 print "Here is the result of transcribing the DNA to RNA : \n\n";
19
20 print"$RNA\n";
21
Esc
$ perl 110301.pl
Here is the starting DNA :
ACGGGAGGACGGGAAAATTACTACGGCATTAGC
Here is the result of transcribing the DNA to RNA :
ACGGGAGGACGGGAAAAUUACUACGGCAUUAGC
Cでキーボードから二つの変数を入力して表示させる♪
$ vim 110301.c
a
1 # include<stdio.h>
2
3 int main(void)
4 {
5 int suti1;
6 int suti2;
7
8 printf("整数を二つ入力してください。>>");
9 scanf("%d %d", &suti1, &suti2);
10
11 printf("入力された値は%dと%dです!!\n", suti1, suti2);
12
13 return 0;
14 }
Esc
$ gcc -o 110301 110301.c
$ ./ 110301
整数を二つ入力してください。>>45
59
入力された値は45と59です!!
$ vim 110301.2.c
a
1 #include<stdio.h>
2
3 int main(void)
4 {
5 int suti1;
6 int suti2;
7
8 printf("一つ目の整数を入力してください。>>");
9 scanf("%d",&suti1);
10
11 printf("二つ目の整数を入力してください。>>");
12 scanf("%d",&suti2);
13
14 printf("入力された値は%dと%dです。\n", suti1, suti2);
15 return 0;
16 }
Esc
$ gcc -o 110301.2 110301.2.c
$ ./110301.2
一つ目の整数を入力してください。>>45
二つ目の整数を入力してください。>>59
入力された値は45と59です。
a
1 # include<stdio.h>
2
3 int main(void)
4 {
5 int suti1;
6 int suti2;
7
8 printf("整数を二つ入力してください。>>");
9 scanf("%d %d", &suti1, &suti2);
10
11 printf("入力された値は%dと%dです!!\n", suti1, suti2);
12
13 return 0;
14 }
Esc
$ gcc -o 110301 110301.c
$ ./ 110301
整数を二つ入力してください。>>45
59
入力された値は45と59です!!
$ vim 110301.2.c
a
1 #include<stdio.h>
2
3 int main(void)
4 {
5 int suti1;
6 int suti2;
7
8 printf("一つ目の整数を入力してください。>>");
9 scanf("%d",&suti1);
10
11 printf("二つ目の整数を入力してください。>>");
12 scanf("%d",&suti2);
13
14 printf("入力された値は%dと%dです。\n", suti1, suti2);
15 return 0;
16 }
Esc
$ gcc -o 110301.2 110301.2.c
$ ./110301.2
一つ目の整数を入力してください。>>45
二つ目の整数を入力してください。>>59
入力された値は45と59です。