---
title: "[Web連携テキスト バイオインフォマティクス](index.html)"
date: '最終更新: `r format(Sys.Date(), "%Y/%m/%d")`'
output: 
  html_document:
    toc: TRUE
    toc_float: TRUE
---
```{r klippy, echo=FALSE, include=TRUE}
klippy::klippy()
```

# 付録

## R関連情報(第1章)  
<p>本書と直接関係のない部分もありますが、初心者の方は基本的にはざっとでよいので順番通りにやっていくとよいと思います。</p>

### R1.010 R本体とRStudioとパッケージのインストール  
- [Windows版(2023/05/31更新)](R/R1.010_win.pdf)  
- [Macintosh版(2023/06/14更新)](R/R1.010_mac.pdf)  

### R1.020 RStudioの基本的な利用法  
- [Windows版(2023/04/06更新)](R/R1.020_win.pdf)  
- [Macintosh版(2023/04/06更新)](R/R1.020_mac.pdf)  

### R1.030 1.2.1項の相補鎖など(基礎)
Rに慣れることを目的として、「[1.2.1 サンガー法の原理](#121_サンガー法の原理)」の冒頭部分で述べた相補鎖や逆相補鎖をRで得るやり方を解説します。まずは[Biostrings](http://bioconductor.org/packages/release/bioc/html/Biostrings.html)パッケージをロードして、このパッケージが提供する`complement`や`reverseComplement`関数を使えるようにします。
```{r, message=FALSE, error=TRUE}
library(Biostrings)
```
<br>
次に、本文中で示したものと同じ塩基配列を`sequenc`という名前のオブジェクトに格納し、その中身を表示させます。背景が灰色になっている部分がコマンド、##から始まっている背景が白くなっている部分が実行結果です。当たり前といえば当たり前ですが、入力と同じ配列が出力されていることがわかります。
```{r, message=FALSE, error=TRUE}
sequenc <- "AGTGACGGTCTT"
sequenc
```
<br>
今取り扱っている`sequenc`オブジェクトの中身は、見た目上明らかに塩基配列です。しかし、[Biostrings](http://bioconductor.org/packages/release/bioc/html/Biostrings.html)パッケージが提供する相補鎖や逆相補鎖を得るための関数を利用するためには、これが塩基配列だと正しく認識させておかねばなりません。そこで次に、`sequenc`を入力として`DNAStringSet`という関数を実行し、その結果を`fasta`という名前のオブジェクトに格納する作業を行います。先ほどと同じく、背景が灰色になっている部分がコマンド、##から始まっている背景が白くなっている部分が実行結果です。背景が白の実行結果部分において、`width`列の12という数値が塩基数、`seq`列が実際の塩基配列になっていることがわかります。左下に見えている`[1]`はまだ気にしなくてもよいです。最初の行で見えている`DNAStringSet object`というのが、(塩基)配列情報を[Biostrings](http://bioconductor.org/packages/release/bioc/html/Biostrings.html)パッケージで取り扱うために定義された特有の形式なのだという理解でよいです。
```{r, message=FALSE, error=TRUE}
fasta <- DNAStringSet(sequenc)
fasta
```
<br>
いよいよ本番1として、相補鎖を得るための`complement`関数を実行します。入力は`DNAStringSet`形式の`fasta`オブジェクトです。確かに実行結果として見えているものは、本文中のものと同じ"TCACTGCCAGAA"です。
```{r, message=FALSE, error=TRUE}
complement(fasta)
```
<br>
本番2として、逆相補鎖を得るための`reverseComplement`関数を実行します。確かに実行結果として見えているものは、本文中のものと同じ"AAGACCGTCACT"です。なお、以下では`reverseComplement(fasta)`としていますが、`reverseComplement(x=fasta)`と書いても構いません。
```{r, message=FALSE, error=TRUE}
reverseComplement(fasta)
```

### R1.040 エラーメッセージに慣れる
次に、存在しないオブジェクト`uge`の中身を見ようとして、どういうエラーが出るかを眺めます。このWebページは[R Markdown](https://rmarkdown.rstudio.com/)で自動的に作成しているため若干エラーメッセージの内容が異なっているかもしれませんが、「`オブジェクト 'uge' がありません(object 'uge' not found`)」のようなメッセージ表示されていることがわかります。このように、**原因が明らかなエラーメッセージを眺めておけば、実際に遭遇する問題に対処しやすい**です。  
```{r, message=FALSE, error=TRUE}
uge
```
<br>
以下では、存在しない`uge`を入力として、逆相補鎖を返す`reverseComplement`関数の実行を試みています。重要なことは、エラーメッセージの中から**対処すべき事柄が「`オブジェクト 'uge' がありません(object 'uge' not found`)」だけ**だと読み解くことです。  
```{r, message=FALSE, error=TRUE}
reverseComplement(uge)
```
<br>
逆相補鎖を返す`reverseComplement`関数は、入力として`DNAStringSet object`を受け付けます。そこで次に、ただの文字列である`sequenc`オブジェクトを入力として実行してみます。エラーメッセージだけを見ると、何をどうすればよいのかさっぱりわかりませんが、原因既知状態で眺めると穏やかに原因が特定できると思います。たとえば「`..."character"...`」が目につくと思いますが、`sequenc <- "AGTGACGGTCTT"`からも明らかなように、`sequenc`の中身はcharacterです。
```{r, message=FALSE, error=TRUE}
sequenc <- "AGTGACGGTCTT"
reverseComplement(sequenc)
```
<br>
任意のオブジェクトを入力として`class`関数を実行すると、そのオブジェクトがどういうものかを調べることができます。`class`実行結果も、エラーへの対応(これを**デバッグ**といいます)を行う上で有用な情報となります。
```{r, message=FALSE, error=TRUE}
class(sequenc)
fasta <- DNAStringSet(sequenc)
class(fasta)
```
<br>
次に、存在しない`uge`オブジェクトを入力として、再度`reverseComplement(uge)`の実行結果を眺めます。よく見ると、自分が指定した覚えのない「引数 'x' の評価中にエラーが起きました」のようなエラーメッセージが出ていることがわかります。これは、この関数の場合は、入力である`DNAStringSet object`を受け付ける引数名がxだからです。それがわかっていれば、「`reverseComplement(uge)`で入力として与えた`uge`の引数名がxだな。なのでこれは引数名を省略しているが、実際には`reverseComplement(x=uge)`だからこのようなメッセージが出ているのだろう」と解釈できるようになります。
```{r, message=FALSE, error=TRUE}
reverseComplement(uge)
```
<br>
相補鎖や逆相補鎖を得る程度なら、入力として与える情報は1つだけです。しかし複雑な関数になってくると関数内で複数のオプション(引数)を与えて実行します。複数オプションがある場合は、その関数が想定するデフォルトの順番通りだと引数名を明示しなくても構いません。しかし、入力する側が理解しやすいように引数名を明示することもあります。このあたりのポリシーはヒトそれぞれですので、ユーザ側がこの段階で慣れておくとよいと思います。`?reverseComplement`は相補鎖を返す`reverseComplement`関数のマニュアルを表示するコマンドです。`help(reverseComplement)`でも同じ挙動になります。
```{r, eval=FALSE}
?reverseComplement
```
若干見栄えは異なりますが、上記実行結果とほぼ同じ情報を以下に示します。①や②で見えている`x`が、これらの関数で入力として与える情報の引数名です。[RStudio](https://www.rstudio.com/)上で上記コマンドを実行すると、RStudio画面右下のHelpタブで小さく表示されて見づらい場合は、たとえば[Biostrings](http://bioconductor.org/packages/release/bioc/html/Biostrings.html)パッケージのページからリンクが張られている[Reference Manual](https://bioconductor.org/packages/release/bioc/manuals/Biostrings/man/Biostrings.pdf)で、調べたい関数名でページ内検索する手もあります。以下の図は「`reverseComplement`」でページ内検索してとったスクショです。
![?reverseComplement実行結果](R/R1.040.png){width="600"}

### R1.050 1.2.1項の相補鎖など(応用)
[RStudio](https://www.rstudio.com/)上で塩基配列を手入力し、実行結果が[RStudio](https://www.rstudio.com/)上に表示されるだけだと実感がわかない方もいらっしゃると思います。ここでは、塩基配列(`AGTGACGGTCTT`)をを[FASTA形式](https://en.wikipedia.org/wiki/FASTA_format)ファイルとして独立に用意したもの([R1.050.fasta](R/R1.050.fasta))を入力として用いるやり方を示します。以下は、入力ファイルとして[R1.050.fasta](R/R1.050.fasta)を指定して、[Biostrings](http://bioconductor.org/packages/release/bioc/html/Biostrings.html)パッケージをロードして、このパッケージが提供する`readDNAStringSet`関数を用いて読み込んだ結果を`fasta`という名前のオブジェクトに格納し、表示させるところまでのスクリプトです。背景が灰色になっている部分がコマンド、##から始まっている背景が白くなっている部分が実行結果です。背景が白の実行結果部分において、`width`列の12という数値が塩基数、`seq`列が実際の塩基配列、そして`names`列の情報がFASTA形式ファイルのdescription行部分に相当することがわかります。  
```{r, message=FALSE, error=TRUE}
in_f <- "R1.050.fasta"
library(Biostrings)
fasta <- readDNAStringSet(in_f, format="fasta")
fasta
```
<br>
次に、逆相補鎖を得たのち、それを`out_f`で指定した出力ファイルに保存するまでのスクリプトを示します。`reverseComplement(fasta)`実行結果を`out`オブジェクトに格納し、それを`writeXStringSet`関数を用いて保存していることが読み解けます。`writeXStringSet`関数実行時に、最初の引数として`out`を指定しているのは納得できると思います(逆相補鎖情報を含む`out`オブジェクトの中身を保存したいからです)。次に出力ファイル名の引数に相当する`file`という引数で、出力ファイル名情報が格納されている`out_f`を与えています。出力形式は、この場合は`"fasta"`形式の一択になりますが念のため明示しています。`width`という引数では、1行あたり何塩基とするかを指定するオプションになります。この場合は12塩基しかないので特に気にする必要はありませんが、染色体など長い配列になると、`width=60`などを指定して1行あたりの塩基数を多めに設定するのを好むヒトもいます。
```{r, message=FALSE, error=TRUE}
in_f <- "R1.050.fasta"
out_f <- "R1.050_out.fasta"
library(Biostrings)
fasta <- readDNAStringSet(in_f, format="fasta")
out <- reverseComplement(fasta)
writeXStringSet(out, file=out_f, format="fasta", width=50)
```
<br>
念のため、以下に逆相補鎖情報を含む`out`オブジェクトの中身を示しておきます。エラーなく実行できたヒトは、出力ファイル名として指定した[R1.050_out.fasta](R/R1.050_out.fasta)が作業ディレクトリ上に保存されているはずです。
```{r, message=FALSE, error=TRUE}
out
```

### R1.060 1.2.2項のクオリティスコア
(1.1)式は、ある塩基のベースコールがエラーとなる確率pを入力として、その塩基のクオリティスコア*q*を返すものです。本文中でいくつかの入力と出力の関係をいくつか例示していますが、簡単ですので[RStudio](https://www.rstudio.com/products/rstudio/)での計算例を示します。以下は、p = 0.01の結果です。確かに出力結果としてq = 20が得られていることがわかります。
```{r, message=FALSE, error=TRUE}
p <- 0.01
q <- -10*log(p, base=10)
q
```
<br>
次は、数値ベクトルを入力として与えるやり方を示します。たとえば、p = 0.0001のときのqは40だと読み解きます。
```{r, message=FALSE, error=TRUE}
p <- c(0.1, 0.01, 0.001, 0.0001)
q <- -10*log(p, base=10)
q
```
<br>
この場合のqは、4つの要素からなる数値ベクトルです。ベクトルの要素数を調べる`length`関数の入力として与えて実行すると、確かに4という結果が返されます。
```{r, message=FALSE, error=TRUE}
length(q)
```
<br>
ベクトルの2番目の要素情報を取り出したい場合は`q[2]`、4番目の要素情報を取り出したい場合は`q[4]`のように指定します。
```{r, message=FALSE, error=TRUE}
q[2]
q[4]
```
<br>
このベクトルqの5番目の要素は存在しないことが分かっていますが、もし`q[5]`とやると、Not Availableを意味する`NA`という結果が返されます。
```{r, message=FALSE, error=TRUE}
q[5]
```

### R1.070 1.7.1項の*k*-mer解析
<p>1.7.1項の終盤から1.7.2項の最初のほうで述べている*k*-mer解析に関する基礎的な事柄に加えて、ショートリードデータを用いたゲノムサイズ推定までの内容を参考資料として以下に提供します。特に、ゲノムサイズ推定は、1.8.2項でも述べたものです。元ネタは、2022年5月17日に行った計3時間分の講義資料です。講義で利用されたい方は、pptx版もありますので、必要に応じてご自由に改変してご利用ください。</p>
  - [*k*-mer解析の基礎と応用(pptx版)](R/R1.070.pptx)
  - [*k*-mer解析の基礎と応用(pdf版)](R/R1.070.pdf)

### R1.080 1.8.4項の例題1.1(基礎)
1.8.4項の最後にある[例題1.1](chapter1/E1.1.pdf)をRでやることもできます。この項目の主目的は、ちょっとややこしいですがそれができるのを実感することです。`readDNAStringSet`, `strsplit`, `nchar`, `rep`, `c`などのRの関数や、Rのノリに慣れてもらうことも兼ねています。ここでは、[図1.8](chapter1/fig1.8.pdf)aの塩基配列を[FASTA形式](https://en.wikipedia.org/wiki/FASTA_format)ファイルとして独立に用意したもの([R1.080.fasta](R/R1.080.fasta))を入力として用いるやり方を示します。以下は、入力ファイルとして[R1.080.fasta](R/R1.080.fasta)を指定して、[Biostrings](http://bioconductor.org/packages/release/bioc/html/Biostrings.html)パッケージをロードして、このパッケージが提供する`readDNAStringSet`関数を用いて読み込んだ結果を`fasta`という名前のオブジェクトに格納するところまでのスクリプトです。  
```{r, message=FALSE, error=TRUE}
in_f <- "R1.080.fasta"
library(Biostrings)
fasta <- readDNAStringSet(in_f, format="fasta")
```
<br>
以下では、入力ファイルの中身に相当する`fasta`オブジェクトを表示させています。背景が灰色になっている部分がコマンド、##から始まっている背景が白くなっている部分が実行結果です。`DNAStringSet`という独特の形式でFASTA形式ファイルの中身が格納されているのだと理解すればよいです。`width`列の60という数値が塩基数、`seq`列が実際の塩基配列の最初と最後の部分、そして`names`列の情報が[FASTA形式](https://en.wikipedia.org/wiki/FASTA_format)のdescription行部分に相当します。随分冗長な表現方法だと思われるかもしれませんが、たとえば入力が数十個の染色体配列から構成されるヒトゲノム配列ファイルを読み込んだ場合をイメージしていただければ納得しやすいかと思います。  
```{r, message=FALSE, error=TRUE}
fasta
```
<br>
以下では、仮想制限酵素の認識配列がACなので、それを用いて`fasta`オブジェクト内の塩基配列(文字列なのでstring)を分割(split)する`strsplit`関数を実行しています。実行結果を`hoge`オブジェクトに格納しています。2行目の`strsplit`関数実行部分の右側に`[[1]]`がついていますが、簡単にいうと話をややこしくしないようにするためのおまじないです。なお、`strsplit`関数実行部分で「`Warning in .normarg_at2(at, x): 'at' names were ignored`」という警告メッセージが表示されていますが、エラーではないので気にする必要はありません。  
**補足：**例題1.1をRでやるのは、あくまでもおまけ程度の位置づけです。例題番号的にこの位置にしましたが、結果的にちょっとややこしい内容のものが含まれています。それゆえ、初心者の方は「Rだとこんな感じでできる」ことを実感する程度の位置づけでよいと思います。もう少し真面目に解説すると、`strsplit`関数の出力がリスト形式とよばれるもので、そこから情報を取り出す際に大括弧(角括弧)が2つ並んだ`[[1]]`のような感じで情報を取り出しています。この場合は配列が1つしかないので、`[[1]]`として1番目の要素を取り出しています。  
```{r, message=FALSE, error=TRUE}
param <- "AC"
hoge <- strsplit(fasta, param)[[1]]
```
<br>
以下では、`hoge`の中身を表示しています。`hoge`の中身と[例題1.1](chapter1/E1.1.pdf)の解答例を見比べればわかりますが、ACで分割後の配列数(=5)が一致していることがわかります。また、2番目以降の配列の冒頭にACが付加されていないことがわかりますが、これは`strsplit`関数の性質上やむを得ません。もしかしたら制限酵素認識配列専用の関数を誰かが作っているのかもしれませんが、実用上は自分が知っている関数の中から近い機能を果たすものをうまく使い、残りは*ad hoc*(その場しのぎ、のようなイメージでよいです)に対応するという感じでよいと思います。
```{r, message=FALSE, error=TRUE}
hoge
```
<br>
*ad hoc*な対応の例として、たとえば`width(hoge)`とやることで、`hoge`オブジェクト中の計5配列の長さ情報を数値ベクトルとして得ることができます。  
```{r, message=FALSE, error=TRUE}
width(hoge)
```
<br>
また、今は仮想制限酵素が2塩基だとわかっていますので、以下のように`width(hoge) + 2`とすることで、ベクトル中の全要素に2を足すことができます。
```{r, message=FALSE, error=TRUE}
width(hoge) + 2
```
<br>
上記は、`param`オブジェクトの中身が2文字から構成されることが分かっているのでそのまま2としました。もちろんもう少し一般性を高めた表現も可能です。例えば、オブジェクト中の文字数をカウントする`nchar`関数を用いて、`nchar(param)`として利用することもできます。こんな感じで少しずつRの関数に慣れていくとよいでしょう。
```{r, message=FALSE, error=TRUE}
nchar(param)
width(hoge) + nchar(param)
```
<br>
ただ上記の結果だと、最初の要素にまで2が足されてしまっています。それゆえ、何を足すかという情報をあらかじめ作成しておいてから、`width(hoge)`の対応する要素と足し合わせるようなことをすればよいです。この場合は、1番目が0、2-5番目の要素が2の計5つの要素からなる数値ベクトルを作成して足すようなイメージです。`rep`関数は、1番目の引数として指定した情報を、2番目の引数で指定した個数だけ複製(replicate)します。以下では、`rep(2, 4)`として`2`を`4`つ作成しています。引数名を明示すると、`rep(x=2, times=4)`と書くこともできます。
```{r, message=FALSE, error=TRUE}
rep(2, 4)
```
<br>
以下は、`rep(nchar(param), 4)`として、実質的に中身が2の`nchar(param)`にしているだけです。どちらでもいいじゃないかと思われるかもしれませんが、こんな感じで一般式で書くことによって手作業で変更せねばならないところを極力減らして汎用性を高めていきます。
```{r, message=FALSE, error=TRUE}
rep(nchar(param), 4)
```
<br>
以下では、ベクトルを連結(combine)する`c`関数を用いて、最初の要素が0、2番目以降が2という数値ベクトルを連結させて得たものを`tasuyou`というオブジェクトに格納しています。ここまでで、`width(hoge)`の対応する要素に足すべき数値が格納された`tasuyou`オブジェクトの完成です。Rの教え方はヒトそれぞれですが、**このあたりを見るだけでも基本テクの習得が重要**であることがわかると思います。  
```{r, message=FALSE, error=TRUE}
tasuyou <- c(0, rep(nchar(param), 4))
tasuyou
```
<br>
残りは、`width(hoge)`と`tasuyou`の対応する要素を以下のように足せばよいです。ここまでで[例題1.1](chapter1/E1.1.pdf)の解答例の最後に示した光学マップ情報が得られたことになります。
```{r, message=FALSE, error=TRUE}
width(hoge) + tasuyou
```

### R1.090 1.8.4項の例題1.1(応用)
上記項目の続きですが、[RStudio](https://www.rstudio.com/products/rstudio/)起動直後の状態でも、以下を実行すれば必要な情報を取得できます。この項目の主目的は、`replaceAt`関数の紹介や、`as.character`関数実行結果オブジェクトの取り扱いに関する小技の伝授です。なお、`strsplit`関数実行部分で「`Warning in .normarg_at2(at, x): 'at' names were ignored`」という警告メッセージが表示されていますが、エラーではないので気にする必要はありません。
```{r, message=FALSE, error=TRUE}
in_f <- "R1.080.fasta"
param <- "AC"
library(Biostrings)
fasta <- readDNAStringSet(in_f, format="fasta")
hoge <- strsplit(fasta, param)[[1]]
```
<br>
`AC`で分割された配列からなる`hoge`オブジェクト中の各配列に対して、`1`番目の塩基の位置に`param`で指定した塩基(つまり`AC`)を挿入することもできます。これも本当は2番目以降の配列のみに挿入したいのですが、ややこしくなるので`replaceAt`という関数を用いればこういうこともできるんだという程度の情報にとどめておきます。以下では、`param`, `hoge`, `replaceAt(hoge, 1, value=param)`の実行結果を一気に表示しています。背景が灰色になっている部分がコマンド、##から始まっている背景が白くなっている部分が実行結果です。
```{r, message=FALSE, error=TRUE}
param
hoge
replaceAt(hoge, 1, value=param)
```
<br>
`strsplit`関数の入力は、入力が塩基配列であることを明示した`DNAStringSet`形式の`fasta`オブジェクトではなく、ただの文字列(as character)としたものを与えることもできます。この場合は、`as.character(fasta)`を入力とすればよいです。
```{r, message=FALSE, error=TRUE}
fasta
as.character(fasta)
```
<br>
以下を実行して得られる、出力結果の`hoge2`を眺めると、`$fig1.8a`というのが最初のほうに見えていて気持ち悪いと思われるかもしれません。慣れてくると、こういうの見た段階で`hoge2$fig1.8a`のような感じで取り扱えばよいと閃くことができるようになります。
```{r, message=FALSE, error=TRUE}
hoge2 <- strsplit(as.character(fasta), param)
hoge2
```
<br>
たとえば、`hoge2`で見えている計5配列の配列長を調べたい場合は、`hoge2$fig1.8a`を入力として、文字数をカウントする`nchar`関数を実行すればよいです。
```{r, message=FALSE, error=TRUE}
nchar(hoge2$fig1.8a)
```
また、見慣れた`DNAStringSet`形式にしたい場合は、`hoge2$fig1.8a`を入力として、`DNAStringSet`関数を実行すればよいです。
```{r, message=FALSE, error=TRUE}
DNAStringSet(hoge2$fig1.8a)
```

### R1.100 1.10.2項のGFF3(ファイル取得)
[1.10.2項](https://www.iu.a.u-tokyo.ac.jp/textbook/#1102_%E3%82%A2%E3%83%8E%E3%83%86%E3%83%BC%E3%82%B7%E3%83%A7%E3%83%B3%E3%81%AE%E5%85%A5%E5%87%BA%E5%8A%9B%E3%81%A8%E5%8F%82%E7%85%A7DB%E3%81%AE%E5%BD%B9%E5%89%B2)で述べた、アノテーション結果を格納する代表的な形式である[GFF3](https://en.wikipedia.org/wiki/General_feature_format)ファイルを[RStudio](https://www.rstudio.com/products/rstudio/)で読み込んで取り扱うやり方を示します。ここではまず、[Ensembl](https://www.ensembl.org/) ([Cunningham et al., Nucleic Acids Res., 2022](http://www.ncbi.nlm.nih.gov/pubmed/34791404))のバクテリア版である[EnsemblBacteria (Release 54)](https://bacteria.ensembl.org/)上で、乳酸菌の1つである[Lactobacillus casei 21/1](https://bacteria.ensembl.org/Lactobacillus_casei_21_1_gca_000309585/Info/Index/)のページを訪れて、全体像を把握します。以下はそのスクショです。まずは、このゲノムの基本情報を把握すべく、①または②のどちらでもよいのでクリックします。  
![?Lactobacillus casei 21/1のトップページ1](R/R1.100.1.png){width="600"}  
<br>
③より、このゲノムサイズは3,215,878 bp (約3.2 MB)だということがわかります。また④より、coding genesが3,080個、non-coding genesが68個、pseudogenesが88個、転写物全体としては3,080 + 68 + 88 = 3,236個だということがわかります。  
![?Lactobacillus casei 21/1の基本情報1](R/R1.100.2.png){width="600"}  
<br>
次に、一旦前のページに戻ります。GFF3ファイルは⑤から取得可能です。ゲノム配列ファイルは⑥から取得できますが、⑦のFASTA形式ファイルの中からも取得できます。なお、ゲノム配列ファイルについては、[1.10.3項](https://www.iu.a.u-tokyo.ac.jp/textbook/#1103_%E3%82%B2%E3%83%8E%E3%83%A0%E4%B8%AD%E3%81%AE%E5%8F%8D%E5%BE%A9%E9%85%8D%E5%88%97)で述べる反復配列をマスクしないバージョン('dna')とマスクしたバージョン(マスクの程度によってさらに'dna_rm'と'dna_sm'があります)など様々なものが用意されています。ここではマスクしないバージョンのものを取り扱います。なお、ファイルはすべて[gzip](https://ja.wikipedia.org/wiki/Gzip)という形式(拡張子は.gz)で圧縮されていますので、圧縮・解凍ソフトで解凍せねばなりません(本当はFASTA形式ファイルは圧縮ファイルでも`readDNAStringSet`関数で読み込めますが、gff3のほうがダメなので全部解凍しといたほうがややこしくないと思います)。  
![?Lactobacillus casei 21/1のトップページ2](R/R1.100.3.png){width="600"}  
<br>
<p>上記の情報は[EnsemblBacteria (Release 54)](https://bacteria.ensembl.org/)のものですが、リリースの違いによってリンク切れになったりします。それゆえ、以降の解析が数年後でも可能となるよう、念のため以下でも解凍後のファイルへのリンクを貼っておきます：</p>
  - [Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.49.gff3](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.49.gff3)  
  - [Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.cdna.all.fa](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.cdna.all.fa)  
全部で3,168配列あります。3,080 coding genes + 88 pseudogenesと一致しています。  
  - [Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.cds.all.fa](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.cds.all.fa)  
全部で3,080配列あります。3,080 coding genesと一致しています。  
  - [Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa)  
全部で75配列あります。総塩基数は3,215,878 bpですので、確かにゲノムサイズと一致しています。  
  - [Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.ncrna.fa](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.ncrna.fa)  
68配列あります。non-coding genesと一致しています。  
  - [Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.pep.all.fa](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.pep.all.fa)  
全部で3,080配列あります。3,080 coding genesと一致しています。ただしこれはアミノ酸配列なので、読み込み時に用いる関数は`readDNAStringSet`ではなく`readAAStringSet`にせねばなりません。  

### R1.110 1.10.2項のGFF3(読込)
[1.10.2項](https://www.iu.a.u-tokyo.ac.jp/textbook/#1102_%E3%82%A2%E3%83%8E%E3%83%86%E3%83%BC%E3%82%B7%E3%83%A7%E3%83%B3%E3%81%AE%E5%85%A5%E5%87%BA%E5%8A%9B%E3%81%A8%E5%8F%82%E7%85%A7DB%E3%81%AE%E5%BD%B9%E5%89%B2)で述べた、アノテーション結果を格納する代表的な形式である[GFF3](https://en.wikipedia.org/wiki/General_feature_format)ファイルを[RStudio](https://www.rstudio.com/products/rstudio/)で読み込むやり方を示します。[rtracklayer](https://bioconductor.org/packages/release/bioc/html/rtracklayer.html)パッケージ([Lawrence et al., Bioinformatics, 2009](https://pubmed.ncbi.nlm.nih.gov/19468054/))が提供する`readGFF`関数や、[GenomicsFeatures](https://bioconductor.org/packages/release/bioc/html/GenomicFeatures.html)([Lawrence et al., PLoS Comput Biol., 2013](https://pubmed.ncbi.nlm.nih.gov/23950696/))が提供する`makeTxDbFromGFF`関数がありますが、ここでは後者のやり方を示します。[GenomicsFeatures](https://bioconductor.org/packages/release/bioc/html/GenomicFeatures.html)を未インストールの方は、まず以下を実行してください。  

```{r, message=FALSE, error=TRUE, warning=FALSE, eval=FALSE}
BiocManager::install("GenomicFeatures", update=F)
```
<br>
次に、[GenomicsFeatures](https://bioconductor.org/packages/release/bioc/html/GenomicFeatures.html)([Lawrence et al., PLoS Comput Biol., 2013](https://pubmed.ncbi.nlm.nih.gov/23950696/))が提供する`makeTxDbFromGFF`関数を用いて、[Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.49.gff3](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.49.gff3)を入力として読み込むやり方を示します。背景が灰色になっている計4行から構成される以下のコマンドでは、まず入力ファイルをin_fというオブジェクト名で取り扱うという宣言をしています。次に、[GenomicsFeatures](https://bioconductor.org/packages/release/bioc/html/GenomicFeatures.html)パッケージを`library`関数を用いてロードしています。これでGenomicFeaturesが提供する関数を利用可能になったので、(入力ファイル名に相当する)`in_f`を入力として`makeTxDbFromGFF`関数を実行した結果を`txdb`というオブジェクトに格納し、その中身を表示しています。  
```{r, message=FALSE, error=TRUE}
in_f <- "Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.49.gff3"
library(GenomicFeatures)
txdb <- makeTxDbFromGFF(in_f, format="auto")
txdb
```
<br>
背景が白の実行結果の冒頭部分に`TxDb object`と書かれていることと、`makeTxDbFromGFF`関数を実行した結果であることから想像できるように、`makeTxDbFromGFF`関数はGFF形式ファイルを読み込んでTxDbという[GenomicsFeatures](https://bioconductor.org/packages/release/bioc/html/GenomicFeatures.html)パッケージ内で有効な独特の形式でGFF情報を格納しているのだと理解すればよいです。重要なことは、無事に読み込めることと、txdbの中身に相当する背景が白の実行結果から見られる情報として「`Nb of transcripts: 3236`」に注目することです。この3236という数値は、確かにダウンロード元である[Lactobacillus casei 21/1](https://bacteria.ensembl.org/Lactobacillus_casei_21_1_gca_000309585/Info/Index/)から辿れる[More information and statistics](https://bacteria.ensembl.org/Lactobacillus_casei_21_1_gca_000309585/Info/Annotation/#assembly)で見られるGene transcriptsの数と一致しています。GFFファイル自体は1つしかないため、すべての情報を含んでないと都合が悪いのだと考えれば納得できると思います。念のために以下にもリンク先のスクショを示しておきます。①の部分の数値と一致しているということです。  
![?Lactobacillus casei 21/1の基本情報2](R/R1.110.1.png){width="600"}  

### R1.120 1.10.2項のGFF3(座標取得)
`TxDb object`である`txdb`オブジェクトだけを見せられても実感がわかないと思います。`txdb`の中に含まれているのは、ゲノム配列中のどこにどのようなフィーチャ(アノテーションがつけられるゲノム上の各領域のこと)があるのかという座標情報です。たとえば`txdb`を入力として`transcripts`関数を実行すると、`GRanges object`という形式の、この場合は`transcripts`なので転写物の座標情報を得ることができます。  
`seqnames`列は、アノテーションに用いた大元のゲノム配列ファイル([Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa))内のdescription行で見られる配列名情報に相当します。1番目の転写物は`ctg01_00128`という配列上、3236番目の転写物は`ctg75_00141`という配列上に存在するものだと解釈します。`ctg`は、コンティグ(contig)を意味します。最もよくつながっているのが染色体レベル、その前段階がスカッフォールドレベル、そしてコンティグレベルとなります。この乳酸菌ゲノムは全部で75配列から構成されますが、バクテリア程度のゲノムサイズで75配列に分断されている程度のアセンブリ結果しか得られていないので、ある程度業界歴が長くなると「確かにコンティグレベルだろうな」というのがわかってきます。  `ranges`列はstartとendの位置情報、`strand`列は+鎖か-鎖のどちらかの情報、`tx_id`列は転写物のシリアル番号、`tx_name`列は転写物名を表しています。  
```{r, message=FALSE, error=TRUE}
in_f <- "Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.49.gff3"
library(GenomicFeatures)
txdb <- makeTxDbFromGFF(in_f, format="auto")
transcripts(txdb)
```
<br>
ノリさえわかってくれば、`txdb`を入力として`promoters`関数を実行することでプロモータ配列を取得することもできます。以下は、転写開始点(TSS)の上流`50`塩基、下流`10`塩基の領域情報を取得するコマンドです。たとえば1番目の転写物のTSSは、`ctg01_00128`という配列の`4769`番目の塩基です(+鎖なので、start位置がTSSです)。この転写物のプロモータ領域は`4719-4778`となっていますが、確かにTSSである4769 - 50 = 4719、そして4769 + 10 - 1 = 4778で妥当な結果です。もう1つの例として、3236番目の転写物物は-鎖上にあるので、転写開始点(TSS)は`ctg75_00141`という配列の`639`番目の塩基になります。それゆえ、上流`50`塩基は639 + 50 = 689番目の塩基、そして下流`10`塩基は639 - 10 + 1 = 630番目の塩基になります。したがって、取得したプロモータ領域は妥当です。  
```{r, message=FALSE, error=TRUE}
in_f <- "Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.49.gff3"
library(GenomicFeatures)
txdb <- makeTxDbFromGFF(in_f, format="auto")
promoters(txdb, upstream=50, downstream=10)
```

### R1.130 1.10.2項のGFF3(配列取得1)
ここでは、実際に取得した座標情報である`GRanges object`を利用して、その領域の塩基配列を取得します。ただし、GFF3ファイルの情報はただの座標情報にすぎないため、大元のGFF3ファイルと対応したゲノム配列も読み込んで利用せねばならない点に注意が必要です。ここではアノテーションに用いた大元のゲノム配列ファイル([Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa))も同時に読み込んでin_f1というオブジェクト名で利用しています。[Rsamtools](https://bioconductor.org/packages/release/bioc/html/Rsamtools.html)パッケージが提供する`getSeq`や`FaFile`関数も利用します。[Rsamtools](https://bioconductor.org/packages/release/bioc/html/Rsamtools.html)を未インストールの方は、まず以下を実行してください。  
```{r, message=FALSE, error=TRUE, warning=FALSE, eval=FALSE}
BiocManager::install("Rsamtools", update=F)
```
<br>
以下は、アノテーションのGFF3形式ファイルとゲノム配列のFASTA形式ファイルの2つを入力として与えて、最終的に`DNAStringSeq`形式の`fasta`オブジェクトを得るところまでのコードです。計3つのパッケージを明示的にロードしていることがわかります。実際問題としては依存関係的に明示する必要がないものもあるかもしれませんが、(ロードしておいて文句を言われることとパッケージをロードしていないことに起因するエラーを天秤にかけて)念のために`library`関数を用いてロードしています。なお、`library`の代わりに`require`関数を使われる方もいらっしゃいますが、実用上はどちらでも構いません。`makeTxDbFromGFF`関数でGFF3ファイル名に相当するin_f2を読み込んだ結果を`txdb`オブジェクトに格納し、`transcripts`関数で転写物の座標情報を`GRanges`形式の`hoge`オブジェクトに格納しています。最後に「ゲノム配列ファイルに相当する`in_f1`」と「先ほどのGFF3形式ファイルの中身に相当する`hoge`オブジェクト」の2つを入力として`getSeq`関数を実行することで、目的の`DNAStringSeq`形式の`fasta`オブジェクトを得ることができます。  
##から始まっている背景が白くなっている実行結果部分を見ると、確かに`fasta`オブジェクトは`DNAStringSeq`形式になっていることがわかります。`width`列は配列長情報で、たとえば1番最初の転写物の配列長は1086です。この転写物の`GRanges`形式のstartとendの位置は`4769-5854`ですので、5,854 - 4769 + 1 = 1,086なので妥当だと判断できます。ここでは表示していませんが、気になるかたは`transcripts(txdb)`実行結果に相当する`hoge`オブジェクトの中身を表示させて確認してみるとよいと思います。  
```{r, message=FALSE, error=TRUE}
in_f1 <- "Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa"
in_f2 <- "Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.49.gff3"
library(Rsamtools)
library(GenomicFeatures)
library(Biostrings)
txdb <- makeTxDbFromGFF(in_f2, format="auto")
hoge <- transcripts(txdb)
fasta <- getSeq(FaFile(in_f1), hoge)
fasta
```
<br>
以下は、得られた転写物配列情報をFASTA形式ファイルとして保存するところまで含めたスクリプトです。上のスクリプトとの違いは、3行目の出力ファイル名情報を格納するための`out_f`オブジェクト作成部分と、最後から2番目の行の`writeXStringSet`関数実行部分です。無事実行完了すると、出力ファイル名として指定したファイル([R1.130.1.fasta](R/R1.130.1.fasta))が作成されているはずです。ここまでで配列取得の基本形はが完成です。なお、最終行の`fasta[1:4]`は、「fastaオブジェクト中の最初の4配列分を表示せよ」という意味です。これは`head(fasta, n=4)`というコマンドと同じです。`head`は最初のn個分(デフォルトは6)を出力する関数であるため、`n=4`として最初の4個分とオプションで指定していることに相当します。  
```{r, message=FALSE, error=TRUE}
in_f1 <- "Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa"
in_f2 <- "Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.49.gff3"
out_f <- "R1.130.1.fasta"
library(Rsamtools)
library(GenomicFeatures)
library(Biostrings)
txdb <- makeTxDbFromGFF(in_f2, format="auto")
hoge <- transcripts(txdb)
fasta <- getSeq(FaFile(in_f1), hoge)
writeXStringSet(fasta, file=out_f, format="fasta", width=50)
fasta[1:4]
```
<br>
さきほどのファイル([R1.130.1.fasta](R/R1.130.1.fasta))のdescription行に相当する情報は、上で見えている`fasta`オブジェクトの`names`列です。"ctg01_00128"は、計75配列からなるゲノム配列ファイル([Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa))の1番最初の配列名です。この情報は、(計3,236配列分あるので画面がざっと流れますが)`names(fasta)`で一望できます。元のゲノム配列は計75配列(75 contigs)から構成されていますので、種類としては計75種類ということになります。どのコンティグ由来配列が何個あったかは、計3,236要素からなる文字列ベクトルである`names(fasta)`を入力として、以下のように`table`関数を実行すればよいです。実行結果は、たとえば"ctg01_00128"というコンティグ由来転写物は38個、"ctg02_00001"由来の転写物は17個という風に解釈します。  
```{r, message=FALSE, error=TRUE}
table(names(fasta))
```
<br>

### R1.140 1.10.2項のGFF3(配列取得2)
配列取得1では、転写物配列取得の基本形を示しました。以下は、先ほどの配列取得1の出力ファイル([R1.130.1.fasta](R/R1.130.1.fasta))を得るスクリプトの再掲です。`hoge`は`transcripts(txdb)`実行結果のGRanges形式のオブジェクト、`fasta`は`hoge`の情報を用いて転写物配列を取得したDNAStringSet形式のオブジェクトです。  
```{r, message=FALSE, error=TRUE}
in_f1 <- "Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa"
in_f2 <- "Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.49.gff3"
out_f <- "R1.130.1.fasta"
library(Rsamtools)
library(GenomicFeatures)
library(Biostrings)
txdb <- makeTxDbFromGFF(in_f2, format="auto")
hoge <- transcripts(txdb)
fasta <- getSeq(FaFile(in_f1), hoge)
writeXStringSet(fasta, file=out_f, format="fasta", width=50)
```
<br>
以下は、`hoge`オブジェクトの最初の3行分を表示させた結果です。最初の2列(`seqnames`と`ranges`列)は、その列名を関数として情報抽出に利用できます。  
```{r, message=FALSE, error=TRUE}
hoge[1:3]
```
<br>
以下は、`fasta`オブジェクトの最初の3行分を表示させた結果です。`names`列の情報に相当する`names(fasta)`は、`hoge`オブジェクトの`seqnames`列の情報と同じであることがわかります。  
```{r, message=FALSE, error=TRUE}
fasta[1:3]
```
<br>
以下は、`seqnames(hoge)`の実行結果です。見慣れないと面食らってしまうかもしれませんが、先程の`table(names(fasta))`実行結果と見比べれば、「ctg01_00128が38個ある」や「ctg02_00001が17個ある」といったことがわかります。ここでは示しませんが、`table(names(fasta))`の実行結果は、`table(seqnames(hoge))`の実行結果と同じです。また、`names(fasta)`と同じ見栄えの情報は、`as.character(seqnames(hoge))`で得られます。さらに、`table(seqnames(hoge))`と同じ結果は、`table(as.character(seqnames(hoge)))`でも得られます。いきなり`table(as.character(seqnames(hoge)))`を見ると難解だと思うかもしれませんが、このように様々な情報と照らし合わせて、要素ごとに分解して眺めると理解しやすいのではと思います。  
```{r, message=FALSE, error=TRUE}
seqnames(hoge)
```
<br>
次は、`ranges(hoge)`の実行結果です。確かに`hoge`オブジェクト中の`ranges`列で見えている情報が含まれています。`start`列と`end`列に加えて`width`列が見えていますが、これは(end - start +1)で算出される配列長情報です。  
```{r, message=FALSE, error=TRUE}
ranges(hoge)
```
<br>
`ranges(hoge)`の実行結果で見えている情報は、その列名を関数として当該列情報を抽出できます。具体的には、`start`列の情報は`start(ranges(hoge))`、`end`列の情報は`end(ranges(hoge))`、そして`width`列の情報は`width(ranges(hoge))`で取得可能です。以下では、(これをそのまま実行してしまうと、それぞれについて計3,236配列分の要素の情報が一気に表示されて見づらいため)さらに`head`関数をつけて最初の6個分に限定して示します。  
```{r, message=FALSE, error=TRUE}
head(start(ranges(hoge)))
head(end(ranges(hoge)))
head(width(ranges(hoge)))
```
<br>
ここで、再度転写物配列情報である`fasta`オブジェクトの最初の5行分を表示させます。`names`列の情報に相当する`names(fasta)`は[R1.130.1.fasta](R/R1.130.1.fasta)のdescription行に相当しますが、現状では由来配列名しか示されていないことがわかります。
```{r, message=FALSE, error=TRUE}
fasta[1:5]
```
<br>
そこで、これまで示してきたテクニック+αで、`names(fasta)`の情報を「*由来配列名*」から「*由来配列名*_*start*_*end*」に変更し、転写物配列がどの配列のどの領域由来かがdescription情報をみればわかるようにします。まず、新たな関数として`paste`を紹介します。`paste`は、複数の数値なり文字なりを、任意の文字や記号で連結する関数です。以下に4つ例を示します。最後のコマンドが実際にやろうとしていることに近いものになります。  
```{r, message=FALSE, error=TRUE}
paste(8, 5, sep="_")
paste("an", 5, "pan", sep="__")
paste("an", "pon", "tan", sep="")
paste(head(start(ranges(hoge))), head(end(ranges(hoge))), sep="_")
```
<br>
ここから本番になっていきますが、まずは変更前の`names(fasta)`の最初の6個分の情報を以下に示します。  
```{r, message=FALSE, error=TRUE}
head(names(fasta))
```
<br>
以下が、`names(fasta)`の情報を「*由来配列名*」から「*由来配列名*_*start*_*end*」に変更する本番のコマンドです。あまり横長になっても見づらいので2つの行にわけていますが、`paste`関数を使って、当該転写物の*由来配列名*に相当する情報の`seqnames(hoge)`、当該転写物の*start*位置に相当する情報の`start(ranges(hoge))`、そして当該転写物の*end*位置に相当する情報の`end(ranges(hoge))`を"\_"で連結した結果を、`names(fasta)`に格納(代入)しています。ここでは右辺の情報を統一的に`hoge`オブジェクトに揃える目的で*由来配列名*に相当する情報として`seqnames(hoge)`を用いましたが、`names(fasta)`で置き換えてもかまいません。その場合は、右辺の`names(fasta)`が「*由来配列名*」のみに相当し、左辺の`names(fasta)`が変更後の「*由来配列名*_*start*_*end*」ということになります。  
```{r, message=FALSE, error=TRUE}
names(fasta) <- paste(seqnames(hoge), start(ranges(hoge)),
                      end(ranges(hoge)), sep="_")
```
<br>
変更後の`names(fasta)`の最初の6個分の情報を以下に示すべく、以下で再度`head(names(fasta))`実行結果を示します。確かに意図通りにstartとendの情報が右側に追加されていることがわかります。  
```{r, message=FALSE, error=TRUE}
head(names(fasta))
```
<br>
注意点としては、`fasta`オブジェクトの中身がこの段階で変更されているということです。つまりこの場合は、`names`列の情報が変更されています。右端が...になっているものがありますが、上のほうで見えている結果を突き合わせることで意図通りに変更されていると確信できると思います。  
```{r, message=FALSE, error=TRUE}
fasta
```
<br>
最後に、`names(fasta)`の情報を更新したものをファイル([R1.130.2.fasta](R/R1.130.2.fasta))に保存するところまでの一連のスクリプトを示します。最終行の`fasta[1:3]`の実行結果を見れば、確かに出力ファイルのdescription行が正しく変更されていると確信できると思います。  
```{r, message=FALSE, error=TRUE}
in_f1 <- "Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa"
in_f2 <- "Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.49.gff3"
out_f <- "R1.130.2.fasta"
library(Rsamtools)
library(GenomicFeatures)
library(Biostrings)
txdb <- makeTxDbFromGFF(in_f2, format="auto")
hoge <- transcripts(txdb)
fasta <- getSeq(FaFile(in_f1), hoge)
names(fasta) <- paste(seqnames(hoge), start(ranges(hoge)),
                      end(ranges(hoge)), sep="_")
writeXStringSet(fasta, file=out_f, format="fasta", width=50)
fasta[1:3]
```

### R1.150 1.10.3項のリピートマスク
ここでは、リピートマスク前後の乳酸菌ゲノム配列を眺めます。[Ensembl](https://www.ensembl.org/) ([Cunningham et al., Nucleic Acids Res., 2022](http://www.ncbi.nlm.nih.gov/pubmed/34791404))のバクテリア版である[EnsemblBacteria (Release 54)](https://bacteria.ensembl.org/)上で、乳酸菌の1つである[Lactobacillus casei 21/1](https://bacteria.ensembl.org/Lactobacillus_casei_21_1_gca_000309585/Info/Index/)のページを再訪します。①の[Download DNA sequence](http://ftp.ensemblgenomes.org/pub/bacteria/release-54/fasta/bacteria_79_collection/lactobacillus_casei_21_1_gca_000309585/dna/)をクリックして、ゲノム配列のダウンロードページに移動します。  
![?Lactobacillus casei 21/1のゲノム配列のダウンロードページ1](R/R1.150.1.png){width="600"}  
<br>
以下がリンク先のスクショです。一番上の[CHECKSUM](http://ftp.ensemblgenomes.org/pub/bacteria/release-54/fasta/bacteria_79_collection/lactobacillus_casei_21_1_gca_000309585/dna/CHECKSUMS)は、ダウンロードしたファイルが壊れていないか(提供元と同一か)を調べるための情報です。一番下の[README](http://ftp.ensemblgenomes.org/pub/bacteria/release-54/fasta/bacteria_79_collection/lactobacillus_casei_21_1_gca_000309585/dna/README)は、計6種類見えているゲノム配列ファイルについての解説が書かれています。  
![?Lactobacillus casei 21/1のゲノム配列のダウンロードページ2](R/R1.150.2.png){width="600"}  
<br>
<p>上記の情報は[EnsemblBacteria (Release 54)](https://bacteria.ensembl.org/)のものですが、リリースの違いによってリンク切れになったりします。それゆえ、以降の解析が数年後でも可能となるよう、念のため以下でもファイルへのリンクを貼っておきます。"Lactobacillus_casei_21_1_gca_000309585.ASM30958v1"までが共通で、それ以降の文字列が異なっていることがわかります。たとえば1-2番目のファイルはは"dna"、3-4番目のファイルは"dna_rm"、そして5-6番目のファイルは"dna_sm"となっていることがわかります。"dna"はマスクされていない配列(unmasked)のことです。"dna_rm"は、[RepeatMasker](https://www.repeatmasker.org/)で同定された反復配列をマスクした(masked)配列のことで、具体的には"N"で置換した配列のことです。"dna_sm"は、[RepeatMasker](https://www.repeatmasker.org/)で同定された反復配列をソフトにマスクした(soft-masked)配列のことです。その次のnonchromosomalは、当該ゲノムの配列決定レベル(解像度)が最高の染色体レベルになったとしても、どの染色体にも割り当てられない断片配列が存在しますので、それらの配列も含めたものです。toplevelは、当該ゲノムの最高の配列決定レベルで得られた配列から構成されるファイルです。この場合は、配列決定精度がコンティグレベル(contig; ctg)ですので、nonchromosomalとtoplevelの中身は同じになります。直感的には、最も"N"が少ないのが"dna"、最も多いのが"dna_rm"です。また、最も"ACGTN"以外の文字が多いのが"dna_sm"ということになります。</p>
<br>
  - [Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.nonchromosomal.fa.gz](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.nonchromosomal.fa.gz)  
  - [Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa.gz](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa.gz)  
  - [Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna_rm.nonchromosomal.fa.gz](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna_rm.nonchromosomal.fa.gz)  
  - [Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna_rm.toplevel.fa.gz](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna_rm.toplevel.fa.gz)  
  - [Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna_sm.nonchromosomal.fa.gz](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna_sm.nonchromosomal.fa.gz)  
  - [Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna_sm.toplevel.fa.gz](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna_sm.toplevel.fa.gz)  
<br>
<p>soft-maskedがわかりづらいかもしれませんが、これは塩基配列の表記法(nomenclature)はA, C, G, T, Nの5種類だけではないという事実を正しく理解するところからスタートせねばなりません。たとえば[プリン塩基](https://ja.wikipedia.org/wiki/%E3%83%97%E3%83%AA%E3%83%B3%E5%A1%A9%E5%9F%BA)(A or G)であることが確かなら、1文字表記のRとします。同様に、[ピリミジン塩基](https://ja.wikipedia.org/wiki/%E3%83%94%E3%83%AA%E3%83%9F%E3%82%B8%E3%83%B3%E5%A1%A9%E5%9F%BA)(C or T)であることが確かならYです。他にもありますので、全体像を以下にまとめておきます。なお、元情報は1984年に勧告された命名法の論文([Cornish-Bowden A., Nucleic Acids Res., 1985](https://pubmed.ncbi.nlm.nih.gov/2582368))中のTable 1です。</p>
| 塩基 | 1文字表記 | 表記の由来 |
|:---:|:---:|:---:|
| A or C | M | a**M**ino |
| A or G | R | pu**R**ine |
| A or T | W | **W**eak interaction (2 H bonds) |
| C or G | S | **S**trong interaction (3 H bonds) |
| C or T | Y | p**Y**rimidine |
| G or T | K | **K**etone |
| A or C or G | V | not-T (not-U), **V** follows U  in the alphabet |
| A or C or T | H | not-G, **H** follows G in the alphabet |
| A or G or T | D | not-C, **D** follows C in the alphabet |
| C or G or T | B | not-A, **B** follows A in the alphabet |
| A or C or G or T | N | a**N**y |
<br>
先程まで[GFF3](https://en.wikipedia.org/wiki/General_feature_format)ファイルである[Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.49.gff3](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.49.gff3)とともに読み込んで、転写物配列ファイル([R1.130.2.fasta](R/R1.130.2.fasta))作成の元データとして用いたのは、[Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa](R/Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa)です。以下ではまず、[Biostrings](http://bioconductor.org/packages/release/bioc/html/Biostrings.html)パッケージが提供する`readDNAStringSet`関数が[gzip](https://ja.wikipedia.org/wiki/Gzip)圧縮ファイルでも読み込めることを示します。`DNAStringSeq`形式の`fasta`オブジェクトの`width`列と`names`列を見比べることで、`width`列の配列長で降順(多-->少)にコンティグ(contig; ctg)が並べられていることがわかります。  
```{r, message=FALSE, error=TRUE}
in_f <- "Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa.gz"
library(Biostrings)
fasta <- readDNAStringSet(in_f, format="fasta")
fasta
```
<br>
ちなみに、`width(fasta)`とやることで、配列(この場合はコンティグ)ごとの配列長情報を取得することができます。計75配列あるので、この出力結果は75個の要素からなる数値ベクトルです。1番左側に[1]や[11]や[21]が見えていますが、これは要素番号情報です。たとえば、表示されている数値ベクトルの1番目の要素が362642、11番目の要素が84732、21番目の要素が54844なのだと解釈します。  
```{r, message=FALSE, error=TRUE}
width(fasta)
```
<br>
[R1.100 1.10.2項のGFF3(ファイル取得)](https://www.iu.a.u-tokyo.ac.jp/textbook/#R1100_1102%E9%A0%85%E3%81%A7%E8%BF%B0%E3%81%B9%E3%81%9FGFF3%E3%83%95%E3%82%A1%E3%82%A4%E3%83%AB%E3%82%92%E5%8F%96%E3%82%8A%E6%89%B1%E3%81%86(%E3%83%87%E3%83%BC%E3%82%BF%E5%8F%96%E5%BE%97%E7%B7%A8))のところで、「③より、このゲノムサイズは3,215,878 bp」と書いています。この数値は、以下のように計75配列の配列長が格納された数値ベクトルである`width(fasta)`を入力として、ベクトルの総和を計算する`sum`関数を実行することで確認できます。  
```{r, message=FALSE, error=TRUE}
sum(width(fasta))
```
<br>
`DNAStringSeq`形式の`fasta`オブジェクトを入力として`alphabetFrequency`関数を実行すると、`fasta`オブジェクトの各配列に対して、塩基の種類ごとの出現頻度を返してくれます。ただ、もし`alphabetFrequency(fasta)`としてしまうと、A, C, G, T, Nだけでなく、上記の表で示したような「プリン塩基(A or G)ならR」、「ピリミジン塩基(C or T)ならY」といったような様々な結果が配列ごとに一気に表示されてしまい見づらくなってしまいます。そのため、以下では`baseOnly=TRUE`というオプションをつけて「A, C, G, Tとそれ以外」という5種類に要約して表示させています。  
```{r, message=FALSE, error=TRUE}
alphabetFrequency(fasta, baseOnly=TRUE)
```
<br>
なお、`baseOnly=TRUE`オプションをつけずに`alphabetFrequency(fasta)`を実行した結果を、そのまま`colSums`の入力として与えることで、入力配列全体のプリン塩基やピリミジン塩基の数がわかります。  
```{r, message=FALSE, error=TRUE}
in_f <- "Lactobacillus_casei_21_1_gca_000309585.ASM30958v1.dna.toplevel.fa.gz"
library(Biostrings)
fasta <- readDNAStringSet(in_f, format="fasta")
colSums(alphabetFrequency(fasta))
```
<br>

## R関連情報(第2章)  

### R2.010 2.1.3項の組合せ計算
page041の1行目で計算している171通りというのは、`choose`関数で実行できます。  
```{r, message=FALSE, error=TRUE}
choose(19, 17)
```

### R2.020 2.2.2項の図2.4
page047の[図2.4](chapter2/fig2.4.pdf)aの再現を行います。まずは「(a)実データ」です。-log(p-value)の値が若干異なります(教科書では0.9248で、以下では0.9249234になっていると思います)。理由は不明ですが(単なる転記ミス?!)、まあ誤差範囲といえるでしょう。  
```{r, message=FALSE, error=TRUE}
data <- matrix(c(222, 86, 149, 228, 105, 124), nrow=3)
data
chisq.test(data)$p.value
-log10(chisq.test(data)$p.value)
```
<br>
次は「(a)仮想データ1」です。  
```{r, message=FALSE, error=TRUE}
data <- matrix(c(222, 106, 129, 228, 105, 124), nrow=3)
data
chisq.test(data)$p.value
-log10(chisq.test(data)$p.value)
```
<br>
次は「(a)仮想データ2」です。  
```{r, message=FALSE, error=TRUE}
data <- matrix(c(202, 86, 169, 228, 105, 124), nrow=3)
data
chisq.test(data)$p.value
-log10(chisq.test(data)$p.value)
```
<br>
次は、表の中にはありませんが、本文中で述べている「実データを基本として全部10倍の罹患群4,570人 vs. 対照群4,570人」です。  
```{r, message=FALSE, error=TRUE}
data <- matrix(c(2220, 860, 1490, 2280, 1050, 1240), nrow=3)
data
chisq.test(data)$p.value
-log10(chisq.test(data)$p.value)
```
<br>
最後に「(b)」です。  
```{r, message=FALSE, error=TRUE}
data <- matrix(c(287, 149, 21, 292, 150, 15), nrow=3)
data
chisq.test(data)$p.value
-log10(chisq.test(data)$p.value)
```

### R2.030 2.4.3項の表2.1(準備)
page055の[表2.1](chapter2/table2.1.xlsx)は、[Bioconductor](https://www.bioconductor.org/)から提供されている[BSgenome.Hsapiens.UCSC.hg19](https://bioconductor.org/packages/release/data/annotation/html/BSgenome.Hsapiens.UCSC.hg19.html)(ヒトゲノム配列のパッケージ)と[TxDb.Hsapiens.UCSC.hg19.knownGene](https://bioconductor.org/packages/release/data/annotation/html/TxDb.Hsapiens.UCSC.hg19.knownGene.html)(上記ヒトゲノム配列パッケージに対応したアノテーション情報のパッケージ)を利用して作成しています。ここではまず、これら2つのパッケージをインストールします。圧縮して提供されているため、ヒトゲノムの一般的なイメージである3GBまではないにせよ、700MB程度はあります。したがって、インストールにはそれなりの時間がかかります(光のWifi環境で約5分とか...)。  
```{r, message=FALSE, error=TRUE, warning=FALSE, eval=FALSE}
BiocManager::install("BSgenome.Hsapiens.UCSC.hg19", update=F)
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene", update=F)
```
<br>

### R2.040 2.4.3項の表2.1a
ここではまず、[表2.1](chapter2/table2.1.xlsx)aの1番左側の列で見えている2-merの並び情報を`narabi`オブジェクトとして定義し、必要なパッケージをロードします。  
```{r, message=FALSE, error=TRUE}
narabi <- c("AA", "AT", "TA", "TT", "CC", "CG", "GC", "GG",
            "AC", "AG", "CA", "CT", "GA", "GT", "TC", "TG")
library(Biostrings)                 #パッケージのロード
library(BSgenome.Hsapiens.UCSC.hg19)#パッケージのロード
```
<br>
以下は、さきほどロードした[BSgenome.Hsapiens.UCSC.hg19](https://bioconductor.org/packages/release/data/annotation/html/BSgenome.Hsapiens.UCSC.hg19.html)と同名のオブジェクトを`genome`という名前に変更したのち、`getSeq`関数で`DNAStringSet object`形式のオブジェクトとして抽出した結果を`fasta`に保存するまでのスクリプトです。最後に`fasta`オブジェクトの配列名情報を与えるべく、`seqnames(genome)`で代入しています。  
```{r, message=FALSE, error=TRUE}
#前処理(指定したパッケージ中のオブジェクト名をgenomeに統一)
genome <- BSgenome.Hsapiens.UCSC.hg19  #renameしているだけです
fasta <- getSeq(genome)                #ゲノム塩基配列情報を抽出した結果をfastaに格納
names(fasta) <- seqnames(genome)       #description情報を追加している
fasta                                  #確認してるだけです
```
<br>
以下は、[表2.1](chapter2/table2.1.xlsx)の2列目の情報を得る部分のスクリプトです。まず、`fasta`を入力として`alphabetFrequency`関数を実行し、配列ごとのACGT各塩基の出現頻度をカウントした結果を`hoge`に格納しています。`hoge`自体は、298行×18列からなる数値行列です。1行目がchr1、2行目がchr2といった感じで、配列ごとの結果が各行に格納されています。各列にはA, C, G, T, …の並びで塩基ごとのカウント結果が格納されています。それゆえ、`colSums`関数を用いて列ごとの総和を計算してやれば、`fasta`オブジェクト全体のAの数、Cの数、Gの数、Tの数の情報が得られます。実行結果である`frequency`オブジェクトはベクトルであり、1番目の要素がAの数、2番目の要素がCの数、といった感じで格納されています。それゆえ、`sum(frequency[2:3])/sum(frequency[1:4])`のような計算をしてあげれば、GC含量の情報(= 0.4101656)が得られるのです。`expected`オブジェクトを得る部分は、はあくまでも理解(説明)のしやすさを優先させています。リアルの研究や投稿論文のSupplementなどで提供する場合は、直前で得た`GCcontent`オブジェクトの情報を使うなどして一般化したコーディングを行いますし、行ってください。
```{r, message=FALSE, error=TRUE}
#本番(expected)
hoge <- alphabetFrequency(fasta)       #A,C,G,T,..の数を配列ごとにカウントした結果をhogeに格納
frequency <- colSums(hoge)             #列ごとの総和をfrequencyに格納
GCcontent <- sum(frequency[2:3])/sum(frequency[1:4])#GC含量を計算(0.4101656)
expected <- c(rep(0.295*0.295, 4), rep(0.205*0.205, 4), rep(0.205*0.295, 8))#べた書きでnarabiの並びで期待値を算出
expected                               #確認してるだけです
```
<br>
以下は、[表2.1](chapter2/table2.1.xlsx)の3列目の情報を得る部分のスクリプトです。まず、`fasta`を入力として`dinucleotideFrequency`関数を実行し、配列ごとの2連続塩基の出現頻度をカウントした結果を`hoge`に格納しています。次に、`colSums`関数を用いて`hoge`オブジェクトの列ごとの総和を計算しています。実行結果である`frequency`オブジェクトはベクトルであり、1番目の要素がAAの数、2番目の要素がACの数、といった感じで格納されています。これは出現頻度情報なので、`frequency/sum(frequency)`の計算をしてあげることで、出現確率情報に変換することができます。それが`probability`オブジェクトを得ている部分です。最後に`narabi`オブジェクトの並びにしたいので、`probability[narabi]`を実行して順番を入れ替えています。
```{r, message=FALSE, error=TRUE}
#本番(observed)
hoge <- dinucleotideFrequency(fasta, as.prob=F)#連続塩基の出現頻度情報をhogeに格納
frequency <- colSums(hoge)             #列ごとの総和をfrequencyに格納
probability <- frequency/sum(frequency)#出現確率を計算した結果をprobabilityに格納
observed <- probability[narabi]        #narabiの順番に変更
observed                               #確認してるだけです
```
<br>
以下は、[表2.1](chapter2/table2.1.xlsx)a全体の情報を得るスクリプトです。シンプルに3つのベクトル(`expected`, `observed`, and `observed/expected`)を`cbind`関数を用いて列方向に結合(column bind)しているだけです。
```{r, message=FALSE, error=TRUE}
table2.1a <- cbind(expected, observed, observed/expected)
head(table2.1a)                        #最初の6行分を表示
```
<br>

### R2.050 2.4.3項の表2.1b
ここではまず、転写開始点(TSS)上流1,000塩基の領域を調べたいので、その数値情報を`param_upstream`オブジェクトに与えています。次に、[TxDb.Hsapiens.UCSC.hg19.knownGene](https://bioconductor.org/packages/release/data/annotation/html/TxDb.Hsapiens.UCSC.hg19.knownGene.html)(アノテーション情報のパッケージ)をロードして、パッケージ名と同じ`TxDb`という形式のオブジェクトを`txdb`という名前に変更(rename)しています。`txdb`の表示結果の詳細を理解する必要はなく、「こんな感じに見えていたらアノテーション情報が`txdb`オブジェクトの中にちゃんと格納されているのだという程度でよい」です。
```{r, message=FALSE, error=TRUE}
param_upstream <- 1000                 #転写開始点上流の塩基配列数を指定
library(TxDb.Hsapiens.UCSC.hg19.knownGene)#パッケージのロード
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene #renameしているだけです
txdb                                   #確認してるだけです
```
<br>
以下は、転写開始点(TSS)上流1,000塩基の配列情報を得る部分のスクリプトです。まず、`txdb`オブジェクトから遺伝子の座標情報を染色体の座標情報順で取得した結果を`gn`オブジェクトに保存しています。次に、`gn`オブジェクトを入力として、5'上流領域(5' flanking region)の座標情報を取得するための`flank`関数を実行することで、上流1,000塩基分の座標情報を取得しています。そして、上流1,000塩基分の範囲情報からなる`hoge`オブジェクトを、`getSeq`関数実行時のオプションとして与えています。`fasta`実行結果の配列長(`width`列)が全て1,000塩基であることから、なんとなくうまくとれてるっぽいことがわかります。
```{r, message=FALSE, error=TRUE}
#本番
gn <- sort(genes(txdb))                #遺伝子の座標情報を取得
hoge <- flank(gn, width=param_upstream)#指定した範囲の座標情報を取得
fasta <- getSeq(genome, hoge)          #指定した範囲の塩基配列情報を取得
fasta                                  #確認してるだけです
```
<br>
以下は、[表2.1](chapter2/table2.1.xlsx)bの1列目の情報を得る部分のスクリプトです。まず、`fasta`を入力として`alphabetFrequency`関数を実行し、配列ごとのACGT各塩基の出現頻度をカウントした結果を`hoge`に格納しています。`hoge`自体は、298行×18列からなる数値行列です。1行目がchr1、2行目がchr2といった感じで、配列ごとの結果が各行に格納されています。各列にはA, C, G, T, …の並びで塩基ごとのカウント結果が格納されています。それゆえ、`colSums`関数を用いて列ごとの総和を計算してやれば、つまり`fasta`オブジェクト全体のAの数、Cの数、Gの数、Tの数の情報が得られます。実行結果である`frequency`オブジェクトはベクトルであり、1番目の要素がAの数、2番目の要素がCの数、といった感じで格納されています。それゆえ、`sum(frequency[2:3])/sum(frequency[1:4])`のような計算をしてあげれば、GC含量の情報(= 0.518723)が得られるのです。`expected`オブジェクトを得る部分は、はあくまでも理解(説明)のしやすさを優先させています。リアルの研究や投稿論文のSupplementなどで提供する場合は、直前で得た`GCcontent`オブジェクトの情報を使うなどして一般化したコーディングを行いますし、行ってください。
```{r, message=FALSE, error=TRUE}
#本番(expected)
hoge <- alphabetFrequency(fasta)       #A,C,G,T,..の数を配列ごとにカウントした結果をhogeに格納
frequency <- colSums(hoge)             #列ごとの総和をfrequencyに格納
GCcontent <- sum(frequency[2:3])/sum(frequency[1:4]) #GC含量を計算(0.518723)
expected <- c(rep(0.24*0.24, 4), rep(0.26*0.26, 4), rep(0.26*0.24, 8))#べた書きでnarabiの並びで期待値を算出
expected                               #確認してるだけです
```
<br>
以下は、[表2.1](chapter2/table2.1.xlsx)bの2列目の情報を得る部分のスクリプトです。まず、`fasta`を入力として`dinucleotideFrequency`関数を実行し、配列ごとの2連続塩基の出現頻度をカウントした結果を`hoge`に格納しています。次に、`colSums`関数を用いて`hoge`オブジェクトの列ごとの総和を計算しています。実行結果である`frequency`オブジェクトはベクトルであり、1番目の要素がAAの数、2番目の要素がACの数、といった感じで格納されています。これは出現頻度情報なので、`frequency/sum(frequency)`の計算をしてあげることで、出現確率情報に変換することができます。それが`probability`オブジェクトを得ている部分です。最後に`narabi`オブジェクトの並びにしたいので、`probability[narabi]`を実行して順番を入れ替えています。
```{r, message=FALSE, error=TRUE}
#本番(observed)
hoge <- dinucleotideFrequency(fasta, as.prob=F)#連続塩基の出現頻度情報をhogeに格納
frequency <- colSums(hoge)             #列ごとの総和をfrequencyに格納
probability <- frequency / sum(frequency)#出現確率を計算した結果をprobabilityに格納
observed <- probability[narabi]        #narabiの順番に変更
observed                               #確認してるだけです
```
<br>
以下は、[表2.1](chapter2/table2.1.xlsx)b全体の情報を得るスクリプトです。シンプルに3つのベクトル(`expected`, `observed`, and `observed/expected`)を`cbind`関数を用いて列方向に結合(column bind)しているだけです。
```{r, message=FALSE, error=TRUE}
table2.1b <- cbind(expected, observed, observed/expected)
head(table2.1b)                        #最初の6行分を表示
```
<br>

### R2.060 2.4.3項の表2.1(出力)
R2.040とR2.050で得た2つのオブジェクト(`table2.1a`と`table2.1b`)を、`narabi`の情報を含めてファイル([table2.1.txt](R/table2.1.txt))に出力するスクリプトです。
```{r, message=FALSE, error=TRUE}
####ファイルに保存
out_f <- "table2.1.txt"               #出力ファイル名を指定してout_fに格納
tmp <- cbind(narabi, table2.1a, table2.1b)#保存したい情報をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
```

### R2.070 2.4.3項の例題2.2
[例題2.2](chapter2/E2.2.pdf)の内容も、表示させたい桁数を指定するやり方の提供を兼ねて、念のためここでも示しておきます。[Saxonov論文](https://pubmed.ncbi.nlm.nih.gov/16432200)の[Table 1](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1345710/table/tbl1/)より、1 kb upstream resionsのGC含量は0.53、Observed CpG fractionは0.042であることがわかります。まずは「(1)このGC含量(0.53)を用いてCCの期待値を小数第4位まで示せ」を行います。  
```{r, message=FALSE, error=TRUE}
# CCの期待値
GCcontent <- 0.53                     #GC含量
C <- GCcontent/2                      #Cの期待値はGC含量の1/2

#結果の基本形
C*C                                   #CCの期待値

#表示桁数を指定
param_digit <- 4                      #表示させたい桁数
round(C*C, digits=param_digit)        #小数点第(param_digit + 1)位で四捨五入した結果
```  

次に、「(2) このGC含量(0.53)を用いてGAの期待値を小数第4位まで示せ」を行います。  
```{r, message=FALSE, error=TRUE}
GCcontent <- 0.53                     #GC含量
G <- GCcontent/2                      #Gの期待値はGC含量の1/2
A <- (1 - GCcontent)/2                #Aの期待値は(1 - GC含量)の1/2

#結果の基本形
G*A                                   #GAの期待値

#表示桁数を指定
param_digit <- 4                      #表示させたい桁数
round(G*A, digits=param_digit)        #小数点第(param_digit + 1)位で四捨五入した結果
```  

次に「(3) 上記1kb upstream regionsの情報(つまりGC含量は0.53、Observed CpG fractionは0.042)を用いてCGの観測値/期待値を少数第3位まで示せ」をやります。  
```{r, message=FALSE, error=TRUE}
param_digit <- 3                      #表示させたい桁数
expected <- C*G                       #CGの期待値
observed <- 0.042                     #問題文から与えられているCGの観測値
round(observed/expected, digits=param_digit)#小数点第(param_digit + 1)位で四捨五入した結果
```  

### R2.080 2.4.4項の性能評価指標
page057の8-9行目で述べている「...混同行列を内部的に作成し、それをもとに上記の感度や特異度などを算出して提示する...」を[CRAN](https://cran.r-project.org/)から提供されている[caret](https://CRAN.R-project.org/package=caret)というRパッケージを用いて実践します。[caret](https://CRAN.R-project.org/package=caret)は、機械学習や回帰系のRパッケージとして有名です。その原著論文([Khun M., J Stat Software, 2008](https://doi.org/10.18637/jss.v028.i05))の引用回数は、[7,000](https://scholar.google.co.jp/scholar?cites=803974403145666124&as_sdt=2005&sciodt=0,5&hl=ja)弱です(2023年9月18日調べ)。ここではまず、パッケージをインストールします。  
```{r, message=FALSE, error=TRUE, warning=FALSE, eval=FALSE}
BiocManager::install("caret", update=F)
```  
 
以下は本番のコードです。まず、ヒトゲノム上のある範囲に計5つの領域が存在し、それぞれの領域はCpGアイランド(CGI)またはそれ以外(non-CGI)のいずれかのラベルが付与されていると仮定します(これが真実の情報)。この範囲に存在する計5つの領域の塩基配列情報をそれぞれ与えて、CGIかどうかを判定させた結果が手元にあったとします(これが予測の結果)。この要素数が同じ2つのベクトル(`truthdayo`と`predicted`)を入力として[caret](https://CRAN.R-project.org/package=caret)パッケージが提供する`confusionMatrix`関数を実行すると、正確度(Accuracy)や感度(Sensitivity)などの結果が得られます。  
```{r, message=FALSE, error=TRUE, warning=FALSE}
truthdayo <- c("CGI", "non-CGI", "CGI", "CGI", "non-CGI")#真実の情報
predicted <- c("CGI", "CGI",     "CGI", "CGI", "CGI")    #予測の結果
library(caret)                                #パッケージのロード
confusionMatrix(data = factor(predicted),     #予測の結果はdataのところ
                reference = factor(truthdayo),#真実の情報はreferenceのところ
                positive = "CGI")             #Positiveを定義
```  

以下は、よりシンプルに書いたスクリプトです。文字列ベクトルではなく、0 or 1の数値ベクトルを入力として与えています。`confusionMatrix`関数は"CGI" or "non-CGI"のラベル情報をアルファベット順にソートするため、対応するラベルを"0" or "1"にして`positive = "0"`をオプションとして与えることで、得られる混同行列の左右の列を意図通りに並べています。
```{r, message=FALSE, error=TRUE, warning=FALSE}
truthdayo <- c(0, 1, 0, 0, 1)                 #真実の情報
predicted <- c(0, 0, 0, 0, 0)                 #予測の結果
library(caret)                                #パッケージのロード
confusionMatrix(data = factor(predicted),     #予測の結果はdataのところ
                reference = factor(truthdayo),#真実の情報はreferenceのところ
                positive = "0")               #Positiveを定義
```  

以下は、感度を下げて特異度を上げるべく、予測結果を変えたスクリプトです。  
```{r, message=FALSE, error=TRUE, warning=FALSE}
truthdayo <- c(0, 1, 0, 0, 1)                 #真実の情報
predicted <- c(1, 1, 1, 0, 1)                 #予測の結果
library(caret)                                #パッケージのロード
confusionMatrix(data = factor(predicted),     #予測の結果はdataのところ
                reference = factor(truthdayo),#真実の情報はreferenceのところ
                positive = "0")               #Positiveを定義
```  

