# 遺伝学研究室実習 2024
2023年度 学生実習 2024/12/12〜2024/12/20
# ＜目的＞
シロイヌナズナの遺伝学とゲノミクス解析により、エピゲノム制御に関わる因子の機能とメカニズムをあきらかにする研究を体験し、遺伝学的手法や遺伝学の考え方、エピゲノムについて学ぶ。
# ＜背景＞
シロイヌナズナのIBM1 (INCREASE IN BONSAI METHYLATION1)遺伝子は、ヒストン脱メチル化酵素をコードし、ヒストンH3の9番目のリジンのメチル化 (H3K9me2) の除去に働いている。通常、H3K9me2はトランスポゾンのような抑制されている配列に存在しているが、IBM1欠損変異体においては、活性に転写されている遺伝子に「異所的に」H3K9me2が蓄積し、同時に非CG型DNAメチル化（non-CGme）も蓄積する。このことは、IBM1は遺伝子領域からH3K9me2とnon-CGmeを排除する働きを持っていることを示している。さらに、IBM1欠損変異体では植物の発生異常、稔性の低下などが引き起こされる。遺伝子領域のH3K9me2/non-CGmeの蓄積による発生異常の仕組みを理解するために、ibm1変異体のサプレッサー変異が同定された。サプレッサー変異とは、ある変異体で起きている現象が、その変異体にさらに変異を加えることでもとに戻る（野生型に近づく）ものであり、遺伝学の分野ではメカニズムの理解に一役買っている。得られたサプレッサー変異の一つは、LDL2（LYSINE SPECIFIC DEMETHYLASE-LIKE2）の機能欠損変異であり、動物のホモログの機能から推測するにLDL2はH3K4の脱メチル化酵素であると考えられた。
本実習では、*ibm1 ldl2*二重変異体の解析を通して、遺伝子内領域におけるエピゲノム修飾間の相互作用と、植物の表現型とのつながりを理解することを目指す。
# ＜実習内容＞
1: PCRによる*IBM1*と*LDL2*の遺伝子型判定と植物の表現型観察
2: メチル化感受性制限酵素を用いた*ibm1 ldl2*二重機能欠損変異体におけるDNAメチル化解析
3: 全ゲノムDNAメチル化データの解析
4: ChIP-seqとRNA-seqデータの解析
5: 実験・解析結果の発表
# ＜日程＞
| 月曜日 | 火曜日| 水曜日| 木曜日| 金曜日|
| -------- | -------- | -------- | -------- | -------- |
||||12/12 イントロ、遺伝子型判定(1.1~1.3)|12/13 DNA抽出(2.1)、DNA濃度測定(2.2)、制限酵素処理(2.3)|
|12/16 メチル化感受性PCR(2.4)、電気泳動(2.5)|12/17 メチル化感受性PCR(2.4)、電気泳動(2.5)、データ解析(3)|12/18 実習おやすみ |12/19 データ解析(3)|12/20 データ解析(3)、発表準備、結果発表会|
# ＜実験1. 遺伝子型判定＞
**<目的>** *ibm1*変異体背景で*ldl2*変異をヘテロで持つ個体(*ibm1*/*ibm1* *LDL2*/*ldl2*)の自殖(自家受粉によって種をとること)によって得られた種をまき、出てくる個体の中で、*LDL2*の遺伝子型(*LDL2*/*LDL2*, *LDL2*/*ldl2*, *ldl2*/*ldl2*の３種類)によって、植物の表現型がどう変わるかを調べる。
*IBM1*と*LDL2*遺伝子の遺伝子型 (genotype)を決定するためのPCR (genotyping; ジェノタイピングという)を行う。シロイヌナズナでは、T-DNAのランダムな挿入により遺伝子が破壊された様々な変異体のストックが存在する。また、研究にはEMS処理により誘導された突然変異体も用いられており、遺伝子型の判定方法は変異の様式によって異なる。
本実習で用いる変異体はT-DNA挿入変異体であるが、代表的な遺伝子型判定方法を以下に示す。
**○PCRによる遺伝子型判定方法**
**・T-DNA挿入変異体**
T-DNAをランダムに挿入した変異体の整備が進められており、シロイヌナズナの研究者コミュニティで共有されている。本実習で用いる変異体はどちらもT-DNA挿入変異体である。
実験1の遺伝子型判定では、複数のプライマーのミックスを用いて、両染色体にT-DNA挿入がある個体(mutant)、片方の染色体のみにT-DNA挿入がある個体(hetero)、T-DNA挿入がない個体(wild type; WT)を判別する。**どのようなプライマーを設計すれば、３つの状態が判別できるか**を考えてみてください。
**・点突然変異体**
EMSなどの変異原処理では、主に点突然変異が誘導される。変異が制限酵素の認識配列上に誘導された場合、PCR産物を制限酵素処理してDNA断片長を区別する(CAPS: cleaved amplified polymorphic sequences)。制限酵素の認識配列以外での点変異については、制限酵素の認識サイトを形成するように改変したプライマーを点変異に隣接するように設計して、制限酵素処理したDNA断片長を区別する(dCAPS: derived CAPS)。
**・挿入/欠失(indel)変異体**
indelを挟むように設計したプライマーでPCRし、増幅したDNAの断片長により区別する。
**◯試薬**
DNA抽出バッファー、TEバッファー、EmeraldAmp (PCR反応試薬、TaKaRa)、50 X TAEバッファー、MidoriGreen (DNA染色試薬)、アガロース
**◯器具**
1.5 mlチューブ、~~爪楊枝~~、チューブラック、アイスボックス、電気泳動装置、ゲル作成トレイとコーム、0.2 ml PCRチューブ、サーマルサイクラー、卓上遠心機
**◯今回使用する変異体**
*ibm1-4*, *ldl2-2* ----全てT-DNA挿入変異体
**◯実験手法**
## 1.1 DNA抽出
シロイヌナズナの葉からDNAを簡易法により抽出する。
本実習では**各自8個体ずつ**(**内1個体は野生型植物をコントロールとして行う**)抽出し、以下の解析に用いる。
1) 0.2 ml PCRチューブに、DNA抽出バッファーを10μlずつ加える。
2) 植物の葉を切り取り(米粒程度)、バッファーに浸るようにチューブに入れる(サンプルの量が多すぎると、不純物が多くなりPCRが回らなくなるので注意)。
~~3) 爪楊枝で植物をすりつぶす。~~
4) PCRチューブのふたをしっかり閉めて軽く遠心したあと、サーマルサイクラーを用いて95℃で10分間処理する。
5) 軽く遠心したあと、TEバッファーを50μlずつ加える。
6) よく混ぜて、軽く遠心する。
## 1.2 PCRによる遺伝子型判定
1.1で抽出したDNAを鋳型としてPCRを行う。本実験ではIBM1とLDL2の2種類の遺伝子について調べる。
酵素の失活を防ぐために、2xEmeraldAmpは氷上に保つ。
1) 以下の組成でプレミックスを作成する。2種類のプライマーセット用に2本用意する。
||使用量|サンプル数分 (x9)|
|---|---|---|
|プライマーmix(各10μM)|0.4 μl||
|[2xEmeraldAmp](https://catalog.takara-bio.co.jp/product/basic_info.php?unitid=U100006935)|5 μl||
|H2O|||
|**total**|**9μl**||
|鋳型DNA (分注後に加える)|1 μl|-|
2) 0.2 ml PCRチューブにプレミックスを9µlずつ分注する。
3) 0.2 ml PCRチューブに鋳型DNA溶液を1μlずつ加える。
4) PCRチューブの蓋をして、スピンダウンする。
5) 以下の条件でPCRを行う。
①98℃ 1分
②98℃ 10秒(熱変性)
58℃ 20秒(アニーリング)
72℃ 60秒(伸長反応)
→35サイクル
③72℃ 1分
使用するプライマーセット
*ibm1-4*判別用プライマーセット
*ldl2-2*判別用プライマーセット
## 1.3 アガロースゲル電気泳動
1) 2% アガロースを作製するため、三角フラスコに1 X TAEバッファー100 mlとアガロース2g、マグネティックスターラーを入れ混和し、電子レンジでアガロースを溶解する。スターラー上で撹拌しながら冷却し、手で持てる程度(50-70℃)に冷めたら、ゲルトレイに5mm厚になるように流し入れ固める。固まったらゲルの上に1 X TAEを少量かけ、コームをゆっくりとはずし、1 X TAEを満たした泳動層に沈める。ゲルが残ったら、蒸発しないように蓋をして恒温槽で保管する。
2) 100bpマーカー2 μlを端のウェルにアプライする。続いて、隣のウェルに漏れないように注意しながら、サンプルを5 μlずつアプライする。（ibm1とldl2のPCR産物をそれぞれ8つ並べてアプライした方があとから結果を確認しやすい）
4) 100Vで20分電気泳動し、MidoriGreen染色液で染色する。
5) UVトランスイルミネーターで確認し、写真を撮る。
# ＜実験2. メチル化感受性酵素を用いたDNAメチル化レベルの解析＞
**DNAメチル化を解析する方法の代表例**
・バイサルファイトシークエンス法 ☆現在はもっぱらこれが使われる.
・メチル化感受性制限酵素を使用した方法 (MSRE-PCR、下図参照) ☆本実験で使用する
・メチル化DNA特異的切断酵素（McrBC）を使用した方法（McrBC-PCR）
・メチル化DNA特異抗体を使用した免疫沈降（MeDIP）
図2.1 制限酵素を用いたDNAメチル化解析の概要 (TaKaRaホームページより)
**◯注意事項**
酵素類を取り扱う際は、以下の点に注意する。
* ボルテックスしない
* 手で酵素部位を持たない
* 低温を保ち、使用後はすぐに冷凍庫に戻す
* 酵素をピペットマンでとるときは、チップの先端だけを酵素液につけ、ゆっくり吸い取る
**◯試薬**
*Msp* I(NEB)、*Bgl* II (NEB)、それぞれの10X制限酵素バッファー、TE、EmeraldAmp (PCR反応試薬、TaKaRa)、50 X TAEバッファー、10X ローディングバッファー、MidoriGreen (DNA染色試薬)、アガロース
**◯器具**
2.0 mlチューブ、チューブラック、アイスボックス、電気泳動装置、ゲル作成トレイとコーム、0.2 ml PCRチューブ、サーマルサイクラー、37℃インキュベーター、 Nanodrop
**◯実験手法**
## **2.1 DNA抽出**
Nucleon Phytopure(GE Healthcare)を用いて、高品質のDNAを抽出する。
1) 3~4枚のロゼット葉を2mlチューブに取り、破砕用ビーズを一粒加えて液体窒素で凍結する。サンプル名は蓋ではなく側面に記入する。
2) Automill(電動破砕機)でサンプルを破砕する。1500rpm、1分を3回繰り返す。
3) Reagent Iを300μl加えて、よく混ぜる。RNase Aを1μl加えて、よく混ぜた後、37℃で30分間インキュベートする。
4) Reagent IIを100μl加えてよく混ぜたあと、65℃で10分間インキュベートする。
5) チューブを氷上に移し、20分間インキュベートする。
6) チューブを軽く遠心し、ビーズを残して溶液をピペットマンで新しい1.5mlエッペンチューブに移す。
7) クロロホルムを250μl加えてボルテックスで激しく混ぜる。室温で10分間振盪する。
8) 10000 x g、4℃で10分間遠心する。
9) 水層を新しいエッペンチューブに回収後、400μlの氷冷イソプロパノールを加えてよく混ぜる.
10) 10000 x g、4℃で5分間遠心する。
11) 上清を除き、70%エタノールを500μl加える。
12) 10000 x g、4℃で5分間遠心する。遠心後、上清をピペットマンで除き、チューブの蓋を開けたまま5分間風乾する。
13) 40μlの滅菌水を加えて、タッピングによりDNAを溶解する。
14) DNAは冷蔵保存する。
## **2.2 濃度測定**
Nanodrop(微量分光光度計)を使用してDNAの濃度を調べる。
1) Nanodropを起動する。電極をSDWで濡らしたキムワイプ、乾いたキムワイプの順で拭う。
2) 設定からDNAを選択し、サンプルの溶媒でブランクをとる(今回は滅菌水)。
下の電極に滅菌水を2 μlのせ、上の電極をゆっくりおろしブランクを取る。
3) 測定したブランク（サンプル）をキムワイプで拭い取り、次のサンプルをのせ吸光度を測定する(DNAの吸光度 50 μg/ mlで、A260=1.0)。
## **2.3 制限酵素処理**
1) 各DNA400ngをPCRチューブに入れ、10 ng/μlの濃度になるようにSDWで希釈する(40μlに調整) 。
||濃度(ng/μl)|400ng DNA(µl)|H2O (µl)|
| --------| -------- | -------- | -------- |
|WT||||
|*ibm1*||||
|*ibm1 ldl2*||||
|*ibm1 ldl2/+*||||
各制限酵素の認識配列は、
*Msp* I: 5'-**C**CGG-3'
*Bgl* II: 5'-AGAT**C**T-3'
であり、どちらもnon-CGメチル化に対して感受性(太字のCがメチル化されていると切断されない)
2) DNA以外の試薬で以下の組成になるようにプレミックスを調製する。酵素は分注の直前に加える(表記は1サンプル分：20 μlスケール。4サンプル+α分をそれぞれ計算する)。
| *Msp* I |使用量| Mixture(x4.5) | | *Bgl* II| 使用量|Mixture(x4.5)|
| -------- | --|-------- | --|-------- |-------- |---|
| CutSmart Buffer|2µl|||NEBuffer 3.1|2µl||
|*Msp* I |0.5µl|||*Bgl* II |0.5µl ||
| H2O | | | |H2O |||
| **pre-mix total** | | | |**pre-mix total** |||
| DNA (50ng) | |プレミックス 分注後に加える| |DNA (50ng) | |プレミックス 分注後に加える|
|**total**|**20µl**|**-**||**total**|**20µl**|**-**|
3) DNA以外のプレミックスをPCRチューブに分注し、DNAを必要量加える。
4) PCRチューブに蓋をし、混合し、スピンダウンする。PCRチューブの蓋がしっかりしまっていることを確認（蓋の相性が悪いとしまらずに蒸発してしまう！）
5) 37℃でO/N(一晩)酵素処理をする。
**注意）酵素使用後はすぐに冷凍庫へ。** 残りのDNAは冷蔵保存。
## **2.4 PCR**
DNAメチル化状態を調べる遺伝子のプライマーセットと酵素処理したDNAを使用して、PCRにより目的断片を増幅する。調べる遺伝子毎にプレミックスを作ると良い。今回はPCRのコントロールとして、DNAが高度にメチル化されているCACTAトランスポゾン配列を増幅させる。
実験方法は、1.2を参照のこと。
今回増幅させる配列は以下の通り
|遺伝子|ID|制限酵素|
|--|--|--|
|BONSAI (BNS)|At1g73177|*Bgl* II|
|ERECTA-LIKE2 (ERL2)|At5g07180|*Msp* I|
|ACTIN2 (ACT2)|At3g18780|*Msp* I|
|CACTA (トランスポゾン)|At2g12210|*Bgl* II, *Msp* I (control)|
1) FとRのプライマー（オリジナルは100 µM）を希釈・混合し、プライマーmix (各10μM)を作る
||使用量|
|---|---|
|Fプライマー(100 µM)|10 μl|
|Rプライマー(100 µM)|10 μl|
|H2O|80 µl |
|**total**|**100μl**|
2) 鋳型DNA以外のPCRプレミックス（4+αサンプル分）を作り、PCRチューブに分注する。
||使用量|Mixture (x4.5)|
|---|---|---|
|プライマーmix (各10μM)|0.4 μl||
|2xEmeraldAmp|5 µl ||
|H2O|3.6 µl||
|**total**|**9μl**||
3) 鋳型DNA 1 µl を加え、蓋をし、スピンダウンする。
4) PCRは以下の条件で行う。
①98℃ 1分
②98℃ 10秒(熱変性)
56℃ 20秒(アニーリング)
72℃ 30秒(伸長反応)
→30サイクル
③72℃ 2分
## **2.5 電気泳動**
1.3と同様の手順で行う。
# ＜実験3. 全ゲノムDNAメチル化データ（BS-seq・WGBS）の解析＞
## 3-1. BS-seqデータのマッピングとDNAメチル化の集計(今回はスキップ)
Bisulfite処理したゲノムDNAを全ゲノムシークエンシング(BS-seq)して得られたリードを、(1)シロイヌナズナのゲノム上にマッピングし、(2)遺伝子ごとのメチル化レベルを集計する。マッピングにはBismarkというプログラムを使用する(Krueger & Andrews (2011) Bioinformatics 27: 1571-1572)。3-1の操作はmacで行うことを前提として書かれている。
### 3-1-0. 下準備
3-1の操作はターミナルを使って、コマンドラインで行う。扱うファイルが大きい（1ファイル数Gb）ので、クラスターコンピューターに接続して計算する。
**ターミナルを使ってみよう**
ターミナルのアイコンをクリックして起動させる(サーバへの接続はこちらでやります)。今回使用するデータは『Practice_2021/BS_seq』内に保存されている。データが格納されているディレクトリへ移動する。
**・ターミナル上での基本的な操作（Unix commands)**
```
ls #(list segments): ディレクトリ内のファイルを表示
cd Practice_2021 #(change directory): 指定したディレクトリ(この場合Practice_2021)に移動
cd .. #ひとつ上のディレクトリに移動（この場合ホームディレクトリ（beowulf））
```
次のように打ち込んで『BS_seq』ディレクトリに移動する。
```
cd Practice_2021
ls
cd BS_seq
```
lessコマンドでファイルの中身が表示される。
『q』を打つと元の画面に戻る。リードデータ(.fastqファイル)の中身を見てみよう。fastqは大規模シークエンス（次世代シーケンス）で得られるリードファイルの一般的な形式で、以下の写真のような形式となっている。
`less xxxx.fastq`
リードデータ(.fastq)の一般的な構成は、
① リードのID (@HWI-ST…)
② 塩基配列 (AAATTGAAATTGTT…)
③ 区切り記号+
④ 信頼性 (CCCFFFFFHHHH…)
### 3-1-1.リードのマッピングと各シトシンのメチル化レベルの計算
☆[ ]の部分は必要に応じて変えてください。
☆[遺伝子名]のところには、変異体の名称に変えてください。
**・リードのマッピング**
Bisulfite処理によって、メチル化されていないシトシンはチミンとして読まれるので、シトシンとチミンを区別せずに（４塩基ではなく、３塩基の情報として）シロイヌナズナゲノムにマッピングする。
https://www.bioinformatics.babraham.ac.uk/projects/bismark/
```bash
nohup ~/source/bismark_v0.15.0/bismark --bowtie1 -n 1 -l 20 -e 90 ~/check_genome/tair10_a_thaliana/ -1 [遺伝子名]_5M_1.fastq -2 [遺伝子名]_5M_2.fastq &
```
と打ち込み、”Bismark”を実行する。進行状況は、『ps –x』と打ち込むと確認出来る。
マッピングが終わると、新しいファイルが出来る(ファイル名は、[gene_name]_1_5M_paired.fastq_bismark_pe.bam)。
**・重複配列の排除**
```bash
~/source/bismark_v0.15.0/deduplicate_bismark [遺伝子名]_5M_1.fastq_bismark_pe.bam &
```
**・シトシンごとのメチル化レベルの計算**
```bash
nohup ~/source/bismark_v0.15.0/bismark_methylation_extractor --bedGraph --CX --cytosine_report --genome_folder ~/check_genome/tair10_a_thaliana/ -s [遺伝子名]_5M_1.fastq_bismark_pe.deduplicated.sam &
```
### 3-1-2. 遺伝子ごとのメチル化レベルの集計
上記のBismark解析によって、xxxx.CX_report.txtというファイルができた。このファイルには、全ゲノム中のシトシン1つ1つが、それぞれどの程度メチル化されているかが書き込まれている。このファイルをもとにして、シロイヌナズナの各遺伝子ごとに、メチル化されたシトシンとされていないシトシンの数を集計する。
まずは、Bismarkの解析の結果、『...CX_report.txt』が作成されていることを確認する。lessコマンドを使って中身を見てみる。
```
ls
less ....CX_report.txt
```
①染色体番号
②染色体上の位置
③ストランド(+ or -)
④シトシンの数(メチル化シトシン)
⑤チミンの数(非メチル化シトシン)
⑥コンテクスト
⑦実際の配列
次に遺伝子の位置情報を格納したファイルを見てみる。
`less tair10_models_startsort-2.txt`
自家製perlプログラム「methyl_C_count.pl」に、CX_reportファイルと遺伝子の位置情報ファイルを与えて実行する。
```bash
perl methyl_C_count.pl [遺伝子名]_1_paired.fastq_bismark_pe.deduplicated.CX_report.txt &
```
計算が終わると、同じディレクトリ内に、『 [遺伝子名]....CX_report.**gene**.txt』という新しいファイルが出来る。
`less [遺伝子名]....CX_report.gene.txt`
このようなファイルが出来上がるはず。
**・できたファイルのダウンロード**
遺伝子ごとにメチル化レベルを集計したファイルをサーバーから自分のコンピュータに移動。TAが操作します。
```bash
scp -P xxxx beowulf@xxx.xx.x.xx:Practice_2019/fastq_5M/[遺伝子名]_5M_1.fastq_bismark_pe.deduplicated.CX_report.gene.txt ./
```
**・ファイル名の変更** (長いので)
```bash
mv [遺伝子名]_5M_1.fastq_bismark_pe.deduplicated.CX_report.gene.txt [遺伝子名]_meC.txt
```
[出来上がりはこちらからダウンロード](https://drive.google.com/drive/folders/11dqyFO8FptAyQeo-OJsJdf8_8UHRLMa7?usp=sharing)
## 3-2.RStudioによるデータ可視化
### 準備
Rをインストール
https://www.r-project.org/
R tips; Rの基本的な使い方を網羅したウェブサイト *閉鎖されてしまった.
http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html
RStudioをインストール
https://rstudio.com/
RStudioを起動する。
まず、File > New projectから新規プロジェクトを作成。プロジェクトを保存したディレクトリ（PC中のフォルダ）が作業ディレクトリとなる。RStudioを閉じた後にプロジェクトに戻りたい場合はFile > Open projectからファイルを選んで開く。
以下の写真のように、Import Dataset > From Text (base)から、ダウンロードしたメチル化ファイルを取り込む。データフレーム(Data frame)と呼ばれるエクセルシートのようなものとして取り込まれる。（"Heading"のところをYesにすることを忘れないように！）
①エディタ画面
②コンソール画面
③オブジェクト一覧
④プロット画面
（配置はRStudio > Preferences > Pane Layoutで変更できる）
### 様々な可視化
**・散布図**
```r
plot(x= [メチル化レベルを集計したData frame名]$[列の名前],y= [Data frame名]$[列の名前])
```
### **オプション**
**・トランスポゾンを赤くする**
```r
plot(x= [xx]$[xx], y=[xx]$[xx], col=ifelse([Data frame名]$te==1,rgb(1,0,0),rgb(0,0,0)))
```
**・色も様々に変えられる**
以下は色指定の例; red, green, blueの順で三原色を0−1の数字で表す。
もしくは"red", "blue", "cyan"などのように色の名前で指定することもできる。ダブルクォーテーションマーク（半角）を忘れずに。
**・シンボルの形を変えるオプション**
```r
pch=20
```
指定する数字により以下のシンボルが使える。
**・シンボルの大きさを変更するオプション**
```r
cex=0.5
```
**・x軸、y軸のラベルはxlab="WT mCG", ylab="ibm1 mCG"のように変えられる**
**・特定の遺伝子を別の色で表す**
```r
plot(x= Col_meC$CHGm, y=ibm1_meC$CHGm,col=rgb(0,0,0,0.3),pch=16,cex=0.5,xlim=c(0,1),ylim=c(0,1),xlab="",ylab="")
par(new=T)
plot(x= Col_meC$CHGm[Col_meC$gene=="AT5G07180"], y=ibm1_meC$CHGm[Col_meC$gene=="AT5G07180"], col=rgb(1,0,0,1),pch=16,cex=1,xlim=c(0,1),ylim=c(0,1))
```
**・散布図の例**
```r
plot(x= Col_meC$CGm, y=ibm1_meC$CGm, col=ifelse(Col_meC$te==1,rgb(1,0,0,0.3),rgb(0,0,0,0.3)),pch=16,cex=0.5)
plot(x= Col_meC$CGm, y=ibm1_ldl2_meC$CGm, col=ifelse(Col_meC$te==1,rgb(1,0,0,0.3),rgb(0,0,0,0.3)),pch=16,cex=0.5)
plot(x= ibm1_meC$CGm, y=ibm1_ldl2_meC$CGm, col=ifelse(Col_meC$te==1,rgb(1,0,0,0.3),rgb(0,0,0,0.3)),pch=16,cex=0.5)
plot(x= Col_meC$CHGm, y=ibm1_meC$CHGm, col=ifelse(Col_meC$te==1,rgb(1,0,0,0.3),rgb(0,0,0,0.3)),pch=16,cex=0.5)
plot(x= Col_meC$CHGm, y=ibm1_ldl2_meC$CHGm, col=ifelse(Col_meC$te==1,rgb(1,0,0,0.3),rgb(0,0,0,0.3)),pch=16,cex=0.5)
plot(x= ibm1_meC$CHGm, y=ibm1_ldl2_meC$CHGm, col=ifelse(Col_meC$te==1,rgb(1,0,0,0.3),rgb(0,0,0,0.3)),pch=16,cex=0.5)
```
**・描いてみよう**
[1] WT（横軸）とibm1（縦軸）でCGメチル化を比較する（TEと遺伝子を別の色で）
[2] WT（横軸）とibm1（縦軸）でCHGメチル化を比較する(TEと遺伝子を別の色で)
[3] ibm1（横軸）とibm1 ldl2（縦軸）でCHGメチル化を比較する(TEと遺伝子を別の色で)
etc...
**・特定の集団を抽出**
```r
[新しいdata frame名] <- subset([抽出元のdata frame名],[条件式])
```
<- (小なりとハイフンで矢印を表現)は右を左に代入するということ。
条件式の例
te == 1
("te"列の値が1に等しい)
mCG >= 0.5
("mCG"列の値が 0.5以上)
参照
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/28.html
**・ヒストグラム**
CGメチル化レベルの分布を見る。breaksを変えると刻みが変えられる。
```r
hist(Col_meC$CGm[Col_meC$te==1], col = "#0000ff99", border = "#0000ff", breaks =seq(0,1,0.1))
```
**・ボックスプロット(箱ひげ図)**
遺伝子 (te列が0)とトランスポゾン遺伝子（te列が1）を分けてCGメチル化レベルのボックスプロットを書く例
```r
boxplot(Col_meC$CGm[Col_meC$te==0],Col_meC$CGm[Col_meC$te==1], names=c("GENE", "TE"), col=c("#993435", "#539952"), main="Boxplot", ylim=c(0, 1), ylab="CG methyl")
```
CGメチル化がある(>1%)とCGメチル化がない(<1%)遺伝子(te=0)を分けてCHGメチル化レベルのボックスプロットを書く例
```r
boxplot(Col_meC$CHGm[Col_meC$CGm>0.01 & Col_meC$te==0],Col_meC$CHGm[Col_meC$CGm<=0.01 & Col_meC$te==0], names=c("CGm high", "CGm low"), col=c("#993435", "#539952"), main="Boxplot", ylim=c(0, 0.1), ylab="CHG methyl")
```
*ibm1*や*ibm1 ldl2*変異体では？
ボックスプロットの各要素（四角の上端、下端、太線、破線など）の意味
https://data-science.gr.jp/implementation/ida_r_boxplot.html
**・要素を分類した上でのボックスプロット**
各遺伝子（TE遺伝子を含む）をCGメチル化レベルにより10個のグループに分類した上で、各グループのCHGメチル化レベルをボックスプロットで示す。CGメチル化レベルとCHGメチル化レベルにどんな関係があるのか？
```r
plot(c(0,9),c(0,1),type="n",xlab="CG_methylation",ylab="CHG_methylation") #c(0,9)->分類するグループの数。c(0,1)->y軸の最小値、最大値
for(i in 0:9){ #分類するグループの数と同じにする。このfor文で各グループのboxplotを順に出力していく。
temp_data<-subset(Col_meC,CGm>0.1*(i)&CGm<0.1*(i+1)) #CG_meレベルを0~0.1,0.1~0.2,0.2~0.3,0.3~0.4…という感じで分けて一時的にグループを作製
boxplot(temp_data$CHGm,outline=F,at = c(i),add=TRUE) #CHGメチル化レベルを出力 #atでplotの場所指定。addで図の中に追加 outline=Fは外れ値を出力しない。
}
rm(temp_data,i) #グループの消去
```
### クラスタリングをしてみる
野生型と*ibm1*, *ibm1 ldl2*変異体でのDNAメチル化データをもとに、挙動が似ている遺伝子同士をまとめるクラスタリングを行う。k-meansクラスタリング法と階層的クラスタリング法を行ってみる。（*階層的クラスタリングはコンピューターのスペックが足りないと動かないので、スキップ可）
**以下 k-meansクラスタリングのやり方**
```r=
#野生型と変異体のWGBSの結果を使って、クラスタリング解析を行う。
#すべての株におけるデータを統合
All_meC <- Col_meC[,c(1,5,6,7,8,9,10)]
All_meC$Col_CGm <- Col_meC$CGm
All_meC$Col_CHGm <- Col_meC$CHGm
All_meC$Col_CHHm <- Col_meC$CHHm
All_meC$ibm1_CGm <- ibm1_meC$CGm
All_meC$ibm1_CHGm <- ibm1_meC$CHGm
All_meC$ibm1_CHHm <- ibm1_meC$CHHm
All_meC$ibm1_ldl2_CGm <- ibm1_ldl2_meC$CGm
All_meC$ibm1_ldl2_CHGm <- ibm1_ldl2_meC$CHGm
All_meC$ibm1_ldl2_CHHm <- ibm1_ldl2_meC$CHHm
#メチル化の列だけを抜き出し。
data<-All_meC[,8:16]
#k-meansクラスタリング。クラスターの数はここでは5に設定。変更可。
km <- kmeans(data,5)
#クラスターの可視化
DC<-data[order(km$cluster),]
heatmap(as.matrix(DC),Rowv = NA,Colv=NA)
#遺伝子、偽遺伝子とトランスポゾンがどのようにクラスタリングされたかをみる
TEorGENE <- All_meC$te #gene,0; TE, 1; pseudogene, 2
result_kmc <- km$cluster
table_kmc <- table(TEorGENE, result_kmc)
table_kmc
#もとデータフレームにクラスター番号列を追加。
All_meC$kmcluster <- result_kmc
#クラスターごとに箱ひげ図を描いてみる。
boxplot(All_meC$Col_CGm[All_meC$kmcluster==1],All_meC$Col_CGm[All_meC$kmcluster==2],All_meC$Col_CGm[All_meC$kmcluster==3],All_meC$Col_CGm[All_meC$kmcluster==4],All_meC$Col_CGm[All_meC$kmcluster==5],
All_meC$Col_CHGm[All_meC$kmcluster==1],All_meC$Col_CHGm[All_meC$kmcluster==2],All_meC$Col_CHGm[All_meC$kmcluster==3],All_meC$Col_CHGm[All_meC$kmcluster==4],All_meC$Col_CHGm[All_meC$kmcluster==5],
All_meC$Col_CHHm[All_meC$kmcluster==1],All_meC$Col_CHHm[All_meC$kmcluster==2],All_meC$Col_CHHm[All_meC$kmcluster==3],All_meC$Col_CHHm[All_meC$kmcluster==4],All_meC$Col_CHHm[All_meC$kmcluster==5],
cex=0,col = c("blue","gray","forest green","purple","magenta"),ylim=c(0,1),ylab="methylation",main="Wild type")
boxplot(All_meC$ibm1_CGm[All_meC$kmcluster==1],All_meC$ibm1_CGm[All_meC$kmcluster==2],All_meC$ibm1_CGm[All_meC$kmcluster==3],All_meC$ibm1_CGm[All_meC$kmcluster==4],All_meC$ibm1_CGm[All_meC$kmcluster==5],
All_meC$ibm1_CHGm[All_meC$kmcluster==1],All_meC$ibm1_CHGm[All_meC$kmcluster==2],All_meC$ibm1_CHGm[All_meC$kmcluster==3],All_meC$ibm1_CHGm[All_meC$kmcluster==4],All_meC$ibm1_CHGm[All_meC$kmcluster==5],
All_meC$ibm1_CHHm[All_meC$kmcluster==1],All_meC$ibm1_CHHm[All_meC$kmcluster==2],All_meC$ibm1_CHHm[All_meC$kmcluster==3],All_meC$ibm1_CHHm[All_meC$kmcluster==4],All_meC$ibm1_CHHm[All_meC$kmcluster==5],
cex=0,col = c("blue","gray","forest green","purple","magenta"),ylim=c(0,1),ylab="methylation",main="ibm1")
boxplot(All_meC$ibm1_ldl2_CGm[All_meC$kmcluster==1],All_meC$ibm1_ldl2_CGm[All_meC$kmcluster==2],All_meC$ibm1_ldl2_CGm[All_meC$kmcluster==3],All_meC$ibm1_ldl2_CGm[All_meC$kmcluster==4],All_meC$ibm1_ldl2_CGm[All_meC$kmcluster==5],
All_meC$ibm1_ldl2_CHGm[All_meC$kmcluster==1],All_meC$ibm1_ldl2_CHGm[All_meC$kmcluster==2],All_meC$ibm1_ldl2_CHGm[All_meC$kmcluster==3],All_meC$ibm1_ldl2_CHGm[All_meC$kmcluster==4],All_meC$ibm1_ldl2_CHGm[All_meC$kmcluster==5],
All_meC$ibm1_ldl2_CHHm[All_meC$kmcluster==1],All_meC$ibm1_ldl2_CHHm[All_meC$kmcluster==2],All_meC$ibm1_ldl2_CHHm[All_meC$kmcluster==3],All_meC$ibm1_ldl2_CHHm[All_meC$kmcluster==4],All_meC$ibm1_ldl2_CHHm[All_meC$kmcluster==5],
cex=0,col = c("blue","gray","forest green","purple","magenta"),ylim=c(0,1),ylab="methylation",main="ibm1 ldl2")
#比較しやすくするためクラスターごとに並べた場合(Cコンテクストごとに3つのパネル)
boxplot(All_meC$Col_CGm[All_meC$kmcluster==1],All_meC$ibm1_CGm[All_meC$kmcluster==1],All_meC$ibm1_ldl2_CGm[All_meC$kmcluster==1],All_meC$Col_CGm[All_meC$kmcluster==2],All_meC$ibm1_CGm[All_meC$kmcluster==2],All_meC$ibm1_ldl2_CGm[All_meC$kmcluster==2],
All_meC$Col_CGm[All_meC$kmcluster==3],All_meC$ibm1_CGm[All_meC$kmcluster==3],All_meC$ibm1_ldl2_CGm[All_meC$kmcluster==3],All_meC$Col_CGm[All_meC$kmcluster==4],All_meC$ibm1_CGm[All_meC$kmcluster==4],All_meC$ibm1_ldl2_CGm[All_meC$kmcluster==4],
All_meC$Col_CGm[All_meC$kmcluster==5],All_meC$ibm1_CGm[All_meC$kmcluster==5],All_meC$ibm1_ldl2_CGm[All_meC$kmcluster==5],
cex=0,col = c("blue","blue","blue","gray","gray","gray","forest green","forest green","forest green","purple","purple","purple","magenta","magenta","magenta"),ylim=c(0,1),ylab="methylation",main="mCG")
boxplot(All_meC$Col_CHGm[All_meC$kmcluster==1],All_meC$ibm1_CHGm[All_meC$kmcluster==1],All_meC$ibm1_ldl2_CHGm[All_meC$kmcluster==1],All_meC$Col_CHGm[All_meC$kmcluster==2],All_meC$ibm1_CHGm[All_meC$kmcluster==2],All_meC$ibm1_ldl2_CHGm[All_meC$kmcluster==2],
All_meC$Col_CHGm[All_meC$kmcluster==3],All_meC$ibm1_CHGm[All_meC$kmcluster==3],All_meC$ibm1_ldl2_CHGm[All_meC$kmcluster==3],All_meC$Col_CHGm[All_meC$kmcluster==4],All_meC$ibm1_CHGm[All_meC$kmcluster==4],All_meC$ibm1_ldl2_CHGm[All_meC$kmcluster==4],
All_meC$Col_CHGm[All_meC$kmcluster==5],All_meC$ibm1_CHGm[All_meC$kmcluster==5],All_meC$ibm1_ldl2_CHGm[All_meC$kmcluster==5],
cex=0,col = c("blue","blue","blue","gray","gray","gray","forest green","forest green","forest green","purple","purple","purple","magenta","magenta","magenta"),ylim=c(0,1),ylab="methylation",main="mCHG")
boxplot(All_meC$Col_CHHm[All_meC$kmcluster==1],All_meC$ibm1_CHHm[All_meC$kmcluster==1],All_meC$ibm1_ldl2_CHHm[All_meC$kmcluster==1],All_meC$Col_CHHm[All_meC$kmcluster==2],All_meC$ibm1_CHHm[All_meC$kmcluster==2],All_meC$ibm1_ldl2_CHHm[All_meC$kmcluster==2],
All_meC$Col_CHHm[All_meC$kmcluster==3],All_meC$ibm1_CHHm[All_meC$kmcluster==3],All_meC$ibm1_ldl2_CHHm[All_meC$kmcluster==3],All_meC$Col_CHHm[All_meC$kmcluster==4],All_meC$ibm1_CHHm[All_meC$kmcluster==4],All_meC$ibm1_ldl2_CHHm[All_meC$kmcluster==4],
All_meC$Col_CHHm[All_meC$kmcluster==5],All_meC$ibm1_CHHm[All_meC$kmcluster==5],All_meC$ibm1_ldl2_CHHm[All_meC$kmcluster==5],
cex=0,col = c("blue","blue","blue","gray","gray","gray","forest green","forest green","forest green","purple","purple","purple","magenta","magenta","magenta"),ylim=c(0,0.3),ylab="methylation",main="mCHH")
```
**最適なクラスター数の決め方(elbow method) for k-means**
``` r=
WCSS<-c()
for(i in 1:10){
WCSS[i]<-kmeans(data,i)$tot.withinss
}
plot(1:10,WCSS)
```
**最適なクラスター数が見つかったら、上記のk-meansクラスタリングのクラスター数を変えてやり直してみる**
**以下、階層的クラスタリングのやり方（階層的クラスタリングはコンピューターのスペックが足りないと動かないので、スキップ可）**
```r=
#階層的クラスタリング。少し時間とコンピューターパワーが必要。
# 非類似度（距離）を計算
distance <- dist(data) #ユークリッド距離を求める
# 樹形図作成
hc <- hclust(distance, "ward.D2")
# デンドログラム プロット。結果をみて、何個のクラスターに分類するかを決める。
plot(hc,label=F)
# 5つのクラスターに分けることにする。
result_hc <- cutree(hc,k=5)
table_hc <- table(TEorGENE, result_hc)
table_hc
#もとデータフレームにクラスター番号列を追加。
All_meC$hcluster <- result_hc
#クラスターごとに箱ひげ図を描いてみる。
boxplot(All_meC$Col_CGm[All_meC$hcluster==1],All_meC$Col_CGm[All_meC$hcluster==2],All_meC$Col_CGm[All_meC$hcluster==3],All_meC$Col_CGm[All_meC$hcluster==4],All_meC$Col_CGm[All_meC$hcluster==5],
All_meC$Col_CHGm[All_meC$hcluster==1],All_meC$Col_CHGm[All_meC$hcluster==2],All_meC$Col_CHGm[All_meC$hcluster==3],All_meC$Col_CHGm[All_meC$hcluster==4],All_meC$Col_CHGm[All_meC$hcluster==5],
All_meC$Col_CHHm[All_meC$hcluster==1],All_meC$Col_CHHm[All_meC$hcluster==2],All_meC$Col_CHHm[All_meC$hcluster==3],All_meC$Col_CHHm[All_meC$hcluster==4],All_meC$Col_CHHm[All_meC$hcluster==5],
cex=0,col = c("blue","gray","forest green","purple","magenta"),ylim=c(0,1),ylab="methylation",main="Wild type")
boxplot(All_meC$ibm1_CGm[All_meC$hcluster==1],All_meC$ibm1_CGm[All_meC$hcluster==2],All_meC$ibm1_CGm[All_meC$hcluster==3],All_meC$ibm1_CGm[All_meC$hcluster==4],All_meC$ibm1_CGm[All_meC$hcluster==5],
All_meC$ibm1_CHGm[All_meC$hcluster==1],All_meC$ibm1_CHGm[All_meC$hcluster==2],All_meC$ibm1_CHGm[All_meC$hcluster==3],All_meC$ibm1_CHGm[All_meC$hcluster==4],All_meC$ibm1_CHGm[All_meC$hcluster==5],
All_meC$ibm1_CHHm[All_meC$hcluster==1],All_meC$ibm1_CHHm[All_meC$hcluster==2],All_meC$ibm1_CHHm[All_meC$hcluster==3],All_meC$ibm1_CHHm[All_meC$hcluster==4],All_meC$ibm1_CHHm[All_meC$hcluster==5],
cex=0,col = c("blue","gray","forest green","purple","magenta"),ylim=c(0,1),ylab="methylation",main="ibm1")
boxplot(All_meC$ibm1_ldl2_CGm[All_meC$hcluster==1],All_meC$ibm1_ldl2_CGm[All_meC$hcluster==2],All_meC$ibm1_ldl2_CGm[All_meC$hcluster==3],All_meC$ibm1_ldl2_CGm[All_meC$hcluster==4],All_meC$ibm1_ldl2_CGm[All_meC$hcluster==5],
All_meC$ibm1_ldl2_CHGm[All_meC$hcluster==1],All_meC$ibm1_ldl2_CHGm[All_meC$hcluster==2],All_meC$ibm1_ldl2_CHGm[All_meC$hcluster==3],All_meC$ibm1_ldl2_CHGm[All_meC$hcluster==4],All_meC$ibm1_ldl2_CHGm[All_meC$hcluster==5],
All_meC$ibm1_ldl2_CHHm[All_meC$hcluster==1],All_meC$ibm1_ldl2_CHHm[All_meC$hcluster==2],All_meC$ibm1_ldl2_CHHm[All_meC$hcluster==3],All_meC$ibm1_ldl2_CHHm[All_meC$hcluster==4],All_meC$ibm1_ldl2_CHHm[All_meC$hcluster==5],
cex=0,col = c("blue","gray","forest green","purple","magenta"),ylim=c(0,1),ylab="methylation",main="ibm1 ldl2")
#比較しやすくするためクラスターごとに並べた場合(Cコンテクストごとに3つのパネル)
boxplot(All_meC$Col_CGm[All_meC$hcluster==1],All_meC$ibm1_CGm[All_meC$hcluster==1],All_meC$ibm1_ldl2_CGm[All_meC$hcluster==1],All_meC$Col_CGm[All_meC$hcluster==2],All_meC$ibm1_CGm[All_meC$hcluster==2],All_meC$ibm1_ldl2_CGm[All_meC$hcluster==2],
All_meC$Col_CGm[All_meC$hcluster==3],All_meC$ibm1_CGm[All_meC$hcluster==3],All_meC$ibm1_ldl2_CGm[All_meC$hcluster==3],All_meC$Col_CGm[All_meC$hcluster==4],All_meC$ibm1_CGm[All_meC$hcluster==4],All_meC$ibm1_ldl2_CGm[All_meC$hcluster==4],
All_meC$Col_CGm[All_meC$hcluster==5],All_meC$ibm1_CGm[All_meC$hcluster==5],All_meC$ibm1_ldl2_CGm[All_meC$hcluster==5],
cex=0,col = c("blue","blue","blue","gray","gray","gray","forest green","forest green","forest green","purple","purple","purple","magenta","magenta","magenta"),ylim=c(0,1),ylab="methylation",main="mCG")
boxplot(All_meC$Col_CHGm[All_meC$hcluster==1],All_meC$ibm1_CHGm[All_meC$hcluster==1],All_meC$ibm1_ldl2_CHGm[All_meC$hcluster==1],All_meC$Col_CHGm[All_meC$hcluster==2],All_meC$ibm1_CHGm[All_meC$hcluster==2],All_meC$ibm1_ldl2_CHGm[All_meC$hcluster==2],
All_meC$Col_CHGm[All_meC$hcluster==3],All_meC$ibm1_CHGm[All_meC$hcluster==3],All_meC$ibm1_ldl2_CHGm[All_meC$hcluster==3],All_meC$Col_CHGm[All_meC$hcluster==4],All_meC$ibm1_CHGm[All_meC$hcluster==4],All_meC$ibm1_ldl2_CHGm[All_meC$hcluster==4],
All_meC$Col_CHGm[All_meC$hcluster==5],All_meC$ibm1_CHGm[All_meC$hcluster==5],All_meC$ibm1_ldl2_CHGm[All_meC$hcluster==5],
cex=0,col = c("blue","blue","blue","gray","gray","gray","forest green","forest green","forest green","purple","purple","purple","magenta","magenta","magenta"),ylim=c(0,1),ylab="methylation",main="mCHG")
boxplot(All_meC$Col_CHHm[All_meC$hcluster==1],All_meC$ibm1_CHHm[All_meC$hcluster==1],All_meC$ibm1_ldl2_CHHm[All_meC$hcluster==1],All_meC$Col_CHHm[All_meC$hcluster==2],All_meC$ibm1_CHHm[All_meC$hcluster==2],All_meC$ibm1_ldl2_CHHm[All_meC$hcluster==2],
All_meC$Col_CHHm[All_meC$hcluster==3],All_meC$ibm1_CHHm[All_meC$hcluster==3],All_meC$ibm1_ldl2_CHHm[All_meC$hcluster==3],All_meC$Col_CHHm[All_meC$hcluster==4],All_meC$ibm1_CHHm[All_meC$hcluster==4],All_meC$ibm1_ldl2_CHHm[All_meC$hcluster==4],
All_meC$Col_CHHm[All_meC$hcluster==5],All_meC$ibm1_CHHm[All_meC$hcluster==5],All_meC$ibm1_ldl2_CHHm[All_meC$hcluster==5],
cex=0,col = c("blue","blue","blue","gray","gray","gray","forest green","forest green","forest green","purple","purple","purple","magenta","magenta","magenta"),ylim=c(0,0.3),ylab="methylation",main="mCHH")
```
### ・Question
・*ibm1*変異体で上昇するCHG/CHHメチル化は*ibm1 ldl2*変異体ではどうなるか？
・*ibm1*変異体でCHG/CHHメチル化が上昇する遺伝子と上昇しない遺伝子ではCGメチル化(野生型および変異体で)はどうか？
## 3-3. WGBSのIGV用のファイル作成
unix commandの一つ、awkを使って、ゲノム上のシトシンの位置と各位置のシトシンのメチル化レベルを書き込んだbedGraphファイルを作製する。IGVに読み込んで、**実験で使用した遺伝子領域**や、解析から興味が持たれる領域についてパターンを見てみる。
・CGメチル化
```bash
echo "track type=bedGraph" > ~/Practice_2021/Col.CG.bedGraph && \
awk -v OFS="\t" '($4 + $5 > 0 && $6=="CG"){print $1,$2-1,$2,$4/($4+$5)*100}' \
ColRL_1_paired.fastq_bismark_pe.deduplicated.CX_report.txt >> ~/Practice_2021/Col.CG.bedGraph &
```
・CHGメチル化
```bash
echo "track type=bedGraph" > ~/Practice_2021/Col.CHG.bedGraph && \
awk -v OFS="\t" '($4 + $5 > 0 && $6=="CHG"){print $1,$2-1,$2,$4/($4+$5)*100}' \
ColRL_1_paired.fastq_bismark_pe.deduplicated.CX_report.txt >> ~/Practice_2021/Col.CHG.bedGraph &
```
・CHHメチル化(ファイルが重いのでスキップ)
```bash
echo "track type=bedGraph" > ~/Practice_2021/Col.CHH.bedGraph && \
awk -v OFS="\t" '($4 + $5 > 0 && $6=="CHH"){print $1,$2-1,$2,$4/($4+$5)*100}' \
ColRL_1_paired.fastq_bismark_pe.deduplicated.CX_report.txt >> ~/Practice_2021/Col.CHH.bedGraph &
```
[bedGraphファイルを圧縮したbigwigファイルはこちらからダウンロード](https://drive.google.com/drive/folders/1-zbE531Z7zruhobOdoCpIWu-2z3IeY63?usp=sharing)
**IGV(Integrated Genome Viewer)**
・様々な用途に使えるgenomeブラウザ。以下のウェブサイトからダウンロードして使う。
http://software.broadinstitute.org/software/igv/download
bigwigファイルをIGVで見てみる。ドラッグ&ドロップ、もしくはFile > Load from Fileからbigwigファイル(xxx.bw)を選択する。
*なぜか、染色体ごとにみようとすると、bwファイルがみられない！！来年までに解決しないと。
# ＜実験4. ChIP-seqとRNA-seqデータの解析＞
## 4-1. ChIP-seqによるヒストン修飾解析
今回の実習に使用したCol(野生型)、*ibm1*、*ibm1 ldl2*、におけるヒストン修飾のChIP-seqデータを解析する。[各サンプルについて、各遺伝子にマッピングされたリード数をリスト化したファイル](https://drive.google.com/drive/folders/1virN1_Damugs7t9CmeD6Grch4ifyR5-V?usp=sharing)を使って解析する。DNAメチル化データと同様、Import Datasetからデータフレームとして取り込む。（"Heading"のところでYesを選ぶことを忘れないように！）。まず、総リード数（今回は、すべての遺伝子にマップされたリード数の総和を使う）及び遺伝子長で補正した値であるRPKM値(Reads Per Kilobase per Million mapped reads)を計算し、サンプル間で比較する。メチル化やRNAの解析と同様に、RStudioで散布図などを書いてみる。DNAメチル化パターンとの相関関係を調べる。
＜以下、来年はbigwigに変更＞
また、気になる遺伝子領域のヒストン修飾パターンを調べるには、各サンプルのマッピングデータから作製した[tdfファイル](https://drive.google.com/drive/folders/1fRhZ34MpNw3xfcvfeZv3gX2ElAxux0-A?usp=sharing) をIGVで見てみる。tdfファイルをIGVで見る際は、ファイルをドラッグアンドドロップでIGVに入れるか、Fileメニューから、Load from File...で指定する。また、Viewメニューから、Preferences...を選び、"Tracks"タブの"Normalize coverage data (.tdf files only)"にチェックを入れてから、一度IGVを再起動すると、tdfファイルでChIP-seqの総リード数補正値が表示できる。
**リード数補正の例 遺伝子ごとのリード数を[列の数字の合計]/100万で割る**
＊このコマンドは下のコマンド群を実行するときに実際にやる必要はない。
```r=
WT_Col_ChIP$WT_H3_RPM <- WT_Col_ChIP$H3_reads*1000000/sum(WT_Col_ChIP$H3_reads)
```
DNAメチル化のクラスタリングで分けたクラスターごとのヒストン修飾パターンを見てみよう。
```r=
#リード数補正 RPM（reads per million reads）の列名に株の名前を入れておく
WT_Col_ChIP$WT_H3_RPM <- WT_Col_ChIP$H3_reads*1000000/sum(WT_Col_ChIP$H3_reads)
WT_Col_ChIP$WT_H3K9me2_RPM <- WT_Col_ChIP$H3K9me2_reads*1000000/sum(WT_Col_ChIP$H3K9me2_reads)
WT_Col_ChIP$WT_H3K4me2_RPM <- WT_Col_ChIP$H3K4me2_reads*1000000/sum(WT_Col_ChIP$H3K4me2_reads)
WT_Col_ChIP$WT_H3K4me1_RPM <- WT_Col_ChIP$H3K4me1_reads*1000000/sum(WT_Col_ChIP$H3K4me1_reads)
ibm1_ChIP$ibm1_H3_RPM <- ibm1_ChIP$H3_reads*1000000/sum(ibm1_ChIP$H3_reads)
ibm1_ChIP$ibm1_H3K9me2_RPM <- ibm1_ChIP$H3K9me2_reads*1000000/sum(ibm1_ChIP$H3K9me2_reads)
ibm1_ChIP$ibm1_H3K4me2_RPM <- ibm1_ChIP$H3K4me2_reads*1000000/sum(ibm1_ChIP$H3K4me2_reads)
ibm1_ChIP$ibm1_H3K4me1_RPM <- ibm1_ChIP$H3K4me1_reads*1000000/sum(ibm1_ChIP$H3K4me1_reads)
ibm1_ldl2_ChIP$ibm1_ldl2_H3_RPM <- ibm1_ldl2_ChIP$H3_reads*1000000/sum(ibm1_ldl2_ChIP$H3_reads)
ibm1_ldl2_ChIP$ibm1_ldl2_H3K9me2_RPM <- ibm1_ldl2_ChIP$H3K9me2_reads*1000000/sum(ibm1_ldl2_ChIP$H3K9me2_reads)
ibm1_ldl2_ChIP$ibm1_ldl2_H3K4me2_RPM <- ibm1_ldl2_ChIP$H3K4me2_reads*1000000/sum(ibm1_ldl2_ChIP$H3K4me2_reads)
ibm1_ldl2_ChIP$ibm1_ldl2_H3K4me1_RPM <- ibm1_ldl2_ChIP$H3K4me1_reads*1000000/sum(ibm1_ldl2_ChIP$H3K4me1_reads)
#遺伝子の長さ補正 RPKM (reads per kilobase per million reads)にする
WT_Col_ChIP$WT_H3_RPKM <- (WT_Col_ChIP$H3_reads*1000000/sum(WT_Col_ChIP$H3_reads))/((WT_Col_ChIP$End-WT_Col_ChIP$Start)/1000)
WT_Col_ChIP$WT_H3K9me2_RPKM <- (WT_Col_ChIP$H3K9me2_reads*1000000/sum(WT_Col_ChIP$H3K9me2_reads))/((WT_Col_ChIP$End-WT_Col_ChIP$Start)/1000)
WT_Col_ChIP$WT_H3K4me2_RPKM <- (WT_Col_ChIP$H3K4me2_reads*1000000/sum(WT_Col_ChIP$H3K4me2_reads))/((WT_Col_ChIP$End-WT_Col_ChIP$Start)/1000)
WT_Col_ChIP$WT_H3K4me1_RPKM <- (WT_Col_ChIP$H3K4me1_reads*1000000/sum(WT_Col_ChIP$H3K4me1_reads))/((WT_Col_ChIP$End-WT_Col_ChIP$Start)/1000)
ibm1_ChIP$ibm1_H3_RPKM <- (ibm1_ChIP$H3_reads*1000000/sum(ibm1_ChIP$H3_reads))/((WT_Col_ChIP$End-WT_Col_ChIP$Start)/1000)
ibm1_ChIP$ibm1_H3K9me2_RPKM <- (ibm1_ChIP$H3K9me2_reads*1000000/sum(ibm1_ChIP$H3K9me2_reads))/((WT_Col_ChIP$End-WT_Col_ChIP$Start)/1000)
ibm1_ChIP$ibm1_H3K4me2_RPKM <- (ibm1_ChIP$H3K4me2_reads*1000000/sum(ibm1_ChIP$H3K4me2_reads))/((WT_Col_ChIP$End-WT_Col_ChIP$Start)/1000)
ibm1_ChIP$ibm1_H3K4me1_RPKM <- (ibm1_ChIP$H3K4me1_reads*1000000/sum(ibm1_ChIP$H3K4me1_reads))/((WT_Col_ChIP$End-WT_Col_ChIP$Start)/1000)
ibm1_ldl2_ChIP$ibm1_ldl2_H3_RPKM <- (ibm1_ldl2_ChIP$H3_reads*1000000/sum(ibm1_ldl2_ChIP$H3_reads))/((WT_Col_ChIP$End-WT_Col_ChIP$Start)/1000)
ibm1_ldl2_ChIP$ibm1_ldl2_H3K9me2_RPKM <- (ibm1_ldl2_ChIP$H3K9me2_reads*1000000/sum(ibm1_ldl2_ChIP$H3K9me2_reads))/((WT_Col_ChIP$End-WT_Col_ChIP$Start)/1000)
ibm1_ldl2_ChIP$ibm1_ldl2_H3K4me2_RPKM <- (ibm1_ldl2_ChIP$H3K4me2_reads*1000000/sum(ibm1_ldl2_ChIP$H3K4me2_reads))/((WT_Col_ChIP$End-WT_Col_ChIP$Start)/1000)
ibm1_ldl2_ChIP$ibm1_ldl2_H3K4me1_RPKM <- (ibm1_ldl2_ChIP$H3K4me1_reads*1000000/sum(ibm1_ldl2_ChIP$H3K4me1_reads))/((WT_Col_ChIP$End-WT_Col_ChIP$Start)/1000)
#すべての株のChIPデータをまとめる
All_ChIP <- WT_Col_ChIP[,c(4,5,15,16,17,18)]
All_ChIP$ibm1_H3_RPKM <- ibm1_ChIP$ibm1_H3_RPKM
All_ChIP$ibm1_H3K9me2_RPKM <- ibm1_ChIP$ibm1_H3K9me2_RPKM
All_ChIP$ibm1_H3K4me2_RPKM <- ibm1_ChIP$ibm1_H3K4me2_RPKM
All_ChIP$ibm1_H3K4me1_RPKM <- ibm1_ChIP$ibm1_H3K4me1_RPKM
All_ChIP$ibm1_ldl2_H3_RPKM <- ibm1_ldl2_ChIP$ibm1_ldl2_H3_RPKM
All_ChIP$ibm1_ldl2_H3K9me2_RPKM <- ibm1_ldl2_ChIP$ibm1_ldl2_H3K9me2_RPKM
All_ChIP$ibm1_ldl2_H3K4me2_RPKM <- ibm1_ldl2_ChIP$ibm1_ldl2_H3K4me2_RPKM
All_ChIP$ibm1_ldl2_H3K4me1_RPKM <- ibm1_ldl2_ChIP$ibm1_ldl2_H3K4me1_RPKM
#DNAメチル化データとChIPデータをまとめる (merge機能を使って、All_meCの"gene"列、All_ChIPの"ID"列を基準に)
All_meC_ChIP<- merge(All_meC,All_ChIP,by.x = "gene",by.y = "ID")
#ChIP RPKMの散布図（k-meansクラスター(All_meC_ChIP$kmclusterの列の数字)で色分けして）H3の場合 (ドットの透明度は各クラスターの遺伝子数に応じて調整するとよい（rgb()の4番目の数字を小さくすると透明度が上がる）)
#WT - ibm1 比較
plot(All_meC_ChIP$WT_H3_RPKM,All_meC_ChIP$ibm1_H3_RPKM,pch=16,cex=0.3,xlim=c(0,40),ylim=c(0,40),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP$kmcluster)],xlab="H3 WT",ylab="H3 ibm1")
#WT - ibm1 ldl2 比較
plot(All_meC_ChIP$WT_H3_RPKM,All_meC_ChIP$ibm1_ldl2_H3_RPKM,pch=16,cex=0.3,xlim=c(0,40),ylim=c(0,40),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP$kmcluster)],xlab="H3 WT",ylab="H3 ibm1 ldl2")
#修飾でも同様の散布図を書いてみよう。x軸、y軸の範囲は適宜、xlim=c(0,40),ylim=c(0,40)の部分を変更して見やすいように変える。
```
クラスターごとに、野生型と変異体間の修飾レベル変化（*ibm1* - WT or *ibm1 ldl2* - WT）の散布図を描いてみる
```r=
#RPKMの野生型変異体間の変化（引き算）を計算
All_meC_ChIP$ibm1_H3K4me1_diff <- All_meC_ChIP$ibm1_H3K4me1_RPKM-All_meC_ChIP$WT_H3K4me1_RPKM
All_meC_ChIP$ibm1_ldl2_H3K4me1_diff <- All_meC_ChIP$ibm1_ldl2_H3K4me1_RPKM-All_meC_ChIP$WT_H3K4me1_RPKM
All_meC_ChIP$ibm1_H3K4me2_diff <- All_meC_ChIP$ibm1_H3K4me2_RPKM-All_meC_ChIP$WT_H3K4me2_RPKM
All_meC_ChIP$ibm1_ldl2_H3K4me2_diff <- All_meC_ChIP$ibm1_ldl2_H3K4me2_RPKM-All_meC_ChIP$WT_H3K4me2_RPKM
All_meC_ChIP$ibm1_H3K9me2_diff <- All_meC_ChIP$ibm1_H3K9me2_RPKM-All_meC_ChIP$WT_H3K9me2_RPKM
All_meC_ChIP$ibm1_ldl2_H3K9me2_diff <- All_meC_ChIP$ibm1_ldl2_H3K9me2_RPKM-All_meC_ChIP$WT_H3K9me2_RPKM
#クラスターで色分けをした修飾変化の散布図 (ドットの透明度は各クラスターの遺伝子数に応じて調整するとよい)
#WT vs ibm1でのH3K9me2とH3K4me1の比較
plot(All_meC_ChIP$ibm1_H3K9me2_diff,All_meC_ChIP$ibm1_H3K4me1_diff,pch=16,cex=0.3,xlim=c(-20,40),ylim=c(-30,20),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP$kmcluster)],xlab="H3K9me2 change (ibm1 - WT)",ylab="H3K4me1 change (ibm1 - WT)")
#WT vs ibm1 ldl2でのH3K9me2とH3K4me1の比較
plot(All_meC_ChIP$ibm1_ldl2_H3K9me2_diff,All_meC_ChIP$ibm1_ldl2_H3K4me1_diff,pch=16,cex=0.3,xlim=c(-20,40),ylim=c(-30,20),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP$kmcluster)],xlab="H3K9me2 change (ibm1 ldl2 - WT)",ylab="H3K4me1 change (ibm1 ldl2 - WT)")
#H3K4me1について、{WT vs ibm1} と {WT vs ibm1 ldl2}の比較
plot(All_meC_ChIP$ibm1_H3K4me1_diff,All_meC_ChIP$ibm1_ldl2_H3K4me1_diff,pch=16,cex=0.3,xlim=c(-30,20),ylim=c(-30,20),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP$kmcluster)],xlab="H3K4me1 change (ibm1 - WT)",ylab="H3K4me1 change (ibm1 ldl2 - WT)")
#H3K4me2でも同様に
plot(All_meC_ChIP$ibm1_H3K9me2_diff,All_meC_ChIP$ibm1_H3K4me2_diff,pch=16,cex=0.3,xlim=c(-20,40),ylim=c(-30,20),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP$kmcluster)],xlab="H3K9me2 change (ibm1 - WT)",ylab="H3K4me2 change (ibm1 - WT)")
plot(All_meC_ChIP$ibm1_ldl2_H3K9me2_diff,All_meC_ChIP$ibm1_ldl2_H3K4me2_diff,pch=16,cex=0.3,xlim=c(-20,40),ylim=c(-30,20),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP$kmcluster)],xlab="H3K9me2 change (ibm1 ldl2 - WT)",ylab="H3K4me2 change (ibm1 ldl2 - WT)")
plot(All_meC_ChIP$ibm1_H3K4me2_diff,All_meC_ChIP$ibm1_ldl2_H3K4me2_diff,pch=16,cex=0.3,xlim=c(-30,20),ylim=c(-30,20),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP$kmcluster)],xlab="H3K4me2 change (ibm1 - WT)",ylab="H3K4me2 change (ibm1 ldl2 - WT)")
```
### ・Question
*ibm1*変異体でH3K9me2が上昇する遺伝子について、*ibm1 ldl2*二重変異体ではH3K9me2変化がみられるか？
*ibm1*変異体でH3K4me1/H3K4me2の変化は？H3K9me2/mCHG変化との関連は？
*ibm1 ldl2*二重変異体におけるH3K4me1/H3K4me2は？
## 4-2. トランスクリプトーム解析
今回の実習に使用したCol(野生型)、*ibm1*、*ibm1 ldl2*、のRNA-seqデータを解析する。各サンプルについて3回 (独立の3実験; biological triplicates)のRNA-seqを実施し、各遺伝子にマッピングされたリード数をまとめた（[こちらからダウンロード](https://drive.google.com/drive/folders/1hku-X9O1MZEjfjEQgHPZCP0QDIbdeoLn?usp=sharing )）。
[参考ページ(東京大学大学院情報学環・学際情報学府 門田先生による「Rでシーケンス解析」)](http://www.iu.a.u-tokyo.ac.jp/~kadota/r_seq.html)
### 4-2-1.発現変動遺伝子の抽出
それぞれの2サンプルの比較（WT vs *ibm1*, WT vs *ibm1 ldl2*, *ibm1* vs *ibm1 ldl2*）で発現量が変化している遺伝子を抽出してみる。[出典](http://www.iu.a.u-tokyo.ac.jp/~kadota/r_seq.html#about_analysis_deg_2_paired)
1.まず最初に、edgeRパッケージをインストール@R studio
*M1 Macの場合は[こちら](https://mac.r-project.org/tools/)を参考に、GNU Fortrun Compiler (Apple Silicon Macs)を入れる必要がある。(もう必要ない？@2024)
```r=
#方法1
install.packages('edgeR')
#方法2
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(version = "3.14")
library("BiocManager")
install('edgeR')
#方法3 (R version 4.1以降はこちらがベター？)
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("edgeR")
```
2.以下のスクリプトをR studioのR scriptに貼り付け、1行ずつ実行。入力ファイル名(in_f)と出力ファイル名(out_f)は適宜変更。
```r=
in_f <- "WT_ibm1_RNAseq.txt" #入力ファイル名を指定してin_fに格納
out_f <- "WT_ibm1_DEG.txt" #出力ファイル名を指定してout_fに格納
param_G1 <- 3 #G1群のサンプル数を指定
param_G2 <- 3 #G2群のサンプル数を指定
param_FDR <- 0.05 #false discovery rate (FDR)閾値を指定
#必要なパッケージをロード
library(edgeR) #パッケージの読み込み
#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
#前処理(1. DGEListオブジェクトの作成、2. 実験デザイン行列の作成など)
d <- DGEList(counts=data) #1. DGEListオブジェクトdを作成
pair <- as.factor(c(1:param_G1, 1:param_G2))#2. 対応関係の情報をpairに格納
cl <- as.factor(c(rep(1, param_G1), rep(2, param_G2)))#2. グループラベル情報をclに格納
design <- model.matrix(~ pair + cl) #2. デザイン行列を作成した結果をdesignに格納
FDR <- 0.1 #DEGES正規化内部でのDEG検出時のFDR閾値を指定
floorPDEG <- 0.05 #上記FDR閾値を満たす遺伝子数がここで指定した値(5%)未満だったときに、最低でも上位5%の遺伝子をDEGとみなして除去するための下限値を指定
#本番(DEGES/edgeR正規化)
### STEP 1 ###
d <- calcNormFactors(d) #TMM正規化を実行した結果をdに格納
d$samples$norm.factors #確認してるだけです
### STEP 2 ###
d <- estimateGLMCommonDisp(d, design) #モデル構築(ばらつきの程度を見積もっている)
d <- estimateGLMTrendedDisp(d, design) #モデル構築(ばらつきの程度を見積もっている)
d <- estimateGLMTagwiseDisp(d, design) #モデル構築(ばらつきの程度を見積もっている)
fit <- glmFit(d, design) #NB GLMへのfitting
out <- glmLRT(fit) #検定
p.value <- out$table$PValue #p値をp.valueに格納
q.value <- p.adjust(p.value, method="BH")#q値をq.valueに格納
if (sum(q.value < FDR) > (floorPDEG * nrow(data))){#FDR閾値を満たす遺伝子数がここで指定した値(5%)よりも多いかどうか?
is.DEG <- as.logical(q.value < FDR) #Yesの場合は、FDR閾値を満たす遺伝子がTRUE、そうでないものがFALSEと判定した結果をis.DEGベクトルに格納
} else { #Noの場合は...
is.DEG <- as.logical(rank(p.value, ties.method="min") <= nrow(data)*floorPDEG)#上位(floorPDEG)%がTRUE、そうでないものがFALSEと判定した結果をis.DEGベクトルに格納
}
### STEP 3 ###
d <- DGEList(counts=data[!is.DEG, ]) #DGEListオブジェクトdを作成(non-DEG)
d <- calcNormFactors(d) #TMM正規化を実行した結果をdに格納
norm.factors <- d$samples$norm.factors*colSums(data[!is.DEG, ])/colSums(data)#non-DEGのみ用に得られたTMM正規化係数を全遺伝子に掛けられるようにするための補正
norm.factors <- norm.factors/mean(norm.factors)#TMM正規化を実行した結果をdに格納
norm.factors #確認してるだけです
#本番(DEG検出)
d <- DGEList(counts=data) #DGEListオブジェクトdを作成
d$samples$norm.factors <- norm.factors #計算した正規化係数を代入
d <- estimateGLMCommonDisp(d, design) #モデル構築(ばらつきの程度を見積もっている)
d <- estimateGLMTrendedDisp(d, design) #モデル構築(ばらつきの程度を見積もっている)
d <- estimateGLMTagwiseDisp(d, design) #モデル構築(ばらつきの程度を見積もっている)
fit <- glmFit(d, design) #NB GLMへのfitting
out <- glmLRT(fit) #検定
#tmp <- topTags(out, n=nrow(data), sort.by="none")#検定結果を抽出
p.value <- out$table$PValue #p値をp.valueに格納
q.value <- p.adjust(p.value, method="BH")#q値をq.valueに格納
ranking <- rank(p.value) #p.valueでランキングした結果をrankingに格納
sum(q.value < param_FDR) #FDR閾値(q.value < param_FDR)を満たす遺伝子数を表示
#ファイルに保存(テキストファイル)
tmp <- cbind(rownames(data), data, p.value, q.value, ranking)#入力データの右側にp.value、q.value、rankingを結合した結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
```
得られたアウトプットファイルをRやエクセルで見てみる。
```r=
tmp[,11]<-apply(tmp[,2:4],1,mean);colnames(tmp)[11]<-'average_of_G1' #G1群の平均
tmp[,12]<-apply(tmp[,5:7],1,mean);colnames(tmp)[12]<-'average_of_G2' #G2群の平均
DEGS<-tmp[tmp$q.value<param_FDR,];nrow(DEGS) #有意に発現変動している遺伝子を含む列を取り出して'DEGS'に格納
DEGS_down<-DEGS[DEGS$average_of_G1>DEGS$average_of_G2,1] #DEGSのうち,G1>G2の遺伝子リストをDEGS_downに格納
DEGS_up<-DEGS[DEGS$average_of_G1<DEGS$average_of_G2,1] #DEGSのうち,G1<G2の遺伝子をDEGS_upに格納
#発現上昇・低下遺伝子をテキストファイルに書き出す。
write.table(DEGS_down,"ibm1_WT_DEGS_down.txt",quote=F,row.names=F)
write.table(DEGS_up,"ibm1_WT_DEGS_up.txt",quote=F,row.names=F)
```
発現上昇遺伝子、低下遺伝子を抽出し、[GO解析](http://geneontology.org/)をやってみる.
*GO（Gene Ontology）解析とは？ 生物学的機能や生化学的機能などによって分類分けされた遺伝子産物が任意の遺伝子集団に濃縮しているかを調べる統計的手法
### ・Question
*ibm1*変異体で発現が上昇・低下する遺伝子はそれぞれどんな遺伝子が多いか？
*ibm1*変異体でH3K9me2やH3K4me1/2が変化する遺伝子(クラスター)では発現量は変動するか？
### 4-2-2. 遺伝子発現とDNAメチル化・ヒストンメチル化の関係性の解析
まず、RNAseq_all.txtをimport datasetから取り込む。HeadingをYesにするのを忘れないように。
```r=
#TPM (Transcript per million reads)の計算
RNAseq_all$WT_1_TPM <- ((RNAseq_all$WT_1/RNAseq_all$cDNA_length*1000)/sum(RNAseq_all$WT_1/RNAseq_all$cDNA_length*1000))*1000000
RNAseq_all$WT_2_TPM <- ((RNAseq_all$WT_2/RNAseq_all$cDNA_length*1000)/sum(RNAseq_all$WT_2/RNAseq_all$cDNA_length*1000))*1000000
RNAseq_all$WT_3_TPM <- ((RNAseq_all$WT_3/RNAseq_all$cDNA_length*1000)/sum(RNAseq_all$WT_3/RNAseq_all$cDNA_length*1000))*1000000
RNAseq_all$ibm1_1_TPM <- ((RNAseq_all$ibm1_1/RNAseq_all$cDNA_length*1000)/sum(RNAseq_all$ibm1_1/RNAseq_all$cDNA_length*1000))*1000000
RNAseq_all$ibm1_2_TPM <- ((RNAseq_all$ibm1_2/RNAseq_all$cDNA_length*1000)/sum(RNAseq_all$ibm1_2/RNAseq_all$cDNA_length*1000))*1000000
RNAseq_all$ibm1_3_TPM <- ((RNAseq_all$ibm1_3/RNAseq_all$cDNA_length*1000)/sum(RNAseq_all$ibm1_3/RNAseq_all$cDNA_length*1000))*1000000
RNAseq_all$ibm1_ldl2_1_TPM <- ((RNAseq_all$ibm1_ldl2_1/RNAseq_all$cDNA_length*1000)/sum(RNAseq_all$ibm1_ldl2_1/RNAseq_all$cDNA_length*1000))*1000000
RNAseq_all$ibm1_ldl2_2_TPM <- ((RNAseq_all$ibm1_ldl2_2/RNAseq_all$cDNA_length*1000)/sum(RNAseq_all$ibm1_ldl2_2/RNAseq_all$cDNA_length*1000))*1000000
RNAseq_all$ibm1_ldl2_3_TPM <- ((RNAseq_all$ibm1_ldl2_3/RNAseq_all$cDNA_length*1000)/sum(RNAseq_all$ibm1_ldl2_3/RNAseq_all$cDNA_length*1000))*1000000
#3回の実験のTPMのlogをとって平均する
RNAseq_all$WT_log_mean <- (log(RNAseq_all$WT_1_TPM+0.005,2)+log(RNAseq_all$WT_2_TPM+0.005,2)+log(RNAseq_all$WT_2_TPM+0.005,2))/3
RNAseq_all$ibm1_log_mean <- (log(RNAseq_all$ibm1_1_TPM+0.005,2)+log(RNAseq_all$ibm1_2_TPM+0.005,2)+log(RNAseq_all$ibm1_2_TPM+0.005,2))/3
RNAseq_all$ibm1_ldl2_log_mean <- (log(RNAseq_all$ibm1_ldl2_1_TPM+0.005,2)+log(RNAseq_all$ibm1_ldl2_2_TPM+0.005,2)+log(RNAseq_all$ibm1_ldl2_2_TPM+0.005,2))/3
#BS-seq, ChIP-seqのデータとまとめる
All_meC_ChIP_RNA<- merge(All_meC_ChIP,RNAseq_all,by.x = "gene",by.y = "Locus")
#サンプル間の比較を散布図にしてみる
#WT vs ibm1
plot(All_meC_ChIP_RNA$WT_log_mean,All_meC_ChIP_RNA$ibm1_log_mean,pch=16,cex=0.2,col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP_RNA$kmcluster)])
#WT vs ibm1 ldl2
plot(All_meC_ChIP_RNA$WT_log_mean,All_meC_ChIP_RNA$ibm1_ldl2_log_mean,pch=16,cex=0.2,col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP_RNA$kmcluster)])
#ibm1 vs ibm1 ldl2
plot(All_meC_ChIP_RNA$ibm1_log_mean,All_meC_ChIP_RNA$ibm1_ldl2_log_mean,pch=16,cex=0.2,col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP_RNA$kmcluster)])
#短い遺伝子はノイズになりやすいので除く (<500bp)
All_meC_ChIP_RNA_500 <- subset(All_meC_ChIP_RNA,length>500)
#mCHGの変化とRNAの変化の関係
#WT vs ibm1
plot(All_meC_ChIP_RNA_500$ibm1_CHGm-All_meC_ChIP_RNA_500$Col_CHGm,All_meC_ChIP_RNA_500$ibm1_log_mean-All_meC_ChIP_RNA_500$WT_log_mean,pch=16,cex=0.2,xlim=c(-0.2,0.6),ylim=c(-6,6),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP_RNA_500$kmcluster)])
#WT vs ibm1 ldl2
plot(All_meC_ChIP_RNA_500$ibm1_ldl2_CHGm-All_meC_ChIP_RNA_500$Col_CHGm,All_meC_ChIP_RNA_500$ibm1_ldl2_log_mean-All_meC_ChIP_RNA_500$WT_log_mean,pch=16,cex=0.2,xlim=c(-0.2,0.6),ylim=c(-6,6),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP_RNA_500$kmcluster)])
#H3K9me2の変化とRNAの変化の関係
#WT vs ibm1
plot(All_meC_ChIP_RNA_500$ibm1_H3K9me2_diff,All_meC_ChIP_RNA_500$ibm1_log_mean-All_meC_ChIP_RNA_500$WT_log_mean,pch=16,cex=0.2,xlim=c(-25,50),ylim=c(-6,6),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP_RNA_500$kmcluster)])
#WT vs ibm1 ldl2
plot(All_meC_ChIP_RNA_500$ibm1_ldl2_H3K9me2_diff,All_meC_ChIP_RNA_500$ibm1_ldl2_log_mean-All_meC_ChIP_RNA_500$WT_log_mean,pch=16,cex=0.2,xlim=c(-25,50),ylim=c(-6,6),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP_RNA_500$kmcluster)])
#H3K4me1の変化とRNAの変化の関係
#WT vs ibm1
plot(All_meC_ChIP_RNA_500$ibm1_H3K4me1_diff,All_meC_ChIP_RNA_500$ibm1_log_mean-All_meC_ChIP_RNA_500$WT_log_mean,pch=16,cex=0.2,xlim=c(-25,25),ylim=c(-6,6),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP_RNA_500$kmcluster)])
#WT vs ibm1 ldl2
plot(All_meC_ChIP_RNA_500$ibm1_ldl2_H3K4me1_diff,All_meC_ChIP_RNA_500$ibm1_ldl2_log_mean-All_meC_ChIP_RNA_500$WT_log_mean,pch=16,cex=0.2,xlim=c(-25,25),ylim=c(-6,6),col=c(rgb(0,0,1,0.5), rgb(0,0,0,0.5), rgb(0,0.5,0,0.5), rgb(0.5,0,1,0.5), rgb(1,0,1,0.5))[unclass(All_meC_ChIP_RNA_500$kmcluster)])
#ibm1でmCHGが上がるクラスターだけ抜き出して散布図を書いてみよう(?にクラスター番号入れる)
plot(All_meC_ChIP_RNA_500$ibm1_H3K4me1_diff[All_meC_ChIP_RNA_500$kmcluster==?],All_meC_ChIP_RNA_500$ibm1_log_mean[All_meC_ChIP_RNA_500$kmcluster==?]-All_meC_ChIP_RNA_500$WT_log_mean[All_meC_ChIP_RNA_500$kmcluster==?],pch=16,cex=0.4,xlim=c(-35,10),ylim=c(-6,6),col=rgb(0,0,0,0.2))
plot(All_meC_ChIP_RNA_500$ibm1_ldl2_H3K4me1_diff[All_meC_ChIP_RNA_500$kmcluster==?],All_meC_ChIP_RNA_500$ibm1_ldl2_log_mean[All_meC_ChIP_RNA_500$kmcluster==?]-All_meC_ChIP_RNA_500$WT_log_mean[All_meC_ChIP_RNA_500$kmcluster==?],pch=16,cex=0.4,xlim=c(-35,10),ylim=c(-6,6),col=rgb(0,0,0,0.2))
```
optional: [多変量解析(PCA)](https://hackmd.io/@KakutaniLab/Sk9J_g5PF)もやってみよう
## 便利機能
### 遺伝子リストの特徴づけツール
* シロイヌナズナの国際データベース[Thalemine](https://bar.utoronto.ca/thalemine/begin.do)
真ん中のAnalyzeの窓に遺伝子リストをコピペしてANALYZEボタンを押す。次のページで遺伝子リストに名前をつけて(そのままでも可)、緑のsave the listボタンを押す-> そのリストの遺伝子に頻出する機能、ドメインなどがわかる
* シロイヌナズナの国際データベース2 [tair10](https://www.arabidopsis.org）
[ここ](https://www.arabidopsis.org/jsp/ChromosomeMap/tool.jsp)に遺伝子リストを入れると、染色体上の位置を図示してくれる
# ＜参考文献＞
[Saze et al 2008](https://www.science.org/doi/10.1126/science.1150987?url_ver=Z39.88-2003&rfr_id=ori:rid:crossref.org&rfr_dat=cr_pub%20%200pubmed)
[Inagaki et al 2010](https://www.embopress.org/doi/full/10.1038/emboj.2010.227)
[Inagaki et al 2017](https://www.embopress.org/doi/full/10.15252/embj.201694983)
[Inagaki 2021 (総説)](https://www.jstage.jst.go.jp/article/ggs/96/5/96_21-00041/_pdf/-char/en)
# ＜レポート課題＞
[こちらのフォームから12月中に提出してください](https://forms.gle/74nqhU4GUY7EiUG86)
Sign in with Wallet
Connect another wallet