# データ解析入門2018

## 0日目

### 解析環境の構築

Windows Subsystem for Linux (WSL)のインストール

以下の点を確認する

1. Windowsを最新版にアップデート (2018年9月だとバージョン1803であること！ 確認方法: キーボードのWindowsキーを押しながら、[R]キーを押し、[ファイル名を指定して実行]画面が表示されるので、[名前]欄に「winver」と入力し、[OK]ボタンを選択します。)

1. WSLのインストール
次のURL等を参照
https://qiita.com/yoshikingt/items/ab86411e6031459db805

1. 日本語化：
[[WSLの日本語化]]を参照

## 1日目

### AWK入門

#### Hello World!
#### 変数
#### if文
#### 演算子
#### FASTQ, FASTAのリード数カウント

資料：{{2018AWK入門1.pptx}}

[[2018AWK1回答例]]


## 2日目

### AWK入門2

#### for文
#### 配列
#### 配列のソート

資料：{{2018AWK入門2.pptx}}

[[2018AWK2回答例]]

## おまけ

### AWKスクリプトの作成を補助してくれるツール
[[AWKスクリプトの作成を補助してくれるツール|詳細]]
### Ubuntuをグラフィカルに使いたい場合(リモートデスクトップでUbuntuに接続)
[[Ubuntuをグラフィカルに使いたい場合|詳細]]

## 3日目

### AWK入門3

#### FASTA -> Tab区切り形式への変換
#### 正規表現

資料：{{2018AWK入門3.pptx}}

[[2018AWK3回答例]]

## 4日目

### AWK入門4

#### 複数ファイルの処理
#### Blast結果を発現情報にアノテーション

資料：{{2018AWK入門4.pptx}}

[[2018AWK4回答例]]

## 5日目

### AWK入門5

#### 正規表現を用いた置換
#### VCFファイル整形

資料：{{2018AWK入門5.pptx}}

[[2018AWK5回答例]]

## 6日目

### AWK入門6

#### エスケープ処理
#### 注意が必要な文字を含むテキストの処理

資料：{{2018AWK入門6.pptx}}

[[2018AWK6回答例]]

## 7日目

### AWK入門7

#### 二次元配列
#### アンプリコンデータ解析

資料：{{2018AWK入門7.pptx}}

[[2018AWK7回答例]]

## 8日目

### AWK入門8

#### 関数
#### GOツリー構築

資料：{{2018AWK入門8.pptx}}

[[2018AWK8回答例]]

## 9日目

### Excel入門

#### 絶対参照・相対参照
https://kokodane.com/kihon19.htm

##### 練習問題1

下記のようにサンプルごとに各遺伝子の発現量の表があるとする。この表はRNA-seqのリード数をそのまま表示したRAWデータだとする。これを補正してみる。

| |sample1|sample2|sample3|sample4|
|gene1|178|199|666|569|
|gene2|7|790|471|17|
|gene3|116|516|216|904|
|gene4|222|720|471|345|

1.まずはサンプル毎にリード数が違うとサンプル間の比較ができないので、サンプル毎に合計リード数を算出する。(sum関数を使用する)

2.次に各サンプルの遺伝子ごとに次のようなRPM (Reads Per Million mapped reads)という値に変換した表を作成せよ。

```
( (各サンプルの遺伝子ごとのリード数) / (各サンプルの合計リード数) ) * 1000 * 1000
```

---

#### テキストファイルをExcelで読み込む
http://kouritu.net/csv-excel-error-prevent/

#### フィルター
https://www.forguncy.com/blog/20170516_filter1

#####練習問題2

4日目に行ったblastの情報を発現量データにアノテーションした結果ファイルを下記のコマンドでダウンロードする。

```
wget http://www.suikou.fs.a.u-tokyo.ac.jp/yosh_data/2018train/take.txt
```

1.フィルターを使用して、「pearl」を含む遺伝子の発現量のみ表示せよ。

2.さらに、その中で遺伝子発現量の平均値が高い順に表示せよ。

---

#### 条件付き書式
http://www4.synapse.ne.jp/yone/excel/excel_syosiki_jyouken_databar.html

http://cblog.crie.jp/excel/290/

#####練習問題3

練習問題2で使用したファイルに条件付き書式を使用して、発現量の大小が一目でわかるようにせよ。

---

#### 関数、マクロ、VBA

https://dekiru.net/article/4429/

https://www.excelspeedup.com/macrosettei/

https://kokodane.com/macro_kouza.htm

#####練習問題4

練習問題2で使用したファイルにマクロを使用して、1～50行までの間、1行ずつ条件付き書式を使用して、発現量の大小が一目でわかるようにせよ。

回答例：{{take-macro.xlsm}}

---

#### 練習問題5

VCFデータの整形

