---
# System prepended metadata

title: 30-01： Qin CUT&RUN 解析
tags: [解析, ノート, Cut&Run]

---

###### tags: `ノート` `解析` `Cut&Run`

# 30-01： Qin CUT&RUN 解析 

<br>

#### 2025/1/20

#### 【 作業場所＆保存場所 】
 //osc-fs_home/s-maeda/30_Qin/Analysis


#### 【rawdata】
//analysisdata/rawseq/fastq/SHARED/000039/30_Qin_MYCN_CUTandRUN/rawdata


#### 【参考】

[0から始めるエピゲノム解析（ChIP-seq）](https://github.com/yuifu/ngsdat2_epigenome_chipseq/blob/master/chipseq.md)
[NGSハンズオン講習会2015](https://github.com/suimye/NGS_handson2015/wiki/NGS%E3%83%87%E3%83%BC%E3%82%BF%E3%81%AE%E6%BA%96%E5%82%99%E3%81%A8mapping)
[平成28年度NGSハンズオン講習会 ChIP-seq](https://www.iu.a.u-tokyo.ac.jp/~kadota/bioinfo_ngs_sokushu_2016/2/20160728_amelieff_20160803.pdf)
[CUT&RUN data analysis](https://hbctraining.github.io/Intro-to-ChIPseq-flipped/CUT&RUN/Intro_to_CUT&RUN.html)
[CUT&RUN Processing Pipeline](https://data.4dnucleome.org/resources/data-analysis/cut-and-run-pipeline)


<br>

---

### 1. fastpでトリミング

fastp.sh
```
#!/bin/bash

#SBATCH --job-name=fastp          
#SBATCH --output=fastp_%A_%a.out         
#SBATCH --error=fastp_%A_%a.err          
#SBATCH --partition=batch                 
#SBATCH --cpus-per-task=8                 
#SBATCH --mem=8G                          
#SBATCH --array=1-12                     

module load fastp

# サンプルリスト
samples=("1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12")
sample_id=${samples[$SLURM_ARRAY_TASK_ID-1]}

# 入力ファイルと出力ディレクトリの設定
input_r1="${sample_id}_L1_1.fq.gz"
input_r2="${sample_id}_L1_2.fq.gz"
output_dir="trimmed_fastq"
mkdir -p ${output_dir}

# fastpコマンド実行
fastp \
    -i ${input_r1} \
    -I ${input_r2} \
    -o ${output_dir}/${sample_id}_L1_1.trimmed.fq.gz \
    -O ${output_dir}/${sample_id}_L1_2.trimmed.fq.gz \
    --detect_adapter_for_pe \
    --trim_poly_g \
    --length_required 30 \
    --thread ${SLURM_CPUS_PER_TASK}

echo "Trimming for sample ${sample_id} completed."

```

run script
```
sbatch qc_trimmed_files.sh
```  
<br/>

*⭐︎ fastpオプションの意味*
 -i ${input_r1}  ： read1 input
 -I ${input_r2}  ： read2 input 
-o ${output_dir}/${sample_id}_L1_1.trimmed.fq.gz   ：  read1 output
-O ${output_dir}/${sample_id}_L1_2.trimmed.fq.gz   ： read2 output
--detect_adapter_for_pe   ： ペアドエンド指定
--trim_poly_g   ： トリミング
--length_required 30   ： 30bpより短いものは破棄
--thread ${SLURM_CPUS_PER_TASK}　  ： スレッド数（今回は8）

<br>

---

### 2. FastQCでクオリティチェック

fastqc.sh
```
#!/bin/bash
#SBATCH --job-name=qc_analysis 
#SBATCH --output=qc_%A_%a.out
#SBATCH --error=qc_%A_%a.err
#SBATCH --partition=batch 
#SBATCH --cpus-per-task=8
#SBATCH --mem=8G 
#SBATCH --array=1-10 


module load fastqc
module load multiqc

# サンプルリスト
files=(
    "10_L1_1.trimmed.fq.gz" "10_L1_2.trimmed.fq.gz"
    "11_L1_1.trimmed.fq.gz" "11_L1_2.trimmed.fq.gz"
    "12_L1_1.trimmed.fq.gz" "12_L1_2.trimmed.fq.gz"
    "1_L1_1.trimmed.fq.gz"  "1_L1_2.trimmed.fq.gz"
    "2_L1_1.trimmed.fq.gz"  "2_L1_2.trimmed.fq.gz"
    "3_L1_1.trimmed.fq.gz"  "3_L1_2.trimmed.fq.gz"
    "4_L1_1.trimmed.fq.gz"  "4_L1_2.trimmed.fq.gz"
    "5_L1_1.trimmed.fq.gz"  "5_L1_2.trimmed.fq.gz"
    "6_L1_1.trimmed.fq.gz"  "6_L1_2.trimmed.fq.gz"
    "7_L1_1.trimmed.fq.gz"  "7_L1_2.trimmed.fq.gz"
    "8_L1_1.trimmed.fq.gz"  "8_L1_2.trimmed.fq.gz"
    "9_L1_1.trimmed.fq.gz"  "9_L1_2.trimmed.fq.gz"
)
input_file=${files[$SLURM_ARRAY_TASK_ID-1]}

# 出力ディレクトリを作成
output_dir="FastQC_reports"
mkdir -p ${output_dir}

# fastqcの実行
fastqc -o ${output_dir} -t ${SLURM_CPUS_PER_TASK} ${input_file}

# 配列ジョブの最終タスクでmultiqcを実行
if [ ${SLURM_ARRAY_TASK_ID} -eq ${SLURM_ARRAY_TASK_MAX} ]; then
    echo "Running multiqc..."
    multiqc ${output_dir} -o ${output_dir}
fi

echo "QC for ${input_file} completed."

```

実行
```
sbatch fastqc.sh
```

<br>

---

### 3. Bowtie2でマッピング
#### Human

bowtie2_hs.sh
```
#!/bin/bash
#SBATCH --job-name=bowtie2_hs
#SBATCH --output=bowtie2_hs%A_%a.out
#SBATCH --error=bowtie2_hs%A_%a.err 
#SBATCH --partition=batch  
#SBATCH --cpus-per-task=8 
#SBATCH --mem=16G 
#SBATCH --array=1-4

module load bowtie2

# インデックス
bowtie2_index="~/osc-fs/REFERENCES/bowtie2/hg38/hg38"

# サンプルリスト
samples=("12_L1" "2_L1" "11_L1" "1_L1")
sample=${samples[$SLURM_ARRAY_TASK_ID-1]}

# 入力ファイルの指定
read1="${sample}_1.trimmed.fq.gz"
read2="${sample}_2.trimmed.fq.gz"

# 出力ディレクトリの作成
output_dir="bowtie2_results"
mkdir -p ${output_dir}

# 出力ファイルの指定
sam_output="${output_dir}/${sample}.sam"
log_output="${output_dir}/${sample}.summary_log.txt"

# Bowtie2の実行
bowtie2 \
    --very-sensitive \
    --local \
    --no-mixed \ 
    --maxins 1000 \
    -k 1 \
    -q \
    -x ${bowtie2_index} \
    -1 ${read1} \
    -2 ${read2} \
    -S ${sam_output} \
    --threads ${SLURM_CPUS_PER_TASK} \
    2> ${log_output}

echo "Mapping completed for ${sample}. Logs saved to ${log_output}."
```
#### マウス

bowtie2_mm.sh

```
#!/bin/bash
#SBATCH --job-name=bowtie2_mm
#SBATCH --output=bowtie2_mm%A_%a.out
#SBATCH --error=bowtie2_mm%A_%a.err 
#SBATCH --partition=batch  
#SBATCH --cpus-per-task=8 
#SBATCH --mem=16G 
#SBATCH --array=1-8

module load bowtie2

# インデックス
bowtie2_index="~/osc-fs/REFERENCES/bowtie2/mm10/mm10"

# サンプルリスト
samples=("10_L1" "3_L1" "4_L1" "5_L1" "6_L1" "7_L1" "8_L1" "9_L1")
sample=${samples[$SLURM_ARRAY_TASK_ID-1]}

# 入力ファイルの指定
read1="${sample}_1.trimmed.fq.gz"
read2="${sample}_2.trimmed.fq.gz"

# 出力ディレクトリの作成
output_dir="bowtie2_results"
mkdir -p ${output_dir}

# 出力ファイルの指定
sam_output="${output_dir}/${sample}.sam"
log_output="${output_dir}/${sample}.summary_log.txt"

# Bowtie2の実行
bowtie2 \
    --very-sensitive \
    --local \
    --no-mixed \ 
    --maxins 1000 \
    -k 1 \
    -q \
    -x ${bowtie2_index} \
    -1 ${read1} \
    -2 ${read2} \
    -S ${sam_output} \
    --threads ${SLURM_CPUS_PER_TASK} \
    2> ${log_output}

echo "Mapping completed for ${sample}. Logs saved to ${log_output}."
```

run script
```
sbatch bowtie2_hs.sh
sbatch bowtie2_mm.sh
```
<br/>

*⭐︎ Bowtie2オプションの意味*
--very-sensitive ：高感度なsensitiveモード（でも遅い）
--local  ：ローカルアライメント（末端のマッチしない部分を切り落とす）
--no-mixed ：ペアのうち一方だけがマップするような「片割れマップ」を無視する
--maxins 1000 ：ペアリードの最大インサートサイズを1000bpに設定（Cut&Runではよく使う）
-k 1 ：各リードに対して最も良いマッピング1本だけを出力する
-q ：FASTQ形式の入力ファイル
2> ${log_output}　：Log出力


<br>

---

### 4. samtoolsでソートとbamCoverage

samtools.sh
```
#!/bin/bash
#SBATCH -n 4

filename=$1

samtools view -bS $filename.sam > $filename.bam

samtools sort $filename.bam -o $filename.sorted.bam

samtools index $filename.sorted.bam

bamCoverage -b $filename.sorted.bam -o $filename.bw --normalizeUsing CPM --binSize 10
```

<br>

---

### 5. NGSplot

テキストファイルの作成

hg38_tss.txt
```
echo 'S01_JHH7_mycn.sorted.bam -1 "JHH7_mycn"'>Human_tss.txt
echo 'S02_JHH7_input.sorted.bam -1 "JHH7_input"'>>Human_tss.txt
echo 'S11_HC_mycn..sorted.bam -1 "HC_mycn"'>>Human_tss.txt
echo 'S12_HC_input.sorted.bam -1 "HC_input"'>>Human_tss.txt
```

mm10_tss.txt
```
echo 'S03_A_mycn.sorted.bam -1 "A_mycn"'>Mouse_tss.txt
echo 'S04_B_mycn.sorted.bam -1 "B_mycn"'>>Mouse_tss.txt
echo 'S05_C_mycn.sorted.bam -1 "C_mycn"'>>Mouse_tss.txt
echo 'S06_D_mycn.sorted.bam -1 "D_mycn"'>>Mouse_tss.txt
echo 'S07_A_input.sorted.bam -1 "A_input"'>>Mouse_tss.txt
echo 'S08_B_input.sorted.bam -1 "B_input"'>>Mouse_tss.txt
echo 'S09_C_input.sorted.bam -1 "C_input"'>>Mouse_tss.txt
echo 'S10_D_input.sorted.bam -1 "D_input"'>>Mouse_tss.txt
```

run script
Human
```
ngs.plot.r -G hg38 -R tss -C Human_tss.txt -O Human_tss_1500 -L 1500 
ngs.plot.r -G hg38 -R genebody -C Human_tss.txt -O Human_genebody_1500 -L 1500
```
![スクリーンショット 2025-02-27 16.27.49](https://hackmd.io/_uploads/Bkpf496qyx.png)


マウス
```
ngs.plot.r -G mm10 -R tss -C Mouse_tss.txt -O Mouse_tss_1500 -L 1500 
ngs.plot.r -G mm10 -R genebody -C Mouse_tss.txt -O Mouse_genebody_1500 -L 1500
```
![スクリーンショット 2025-02-27 16.40.53](https://hackmd.io/_uploads/ryhGvqTc1e.png)


<br>

---

### 6. macs2 callpeak 

macs2_peakcall.sh
```
#!/bin/bash
#SBATCH -n 8

treatment=$1
control=$2
genome=$3
filename=$4

macs2 callpeak -t $treatment -c $control -f BAM -g $genome -n $filename --outdir macs2 --keep-dup auto -q 0.05 -B
```
run script (B-mycnとC-mycnに絞って)

```
sbatch -J mp0408 macs2_peakcall.sh S04_B_mycn.sorted.bam S08_B_input.sorted.bam mm B_mycn_vs_input

sbatch -J mp0509 macs2_peakcall.sh S05_C_mycn.sorted.bam S09_C_input.sorted.bam mm C_mycn_vs_input
```

<br>

---

### 7. HOMER でモチーフ検索

homer_motif.sh
```
#!/bin/bash
#SBATCH -n 8

file=$1
genome=$2
output=$3

findMotifsGenome.pl $file $genome $output -size 200 -mask
```
run script (B-mycnとC-mycnに絞って)
```
sbatch -J hoB homer_motif.sh B_mycn_vs_input_summits.bed mm10 homer_motif_B_mycn_vs_input

sbatch -J hoC homer_motif.sh C_mycn_vs_input_summits.bed mm10 homer_motif_C_mycn_vs_input
```
<br>

---

### 8. ChIPpeakAnnoでアノテーション (RStudio)


```r 
library(ChIPpeakAnno)

#各サンプル読み込み
JHH7 <- toGRanges("JHH7_mycn_vs_input_peaks.narrowPeak", format="narrowPeak", header=FALSE)
Bmycn <- toGRanges("B_mycn_vs_B_input_peaks.narrowPeak", format="narrowPeak", header=FALSE)
Cmycn <- toGRanges("C_mycn_vs_C_input_peaks.narrowPeak", format="narrowPeak", header=FALSE)

# ピーク同士の重なりを調査する
ol <- findOverlapsOfPeaks(JHH7, Bmycn, Cmycn, connectedPeaks="keepAll") 

# 重なりをベン図として可視化する
makeVennDiagram(ol, connectedPeaks="keepFirstListConsistent",  fill=c("#CC79A7", "#56B4E9", "#F0E442"), col=c("#D55E00", "#0072B2", "#E69F00"), cat.col=c("#D55E00", "#0072B2", "#E69F00"))

---
                
>$p.value
     JHH7 Bmycn Cmycn         pval
[1,]    0     1     1 4.421496e-19
[2,]    1     0     1 1.000000e+00
[3,]    1     1     0 1.663610e-02

>$vennCounts
     JHH7 Bmycn Cmycn Counts count.JHH7 count.Bmycn count.Cmycn
[1,]    0     0     0      0          0           0           0
[2,]    0     0     1   1603          0           0        1603
[3,]    0     1     0     40          0          40           0
[4,]    0     1     1     14          0          14          14
[5,]    1     0     0     78         78           0           0
[6,]    1     0     1      0          0           0           0
[7,]    1     1     0      1          7           1           0
[8,]    1     1     1      0          0           0           0
attr(,"class")
[1] "VennCounts"
```
![スクリーンショット 2025-03-11 15.26.32](https://hackmd.io/_uploads/BJ2nDUaj1g.png)



```r 
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(TxDb.Hsapiens.UCSC.hg38.knownGene) #どっちかでいいかも

# GTFファイルから遺伝子アノテーションを取得
txdb_mm <- makeTxDbFromGFF("mm10.refGene.gtf", format="gtf") # マウス
txdb_hs <- makeTxDbFromGFF("hg38.refGene.gtf", format="gtf") # ヒト

# 遺伝子情報を取得
annoData_mm <- genes(txdb_mm)
annoData_hs <- genes(txdb_hs) 

# Annotate　JHH7
anno_JHH7 <- annotatePeakInBatch(JHH7, AnnotationData=annoData_hs)
# JHH7 割合
pie1(table(anno_JHH7$insideFeature), main="JHH7")
#　save
write.table(anno_JHH7, "JHH7_annotate.txt", sep="\t", quote=F)

# Annotate　B-mycn
anno_B <- annotatePeakInBatch(Bmycn, AnnotationData=annoData_mm)
# B-mycn 割合
pie1(table(anno_B$insideFeature), main="B_mycn")
#　save
write.table(anno_B, "B_mycn_annotate.txt", sep="\t", quote=F)

# Annotate　C-mycn
anno_C <- annotatePeakInBatch(Cmycn, AnnotationData=annoData_mm)
# C-mycn 割合
pie1(table(anno_C$insideFeature), main="C_mycn")
#　save
write.table(anno_C, "C_mycn_annotate.txt", sep="\t", quote=F)
```
![スクリーンショット 2025-03-11 15.26.06](https://hackmd.io/_uploads/ryapvI6okl.png)

<br>

---

### 9. IGV

Mycn
![igv_snapshot_mycn](https://hackmd.io/_uploads/Byqbmf-See.png)

Tfrc
![igv_snapshot_Tfrc](https://hackmd.io/_uploads/S1xmQM-rle.png)

Cyp7b1
![igv_snapshot_Cyp7b1](https://hackmd.io/_uploads/HkZVmzWrel.png)

Lgalsl
![igv_snapshot_Lgalsl](https://hackmd.io/_uploads/B1xIXzbBgl.png)

Zdhhc2
![igv_snapshot_Zdhhc2](https://hackmd.io/_uploads/rk5LQz-Sxe.png)

Rsbn1l
![igv_snapshot_Rsbm1l](https://hackmd.io/_uploads/rkNv7GZBel.png)

En2
![igv_snapshot_En2](https://hackmd.io/_uploads/B1kuXGWrgg.png)

Dap
![igv_snapshot_Dap](https://hackmd.io/_uploads/B1dOQfbHll.png)













