不過(guò)濾，用測(cè)序公司直接給的BAM文件callpeak
#MACS2 軟件
環(huán)境還是 conda activate m6a
~/m6a ./singlepeakcall.sh
singlepeakcall.sh：
#!/bin/bash
# 20250419 m6Am-seq peak calling
input_dir="/home/data/t210424/m6a/"
output_dir="/home/data/t210424/m6a/peakcall250419"
mkdir -p "$output_dir"
groups=("sh" "Sham" "SNI")
for group in "${groups[@]}"; do
? ? group_output_dir="$output_dir/${group}"
? ? mkdir -p "$group_output_dir"
? ? for ((i=1; i<=3; i++)); do
? ? ? ? ip_file="${input_dir}/RIP_${group}_${i}.bam"
? ? ? ? input_file="${input_dir}/In_${group}_${i}.bam"
? ? ? ? macs2 callpeak \
? ? ? ? ? ? -t "$ip_file" \
? ? ? ? ? ? -c "$input_file" \
? ? ? ? ? ? -f BAM \
? ? ? ? ? ? -g 1000000000 \?# 小鼠轉(zhuǎn)錄組大?。?e9，人類(lèi)用3e9）
? ? ? ? ? ? --nomodel \? # 禁用自動(dòng)預(yù)測(cè)片段長(zhǎng)度（RNA固定長(zhǎng)度）
? ? ? ? ? ? --extsize 50 \? # 匹配m6Am-exo片段長(zhǎng)度（根據(jù)實(shí)驗(yàn)調(diào)整）
? ? ? ? ? ? --shift -10 \? # 5'端偏移10bp，對(duì)準(zhǔn)Cap m6Am位點(diǎn)
? ? ? ? ? ? --keep-dup all \? # 保留所有重復(fù)（RNA數(shù)據(jù)可能高重復(fù)）
? ? ? ? ? ? -q 0.05 \?# FDR閾值（默認(rèn)5%，可調(diào)至0.01更嚴(yán)格）
? ? ? ? ? ? --call-summits \?# 檢測(cè)峰頂（提高分辨率）
????????????--bdg \
????????????--SPMR \
? ? ? ? ? ? -n "${group}_${i}" \?# 輸出前綴
? ? ? ? ? ? --outdir "$group_output_dir"
? ? done
done
轉(zhuǎn)換 bedGraph → BigWig
# 安裝 UCSC 工具
conda install -c bioconda ucsc-bedgraphtobigwig
# 轉(zhuǎn)換處理組信號(hào)
bedGraphToBigWig sh_1_treat_pileup.bdg mm10.chrom.sizes sh_1_signal.bw
# 轉(zhuǎn)換對(duì)照組信號(hào)（可選）
bedGraphToBigWig sh_1_control_lambda.bdg mm10.chrom.sizes sh_1_background.bw
更寬松的篩選標(biāo)準(zhǔn)
./peakcall0419-2.sh
#!/bin/bash
# 20250419 m6Am-seq peak calling（寬松參數(shù)版）
input_dir="/home/data/t210424/m6a/"
output_dir="/home/data/t210424/m6a/peakcall250419_lenient"
chrom_sizes="/path/to/mm10.chrom.sizes"? # 需替換為實(shí)際路徑
mkdir -p "$output_dir"
groups=("sh" "Sham" "SNI")
for group in "${groups[@]}"; do
? ? group_output_dir="$output_dir/${group}"
? ? mkdir -p "$group_output_dir"
? ? for ((i=1; i<=3; i++)); do
? ? ? ? ip_file="${input_dir}/RIP_${group}_${i}.bam"
? ? ? ? input_file="${input_dir}/In_${group}_${i}.bam"
? ? ? ? # 運(yùn)行 MACS2（參數(shù)全面寬松化）
? ? ? ? macs2 callpeak \
? ? ? ? ? ? -t "$ip_file" \
? ? ? ? ? ? -c "$input_file" \
? ? ? ? ? ? -f BAM \
? ? ? ? ? ? -g 1000000000 \
? ? ? ? ? ? --nomodel \
? ? ? ? ? ? --extsize 30 \? ? ? ? ? # 更短延伸（適配小片段）
? ? ? ? ? ? --shift -5 \? ? ? ? ? ? # 減少偏移量
? ? ? ? ? ? --keep-dup all \
? ? ? ? ? ? -p 0.1 \? ? ? ? ? ? ? ? # 寬松p值閾值（代替q值）
? ? ? ? ? ? --min-length 20 \? ? ? # 允許20bp短峰
? ? ? ? ? ? --max-gap 30 \? ? ? ? ? # 合并間隔≤30bp的峰
? ? ? ? ? ? --call-summits \
? ? ? ? ? ? --bdg \
? ? ? ? ? ? --SPMR \
? ? ? ? ? ? -n "${group}_${i}" \
? ? ? ? ? ? --outdir "$group_output_dir"
done
有幾組報(bào)錯(cuò)了，手動(dòng)作的 SNI2和3
(m6a) t210424@stark11-server:~/m6a/peakcall250419_lenient$ macs2 callpeak -t RIP_SNI_3.bam -c In_SNI_3.bam -f BAM -g 1000000000 --nomodel --extsize 30 --shift -5 --keep-dup all -p 0.1 --min-length 20 --max-
gap 30? ? --call-summits? ? --bdg? ? --SPMR? --name SNI_3? --outdir /home/data/t210424/m6a/peakcall250419_lenient? ? --seed 43
然后bdg文件和染色質(zhì)坐標(biāo)不一致，一個(gè)是1一個(gè)是chr1，把bdg的坐標(biāo)加前綴
(m6a) t210424@stark11-server:~/m6a/peakcall250419_lenient$ sed 's/^/chr/' /home/data/t210424/m6a/peakcall250419_lenient/SNI_2_treat_pileup.bdg > SNI_2_treat_pileup_chr.bdg
(m6a) t210424@stark11-server:~/m6a/peakcall250419_lenient$ bedGraphToBigWig /home/data/t210424/m6a/peakcall250419_lenient/SNI_2_treat_pileup_chr.bdg /home/data/t210424/Genome/Mouse/mm10.chrom.sizes SNI_2_signal.bw
發(fā)現(xiàn)13號(hào)染色體超界
End coordinate 120423434 bigger than chr13 size of 120421639 line 3495526 of /home/data/t210424/m6a/peakcall250419_lenient/SNI_2_treat_pileup_chr.bdg
修剪
awk '
BEGIN {OFS="\t"}
NR==FNR {size[$1]=$2; next}
$1 in size {
? ? if ($3 > size[$1]) $3 = size[$1];
? ? if ($2 < size[$1] && $3 <= size[$1]) print $0
}' /home/data/t210424/Genome/Mouse/mm10.chrom.sizes \
? SNI_2_treat_pileup_chr.bdg > SNI_2_treat_pileup_chr_fixed.bdg
修剪完再次轉(zhuǎn)換
(m6a) t210424@stark11-server:~/m6a/peakcall250419_lenient$ awk ' BEGIN {OFS="\t"}
NR==FNR {size[$1]=$2; next}
$1 in size {
? ? if ($3 > size[$1]) $3 = size[$1];
? ? if ($2 < size[$1] && $3 <= size[$1]) print $0
}' /home/data/t210424/Genome/Mouse/mm10.chrom.sizes? ? SNI_3_treat_pileup_chr.bdg > SNI_3_treat_pileup_chr_fixed.bdg
(m6a) t210424@stark11-server:~/m6a/peakcall250419_lenient$ bedGraphToBigWig SNI_3_treat_pileup_chr_fixed.bdg /home/data/t210424/Genome/Mouse/mm10.chrom.sizes SNI_3_signal.bw
