不过滤，用测序公司直接给的BAM文件callpeak
#MACS2 软件
环境还是 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 \ # 小鼠转录组大小（1e9，人类用3e9）
--nomodel \ # 禁用自动预测片段长度（RNA固定长度）
--extsize 50 \ # 匹配m6Am-exo片段长度（根据实验调整）
--shift -10 \ # 5'端偏移10bp，对准Cap m6Am位点
--keep-dup all \ # 保留所有重复（RNA数据可能高重复）
-q 0.05 \ # FDR阈值（默认5%，可调至0.01更严格）
--call-summits \ # 检测峰顶（提高分辨率）
--bdg \
--SPMR \
-n "${group}_${i}" \ # 输出前缀
--outdir "$group_output_dir"
done
done
转换 bedGraph → BigWig
# 安装 UCSC 工具
conda install -c bioconda ucsc-bedgraphtobigwig
# 转换处理组信号
bedGraphToBigWig sh_1_treat_pileup.bdg mm10.chrom.sizes sh_1_signal.bw
# 转换对照组信号（可选）
bedGraphToBigWig sh_1_control_lambda.bdg mm10.chrom.sizes sh_1_background.bw
更宽松的筛选标准
./peakcall0419-2.sh
#!/bin/bash
# 20250419 m6Am-seq peak calling（宽松参数版）
input_dir="/home/data/t210424/m6a/"
output_dir="/home/data/t210424/m6a/peakcall250419_lenient"
chrom_sizes="/path/to/mm10.chrom.sizes" # 需替换为实际路径
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（参数全面宽松化）
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
有几组报错了，手动作的 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文件和染色质坐标不一致，一个是1一个是chr1，把bdg的坐标加前缀
(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
发现13号染色体超界
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
修剪完再次转换
(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
