医科学のためのデータシミュレーション、のための予備知識 補遺

  • こちらの続き
  • 医学・生物学データの扱いに慣れるための、データシミュレーションセミナーの資料(医学研究科大学院コース ゲノム遺伝学@京大 2014:こちら)
  • htmlファイルはこちら(Login as a guest でログインすると見られます)
  • こちらを参考にすればhtml化可能です
# 関数を作る Make functions

## 関数と入出力 Functions and their inputs and outputs

関数は仕事をします。
Functions do jobs.

関数は
* 名前を持っていて
* 入力を受け取り
* 出力を返します

Functions
* have their name
* take inputs
* returns outputs.

```{r}
# "my1stFx" is "name of function" you make
# x is input(s)
# y is output(s)
my1stFx <- function(x){
  y <- 2*x
  return(y)
}
# 使ってみます。Usage
x <- 1:3
x
my1stFx(x)
# 出力は格納できます You can put the output into an object
my1stOut <- my1stFx(x)
my1stOut
```

入出力は複数にもできる。
Inputs and outputs can be multiple.

```{r}
my2ndFx <- function(x,a){
  y <- x+1
  b <- a + rnorm(length(a))
  return(list(xAdd1=y,aAddz=b))
}
x <- 1:10
a <- x*2+3
my2ndFx(x,a)
my2ndOut <- my2ndFx(x,a)
my2ndOut
my2ndOut$xAdd1
my2ndOut$aAddz
plot(x,a)
plot(my2ndOut$xAdd1,my2ndOut$aAddz,col=2)
```

## if条件分岐とループ If conditions and Looping

### if
```{r}
# 多すぎるデータの場合にだけ一部を取り出す
# When number of records is too many, sample the elements.
my3rdFx <- function(x,n){
  if(length(x)>n){
    y <- sample(x,n)
  }else{
    y <- x
  }
  return(y)
}

x <- rnorm(50)
n <- 100
my3rdOut <- my3rdFx(x,n)
length(my3rdOut)

x.2 <- rnorm(110)
my3rdOut.2 <- my3rdFx(x.2,n)
length(my3rdOut.2)
```

### ループ Looping

A,T,G,Cの塩基配列をn.iter個作る。
長さは20塩基前後とし、A,T,G,Cは等確率でランダムに選ぶとする。

Let's make n.iter base-sequences with A,T,G,and C.
Their length should be around 20 and bases should be chosen with eaqual probability at random.
```{r}
n.iter <- 5
mean.len <- 20
Seqs <- list()
ATGC <- c("A","T","G","C")
pr <- rep(1/4,4)
for(i in 1:n.iter){
  len <- rpois(1,mean.len)
  Seqs[[i]] <- sample(ATGC,len,replace=TRUE,prob=pr)
}
Seqs
```

ランダムな塩基配列を数を増やして作成し、どんな長さ分布かを見てみることとする。
Let's generate much more sequences and check the distribution of length of generated sequences.

```{r}
n.iter <- 1000
mean.len <- 20
Seqs <- list()
ATGC <- c("A","T","G","C")
pr <- rep(1/4,4)
for(i in 1:n.iter){
  len <- rpois(1,mean.len)
  Seqs[[i]] <- sample(ATGC,len,replace=TRUE,prob=pr)
}

len.check <- lapply(Seqs,length)
hist(unlist(len.check))
len.check.s <- sapply(Seqs,length)
hist(len.check.s)
mean(unlist(len.check))
mean(len.check.s)
```

### if とループの組み合わせ Combine if and looping

上記でやったように適当な配列を作りつつ、読み取りエラーを想定に加える。
うまく読めなかった塩基を"N"として加え、"N"の入った配列は捨て、その上で配列を指定の本数だけ作成する。

As we did above, we generate random sequences but this time, bases might not be read well. "N" represents the unread bases and sequences that has one or more "N" should be discarded. We want to have n.iter sequences without "N".

まずは"N"で排除しない場合。
"N" is allowed at first.
```{r}
n.iter <- 10
mean.len <- 20
Seqs <- list()
ATGCN <- c("A","T","G","C","N")
pr <- c(rep(1/4,4),0.05)
pr <- pr/sum(pr)
pr
for(i in 1:n.iter){
  len <- rpois(1,mean.len)
  Seqs[[i]] <- sample(ATGCN,len,replace=TRUE,prob=pr)
}
Seqs
```

"N"があったら登録しない。
When "N" is included, a vacant sequence is registered.

```{r}
n.iter <- 10
mean.len <- 20
Seqs <- list()
ATGCN <- c("A","T","G","C","N")
pr <- c(rep(1/4,4),0.05)
pr <- pr/sum(pr)
pr
for(i in 1:n.iter){
  len <- rpois(1,mean.len)
  tmp.seq <- sample(ATGCN,len,replace=TRUE,prob=pr)
  checker <- which(tmp.seq=="N")
  if(length(checker)>0){
    Seqs[[i]] <- ""
  }else{
    Seqs[[i]] <- tmp.seq
  }
}
Seqs
```

"N"が含まれているときには登録自体を飛ばすことにする。
When "N" is included, no registration.

```{r}
n.iter <- 10
mean.len <- 20
Seqs <- list()
ATGCN <- c("A","T","G","C","N")
pr <- c(rep(1/4,4),0.05)
pr <- pr/sum(pr)
pr
# カウンタを導入する
cnt <- 1
for(i in 1:n.iter){
  len <- rpois(1,mean.len)
  tmp.seq <- sample(ATGCN,len,replace=TRUE,prob=pr)
  checker <- which(tmp.seq=="N")
  if(length(checker)>0){
    #Seqs[[i]] <- ""
  }else{
    Seqs[[cnt]] <- tmp.seq
    cnt <- cnt + 1
  }
}
Seqs
```

指定の数の配列を作る。
Let's generate given-number of sequences.

```{r}
n.iter <- 10
mean.len <- 20
Seqs <- list()
ATGCN <- c("A","T","G","C","N")
pr <- c(rep(1/4,4),0.05)
pr <- pr/sum(pr)
pr
# カウンタを導入する
cnt <- 1
# ループを続けるかどうかの判断のため
# Continue/discontinue looping?
continueLoop <- TRUE
while(continueLoop){
  len <- rpois(1,mean.len)
  tmp.seq <- sample(ATGCN,len,replace=TRUE,prob=pr)
  checker <- which(tmp.seq=="N")
  if(length(checker)>0){
    #Seqs[[i]] <- ""
  }else{
    Seqs[[cnt]] <- tmp.seq
    cnt <- cnt + 1
  }
  # 登録本数が十分になったらループをやめる
  # When adequate sequences are registered, stop looping
  if(cnt==n.iter+1){
    continueLoop <- FALSE
  }
}
Seqs
```