---
# System prepended metadata

title: 21： Jing-ru EM-seq 解析 (Lib1〜8)
tags: [解析, JL, EM-seq, ノート]

---

###### tags:`ノート` `JL` `解析` `EM-seq` 

21： Jing-ru EM-seq 解析 (Lib1〜8)
===

<br>

#### 【 作業場所＆保存場所 】
 //osc-fs_home/s-maeda/21_JL_EM-seq_20220316


#### 【rawdata】
//analysisdata/rawseq/fastq/SHARED/000037/21_JL_EM-seq_20220316_rawdata


#### 【参考】

bismarkでのメチル化解析 HackMD
https://hackmd.io/@takahirosuzuki/B1VFrdwoH

<br>

---
#### 2022/03/16 

### 1. Adapter Trimming

今回の配列

Read 1: AGATCGGAAGAGCACACGTCTGAACTCCAGTCA
Read 2: AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT


#### 1_cutadapt.sh
```
#!/bin/bash

output1=$1
output2=$2
read1=$3
read2=$4

cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT -o $output1 -p $output2 $read1 $read2
```




<br>

### 2. クオリティコントロール

#### 2_prinseqpp.sh
```
#!/bin/bash

read1=$1
read2=$2
output=$3

conda activate prinseq-plus-plus
prinseq++ -fastq $read1 -fastq2 $read2 -out_name $output -trim_tail_left 5 -trim_tail_right 5 -trim_qual_left 30 -trim_qual_right 30 -ns_max_n 20 -min_len 50 -threads 8
```


<br>

### 3. FastQC

#### 3_fastqc.sh
```
#!/bin/bash

read=$1
fastqc --nogroup -t 8 -o FastQC_reports $read
```

FastQCの結果は
//osc-fs_home/s-maeda/21_JL_EM-seq_20220316/FastQC_reports
にある。

<br>

### 4. Bismark マッピング

#### 4_bismark.sh
```
#!/bin/bash

read1=$1
read2=$2

bismark --bowtie2 --fastq --parallel 8 -N 1 bisulfite -1 $read1 -2 $read2
```


Lib1 Mapping efficiency: 70.8%
Lib2 Mapping efficiency: 53.0%
Lib3 Mapping efficiency: 72.4%
Lib4 Mapping efficiency: 74.0%
Lib5 Mapping efficiency: 76.9%
Lib6 Mapping efficiency: 76.9%
Lib7 Mapping efficiency: 76.8%
Lib8 Mapping efficiency: 76.4%

Lib2のマッピング効率が低い。
アダプター少し変えて再度やってみたが変わらず...

<br>

### 5. duplicateの除去
```
#!/bin/bash

bamfile=$1
deduplicate_bismark --bam -p $bamfile
```



<br>

### 6. Methylation call
```
#!/bin/bash

bamfile=$1
bismark_methylation_extractor --gzip --bedGraph --parallel 4  --buffer_size 10G --cytosine_report --genome_folder //osc-fs_home/s-maeda/21_JL_EM-seq_20220316/bisulfite $bamfile
```


<br>

### 7. IGVで可視化用の材料つくる
ソートしてindexを作る
```
#!/bin/bash
#SBATCH -p amd
#SBATCH -n 4

bamfile=$1
outputbam=$2

samtools sort $bamfile -o $outputbam
samtools index $outputbam
```



*.bedGraph, ソートしindexを作ったbamファイルをIGVなどのgenome browserで読み込むことで可視化できる。

<br>

### 8. MethylKitでそれぞれ比較

#### Sample情報

| Number | Name |
| ---- | ---- |
| Lib1 | WT_Day6　　|
| Lib2 | WT_Day9 |
| Lib3 | KO_Day6 |
| Lib4 | KO_Day9 |
| Lib5 | VP64_Day6   |
| Lib6 | VP64_Day6   |
| Lib7 | WT_Day9_P6  |
| Lib8 | VP64_Day9_P6 |


#### 組み合わせ

* WT_Day9　vs　WT_Day6
* KO_Day9　vs　KO_Day6
* VP64_Day9　vs　VP64_Day6
* WT_Day9_P6　vs　WT_Day6
* VP64_Day9_P6　vs　VP64_Day6

<br>

### bam fileの読み込み
 
8_methylKit_read2.sh

```
#!/bin/bash
#SBATCH -p amd
#SBATCH -n 4

source ~/.bash_profile
R --vanilla --slave --args $1 $2 $3 $4 $5 < 8_methylKit_read2.R
```

8_methylKit_read2.R
```
args <- commandArgs(trailingOnly = T)

library(methylKit)

file1 <- args[1]
file2 <- args[2]
sample1 <- args[3]
sample2 <- args[4]
rdata <- args[5]


files <- list(file1, file2)
sample.id <- list(sample1, sample2)
group_vct <- c(0,1)

my.methRaw.list=processBismarkAln( location = files, sample.id=sample.id, assembly="hg38", read.context="CpG", save.folder=getwd(), mincov = 4, treatment=group_vct)

save(list=ls(), file=rdata)
```