```
wget http://www.suikou.fs.a.u-tokyo.ac.jp/yosh_data/2018train/sample.txt
```

上記ファイルは、黒色系統の生物(b)と黒色系統の中から生まれた白色の劣性突然変異体(w)のジェノタイプ情報である。

1.まずは上記ファイルをExcelで読み込む。

2.ホモ(1/1)、ヘテロ(0/1)、変異なし(0/0)、デプスが低くてジェノタイプ不明(./.)を分かりやすく色付けせよ。

3.白色の遺伝様式は劣性の遺伝様式を示すことが知られている。そのため、白色になる原因変異の近傍は、ホモ接合となっていることが想定される。そこで白色系統の個体でホモ(1/1および0/0)、黒色系統の個体でヘテロ(0/1)となる割合をそれぞれ計算せよ。(ヒント：countif関数を使うと便利)

4.候補変異がありそうな領域はどこか？

回答例：{{sample.xlsx}}

## 10日目

### バイオインフォ入門(マッピング編)

|FASTAファイル相同性検索|[[https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=Download|blast]], [[http://last.cbrc.jp/|last]]|
|ゲノム(ショートリード)|[[https://github.com/lh3/bwa|bwa]], [[http://bowtie-bio.sourceforge.net/bowtie2/index.shtml|bowtie2]]|
|ゲノム(ロングリード)|[[https://github.com/lh3/minimap2|minimap2]], [[https://github.com/ocxtal/minialign|minialign]]| 
|RNA-seq|[[https://ccb.jhu.edu/software/tophat/index.shtml|tophat]], [[http://ccb.jhu.edu/software/hisat2/index.shtml|hisat]]|

#### ツールのインストール 

上記のツールは比較的メジャーなため、実行ファイルをダウンロードするだけで良いときもあるが、Perl, Python, AWKなどのスクリプトファイルになっていることも多く、それらのインタープリタをインストールする必要がある。また、自分でソースファイルから実行ファイルをビルドする(作成する)必要があることもある。今回ツールとは別に必要な環境は下記のコマンドでインストール出来る。

```
sudo apt update
sudo apt install build-essential unzip python zlib1g-dev
```

ツールごとにインストール方法は異なるので、ツールのマニュアルをよく読んでインストールすること。

マニュアルはツールのHPに纏められていたり、ダウンロードしたファイルの中にREADME, MANUALといった名前のファイルに書いてあったりする。もしくは、ツールのオプションとして、「-h」、「--help」などをつけて実行したら出てくる説明文を読んだほうが早い場合もある。マニュアルを探す順番は、本家WEBサイト→ダウンロードしたファイルの中→とりあえず実行して表示されるメッセージ→本家以外のサイトをググる　という順番が無難。ただ慣れないうちは、日本語で使い方が書いてあるページをググってみるのも良いかも。本家サイト以外では情報が古いことが多く、現在のバージョンと違いがある可能性を考慮すること。

ダウンロードしたファイルは基本的には圧縮されているため、フォーマットに適した解凍を行う必要がある。

|ファイル名|解凍方法|
|xxx.tar.gz|tar zvxf xxx.tar.gz|
|xxx.gz|gzip -d xxx.gz|
|xxx.tar.bz2|tar jvxf xxx.tar.bz2|
|xxx.bz2|bzip2 -d xxx.bz2|
|xxx.tar|tar vxf xxx.tar|
|xxx.zip|unzip xxx.zip|

macの人は、自分でビルドすると苦労することが多いので、https://bioconda.github.io/index.html#set-up-channels
の手順でbiocondaをインストールしてから、https://bioconda.github.io/recipes.html
で希望のパッケージを見つけてインストールするのが無難。

---

#### ツールの比較

マッピングツールを動かし、特徴をつかむ練習。
下記の配列データがショウジョウバエのゲノムのどこに由来する配列なのか、正しく結論を出してください。

FASTQ: {{read.fq}}

FASTA: {{read.fa}} (blast, last用)

ショウジョウバエのゲノム配列を下記からダウンロードする。

http://hgdownload.cse.ucsc.edu/goldenPath/dm6/bigZips/dm6.fa.gz

準備した配列は次のようになっています。
|配列名|ヒント|見てほしいところ|
|unique-seq1|ゲノム中に1か所だけ||
|low-complexity-seq1|かなりATリッチ|blastでは・・・|
|centromere-seq1|セントロメアの一部配列（かなり多コピー)||
|unique-seq1-copy, low-complexity-seq1-copy, centromere-seq1-copy|上記配列と全く同じ配列|配列は同じでもヒットする箇所は・・・|
|rnaseq-no-intron|RNA-seqでイントロンを挟まない配列||
|rnaseq-intron|RNA-seqでイントロンを挟む配列||
|p53_partial|ショウジョウバエのp53遺伝子|イントロンがあります|
|nanopore|ナノポアのデータ(低クオリティ、ロングリード)||

大まかには次のキーワードに注意して挙動の違いを見て下さい。

- 同じスコアで複数箇所にヒットした場合

- 何回か実行してみて同じ結果が得られるか

- グローバルアライメントとローカルアライメント

- イントロンを含むリード(大きなdeletionがあるように見える）

####出力ファイル形式

SAMフォーマット

http://crusade1096.web.fc2.com/sam.html

tophat2ではそれをSAMを圧縮したBAMファイルが出力される。それを見るにはsamtoolsというツールが必要だが、

```
sudo apt install samtools
```

とやればインストール可能で、使用するときは、

```
samtools view tophat_out/accepted_hits.bam
```

などとすればよい。またFASTAファイル中の特定の範囲を抜き出して見たい場合、例えばchr2Lの1-10番目の塩基を表示させるには、

```
samtools faidx dm6.fa chr2L:1-10
```

[[2018mapping例]]

## 11日目

### R入門

#### Rのインストール

1. R本体のインストール
 
 https://cran.r-project.org/

2. RStudioのインストール
 
 https://www.rstudio.com/products/rstudio/download/

#### ライブラリのインストール

Rの本体だけではあまり大したことは出来ないので、研究者の人たちが作成したライブラリをインストールして使用することになる。

ライブラリのインストールは各ライブラリのHPの通りに行えばよいが、基本的には下記の2通りで行うことになる。

1. CRANから(例：veganをインストールする場合)
 
```
install.packages("vegan")
```

2. Bioconductorから(例：DESeq2をインストールする場合)

```
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DESeq2", version = "3.8")
```

インストールが終わったら、いずれの方法でインストールしても、

```
library("DESeq2")
```

などとライブラリ名を打てばライブラリがロードされる。

#### Rの用途

データ解析におけるRの主な用途は、統計検定やグラフ作成になる。グラフ作成時には、どのようなグラフをRで描くことが可能か、下記のようなサイトで調べると良い。

- https://www.r-graph-gallery.com/

- https://www.imsbio.co.jp/RGM/top

- google画像検索

- https://github.com/rstudio/cheatsheets/blob/master/data-visualization-2.1.pdf

#### R言語入門

R言語をざっくり読めるようになるには、下記のサイトなどが良さそう。

https://qiita.com/stkdev/items/6aba2c1db2fa056170ae

#### Rで変数の中身を理解するには

- RStudioで変数名をクリック

- class() 関数

- str() 関数

#### ①普通に作図：散布図の作成

サンプルデータ{{mariom.tpm.txt}}は、0h vs 24hの比較で有意に発現上昇or減少した遺伝子の0h,24h,48h,・・・のデータ。まずは普通に散布図を描いてみる。

```
data=read.table("mariom.tpm.txt",header=T,row.names=1,sep="\t")
plot(x=data[,1],y=data[,2])
```

色々なパラメータ例：http://cse.naro.affrc.go.jp/takezawa/r-tips/r/53.html

練習

・data[,1], data[,2]をlog2()で変換して対数の値で表示

・RNA-seqで良く見るMAプロットを描いてみる
```
MAプロットはX軸にA, Y軸にMを描いたもの
M = log2(サンプル1) – log2(サンプル2)
A =  { log2(サンプル1) + log2(サンプル2) } / 2
```

・0h、24hのサンプルをそれぞれ1つ選んで描いてみる

#### ②ggplot2を使って作図：BOX AND SCATTER PLOT

http://www.r-graph-gallery.com/89-box-and-scatter-plot-with-ggplot2/

```
# Ggplot2 library
library(ggplot2)
 
# Data 
names=c(rep("A", 80) , rep("B", 50) , rep("C", 70))
value=c( sample(2:5, 80 , replace=T) , sample(4:10, 50 , replace=T), 
       sample(1:7, 70 , replace=T) )
data=data.frame(names,value)
 
#Graph
qplot( x=names , y=value , data=data , fill=names) + geom_boxplot() + geom_jitter()
```

・①のmariom.tpm.txtでサンプルを3つ、遺伝子を100個選んで、そのboxplotを描いてみる。

#### ③ヒートマップ作成

```
data=read.table("mariom.tpm.txt",header=T,row.names=1,sep="\t")
library(Heatplus)
data=as.matrix(data)
corrdist = function(x) as.dist(1-cor(t(x)))
hclust.avl = function(x) hclust(x, method="average")
plot(regHeatmap(data, legend=2, dendrogram = list(clustfun=hclust.avl, distfun=corrdist)))
```

・遺伝子の数が多すぎてヒートマップが見えにくいので、20個だけ使用してヒートマップを描く