###### 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
```
マウス
```
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
```
<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"
```
```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)
```
<br>
---
### 9. IGV
Mycn
Tfrc
Cyp7b1
Lgalsl
Zdhhc2
Rsbn1l
En2
Dap
Sign in with Wallet
Connect another wallet