文件夹/home/data/t210424/m
1.QC是传统QC
2.trim
/home/data/t210424/m/trim.sh
# 定义路径和参数
raw_dir="/home/data/t210424/m"
trimmed_dir="/home/data/t210424/m/trimmed_data"
adapter_file="/home/data/t210424/miniconda3/envs/m6a/share/trimmomatic-0.39-2/adapters/TruSeq3-PE-2.fa"
for i in {1..9}; do
# RIP 样本
trimmomatic PE -threads 8 -phred33 \
"${raw_dir}/RIP_${i}.R1.raw.fastq.gz" \
"${raw_dir}/RIP_${i}.R2.raw.fastq.gz" \
"${trimmed_dir}/RIP_${i}.R1.trim.fq.gz" \
"${trimmed_dir}/RIP_${i}.R1.unpaired.fq.gz" \
"${trimmed_dir}/RIP_${i}.R2.trim.fq.gz" \
"${trimmed_dir}/RIP_${i}.R2.unpaired.fq.gz" \
"ILLUMINACLIP:${adapter_file}:2:30:10" \
LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:50
trimmomatic PE -threads 8 -phred33 \
"${raw_dir}/In_${i}.R1.raw.fastq.gz" \
"${raw_dir}/In_${i}.R2.raw.fastq.gz" \
"${trimmed_dir}/In_${i}.R1.trim.fq.gz" \
"${trimmed_dir}/In_${i}.R1.unpaired.fq.gz" \
"${trimmed_dir}/In_${i}.R2.trim.fq.gz" \
"${trimmed_dir}/In_${i}.R2.unpaired.fq.gz" \
"ILLUMINACLIP:${adapter_file}:2:30:10" \
LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:50
done
参考基因组和注释都是新的
# 创建参考基因组目录
mkdir -p /path/to/genome/mm10
cd /path/to/genome/mm10
# 下载基因组 FASTA 文件
wget ftp://ftp.ensembl.org/pub/release-99/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
gunzip Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
# 下载基因注释 GTF 文件
wget ftp://ftp.ensembl.org/pub/release-99/gtf/mus_musculus/Mus_musculus.GRCm38.99.gtf.gz
gunzip Mus_musculus.GRCm38.99.gtf.gz
Hisat索引
hisat2-build -p 8 \
Mus_musculus.GRCm38.dna.primary_assembly.fa \
mm10_hisat2_index
比对
#!/bin/bash
# 定义路径
trimmed_dir="/home/data/t210424/m/trimmed_data"
align_dir="/home/data/t210424/m/alignment"
genome_index="/home/data/t210424/Genome/Mouse/mm10_hisat2_index"
# 创建输出目录
mkdir -p "${align_dir}"
# 定义样本前缀列表（RIP 和 Input）
samples=("RIP" "In")
# 循环处理所有样本（1-9）
for prefix in "${samples[@]}"; do
for i in {1..9}; do
# 定义输入文件和输出前缀
input_r1="${trimmed_dir}/${prefix}_${i}.R1.trim.fq.gz"
input_r2="${trimmed_dir}/${prefix}_${i}.R2.trim.fq.gz"
output_prefix="${align_dir}/${prefix}_${i}"
# 检查输入文件是否存在
if [[ ! -f "${input_r1}" || ! -f "${input_r2}" ]]; then
echo "错误: 输入文件缺失 ${prefix}_${i}"
exit 1
fi
# 比对命令（HISAT2）
hisat2 -p 8 --dta -x "${genome_index}" \
-1 "${input_r1}" \
-2 "${input_r2}" \
-S "${output_prefix}.sam" 2> "${output_prefix}.hisat2.log"
# 转换为排序后的 BAM 文件
samtools sort -@ 8 -o "${output_prefix}.sorted.bam" "${output_prefix}.sam"
samtools index -@ 8 "${output_prefix}.sorted.bam"
# 清理临时文件
rm -f "${output_prefix}.sam"
echo "完成处理样本: ${prefix}_${i}"
done
done
echo "所有样本比对完成！结果保存在: ${align_dir}"
callpeak
传统版
#!/bin/bash
align_dir="/home/data/t210424/m/alignment"
peak_dir="/home/data/t210424/m/peakcall"
chrom_sizes="/home/data/t210424/Genome/Mouse/mm10.chrom.sizes"
mkdir -p "${peak_dir}"
for i in {1..9}; do
macs2 callpeak \
-t "${align_dir}/RIP_${i}.sorted.bam" \
-c "${align_dir}/In_${i}.sorted.bam" \
-f BAMPE -g mm -q 0.05 \
--nomodel --extsize 200 \
--outdir "${peak_dir}" \
-n "RIP_${i}_peaks" 2>&1 | tee "${peak_dir}/RIP_${i}_callpeak.log
# 生成标准化后的 bedGraph 用于 IGV
macs2 bdgcmp -t ${peak_dir}/RIP_${i}_peaks_treat_pileup.bdg \
-c ${peak_dir}/RIP_${i}_peaks_control_lambda.bdg \
-m FE --o-prefix ${peak_dir}/RIP_${i}_FE
# 转换为 bigWig
bedGraphToBigWig ${peak_dir}/RIP_${i}_FE.bdg ${chrom_sizes} ${peak_dir}/RIP_${i}_FE.bw
done
没有峰
宽松版p=0.1
#!/bin/bash
align_dir="/home/data/t210424/m/alignment"
peak_dir="/home/data/t210424/m/peakcallmore"
chrom_sizes="/home/data/t210424/Genome/Mouse/mm10.chrom.sizes"
mkdir -p "${peak_dir}"
for i in {1..9}; do
macs2 callpeak \
-t "${align_dir}/RIP_${i}.sorted.bam" \
-c "${align_dir}/In_${i}.sorted.bam" \
-f BAMPE \
-g mm \
--nolambda \
--keep-dup all \
--bw 500 \
--slocal 5000 \
--llocal 10000 \
--nomodel \
--min-length 20 \
--max-gap 30 \
--extsize 100 \
--bdg \
--SPMR \
--outdir "${peak_dir}" \
-n "RIP_${i}_peaks" 2>&1 | tee "${peak_dir}/RIP_${i}_callpeak.log"
# 生成标准化后的 bedGraph 用于 IGV
macs2 bdgcmp -t ${peak_dir}/RIP_${i}_peaks_treat_pileup.bdg \
-c ${peak_dir}/RIP_${i}_peaks_control_lambda.bdg \
-m FE --o-prefix ${peak_dir}/RIP_${i}_FE
# 转换为 bigWig
bedGraphToBigWig ${peak_dir}/RIP_${i}_FE.bdg ${chrom_sizes} ${peak_dir}/RIP_${i}_FE.bw
done
还是没有想要的，不行就加shift参数手动偏移--shift -5 --extsize 30
#!/bin/bash
align_dir="/home/data/t210424/m/alignment"
peak_dir="/home/data/t210424/m/peakcallmoreandmore"
chrom_sizes="/home/data/t210424/Genome/Mouse/mm10.chrom.sizes"
mkdir -p "${peak_dir}"
for i in {1..9}; do
macs2 callpeak \
-t "${align_dir}/RIP_${i}.sorted.bam" \
-c "${align_dir}/In_${i}.sorted.bam" \
-f BAMPE \
-g mm \
--nolambda \
--keep-dup all \
--bw 500 \
--nomodel \
-p 0.1 \
--min-length 20 \
--max-gap 30 \
--extsize 30 \
--shift -5 \
--bdg \
--SPMR \
--outdir "${peak_dir}" \
-n "RIP_${i}_peaks" 2>&1 | tee "${peak_dir}/RIP_${i}_callpeak.log"
# 生成标准化后的 bedGraph 用于 IGV
macs2 bdgcmp -t ${peak_dir}/RIP_${i}_peaks_treat_pileup.bdg \
-c ${peak_dir}/RIP_${i}_peaks_control_lambda.bdg \
-m FE --o-prefix ${peak_dir}/RIP_${i}_FE
# 转换为 bigWig
bedGraphToBigWig ${peak_dir}/RIP_${i}_FE_FE.bdg ${chrom_sizes} ${peak_dir}/RIP_${i}_FE.bw
done
