Bowtie2 比对 参考基因集路径
FA=/refdir/database/server/reference/genome/mm10/mm10.fa
INDEX=/refdir/database/server/reference/index/bowtie/mm10
CHROMSIZES=/refdir/database/server/reference/genome/mm10/mm10.chrom.sizes
GTF=/home/data/t210424/Genome/Mouse/Mus_musculus.GRCm38.99.ucsc.gtf
工作路径/home/data/t210424/chipseq
样本信息samples.csv
环境检查
/home/data/t210424/chipseq/scripts/check_chip_env.sh
#!/usr/bin/env bash
set -euo pipefail
FA="/refdir/database/server/reference/genome/mm10/mm10.fa"
INDEX="/refdir/database/server/reference/index/bowtie/mm10"
CHROMSIZES="/refdir/database/server/reference/genome/mm10/mm10.chrom.sizes"
GTF="/home/data/t210424/Genome/Mouse/Mus_musculus.GRCm38.99.ucsc.gtf"
SAMPLE_SHEET="/home/data/t210424/chipseq/samples.csv"
echo "=== 检查软件 ==="
for x in bowtie2 samtools macs3 bamCoverage awk; do
if command -v "$x" >/dev/null 2>&1; then
echo "[OK] $x -> $(command -v "$x")"
else
echo "[MISSING] $x"
fi
done
echo
echo "=== 检查参考文件 ==="
for f in \
"$FA" \
"$CHROMSIZES" \
"$GTF" \
"${INDEX}.1.bt2" \
"${INDEX}.2.bt2" \
"${INDEX}.3.bt2" \
"${INDEX}.4.bt2" \
"${INDEX}.rev.1.bt2" \
"${INDEX}.rev.2.bt2" \
"$SAMPLE_SHEET"
do
if [[ -s "$f" ]]; then
echo "[OK] $f"
else
echo "[MISSING OR EMPTY] $f"
fi
done
echo
echo "=== 检查 samples.csv 列数 ==="
awk -F',' '{print NF}' "$SAMPLE_SHEET" | sort | uniq -c
echo
echo "=== 检查 samples.csv 前几行 ==="
head "$SAMPLE_SHEET"
echo
echo "=== 检查 samples.csv 里的 fastq 是否存在 ==="
awk -F',' 'NR>1 {print $5"\n"$6}' "$SAMPLE_SHEET" | while read -r f; do
[[ -z "$f" ]] && continue
if [[ -s "$f" ]]; then
echo "[OK] $f"
else
echo "[MISSING OR EMPTY] $f"
fi
done
echo
echo "=== 检查完成 ==="
QC
#Bowtie2比对参数→Samtools过滤 → 去重Markdup → 去除线粒体噪音 → Peak Calling bigWig信号归一
#H3K36me3所以写的board宽峰
/home/data/t210424/chipseq/scripts/run_chipseq_mm10_trim_fastp.sh
#!/usr/bin/env bash
set -euo pipefail
shopt -s nullglob
SAMPLE_SHEET="${1:-/home/data/t210424/chipseq/samples.csv}"
PROJECT_DIR="${2:-/home/data/t210424/chipseq/run_trimmed}"
FA="/refdir/database/server/reference/genome/mm10/mm10.fa"
INDEX="/refdir/database/server/reference/index/bowtie/mm10"
CHROMSIZES="/refdir/database/server/reference/genome/mm10/mm10.chrom.sizes"
GTF="/home/data/t210424/Genome/Mouse/Mus_musculus.GRCm38.99.ucsc.gtf"
THREADS=8
SORT_THREADS=4
SORT_MEM="1G"
GENOME_SIZE="mm"
PEAK_MODE="broad"
MACS_QVALUE="0.01"
MAPQ=30
REMOVE_CHRM=1
MAKE_BIGWIG=1
LOG_DIR="${PROJECT_DIR}/00_logs"
FASTP_DIR="${PROJECT_DIR}/01_fastp"
BAM_DIR="${PROJECT_DIR}/02_bam"
QC_DIR="${PROJECT_DIR}/03_qc"
PEAK_DIR="${PROJECT_DIR}/04_peaks"
BW_DIR="${PROJECT_DIR}/05_bigwig"
TMP_DIR="${PROJECT_DIR}/06_tmp"
mkdir -p "${PROJECT_DIR}" "${LOG_DIR}" "${FASTP_DIR}" "${BAM_DIR}" "${QC_DIR}" "${PEAK_DIR}" "${BW_DIR}" "${TMP_DIR}"
log() {
echo "[$(date '+%F %T')] $*"
}
die() {
echo "[ERROR] $*" >&2
exit 1
}
check_file() {
[[ -s "$1" ]] || die "文件不存在或为空 $1"
}
trim_spaces() {
local x="$1"
x="${x#"${x%%[![:space:]]*}"}"
x="${x%"${x##*[![:space:]]}"}"
echo "$x"
}
declare -A FINAL_BAM
declare -A GROUP_OF
declare -A KIND_OF
declare -A INPUT_OF
process_one_sample() {
local sample="$1"
local group="$2"
local kind="$3"
local input_id="$4"
local fq1="$5"
local fq2="$6"
local sample_dir="${BAM_DIR}/${sample}"
mkdir -p "${sample_dir}"
local fastp_log="${LOG_DIR}/${sample}.fastp.log"
local fastp_html="${FASTP_DIR}/${sample}.fastp.html"
local fastp_json="${FASTP_DIR}/${sample}.fastp.json"
local trim_r1="${TMP_DIR}/${sample}.trim.R1.fastq.gz"
local trim_r2="${TMP_DIR}/${sample}.trim.R2.fastq.gz"
local bowtie_log="${LOG_DIR}/${sample}.bowtie2.log"
local final_bam="${sample_dir}/${sample}.final.bam"
local final_bai="${final_bam}.bai"
local flagstat_txt="${QC_DIR}/${sample}.flagstat.txt"
local idxstats_txt="${QC_DIR}/${sample}.idxstats.txt"
local namesort_bam="${sample_dir}/${sample}.namesort.bam"
local fixmate_bam="${sample_dir}/${sample}.fixmate.bam"
local coordsort_bam="${sample_dir}/${sample}.coordsort.bam"
local rmdup_bam="${sample_dir}/${sample}.rmdup.bam"
local keep_chrom_txt="${sample_dir}/${sample}.keep_chroms.txt"
if [[ -s "${final_bam}" && -s "${final_bai}" ]]; then
log "跳过 ${sample}，已存在最终 BAM"
FINAL_BAM["${sample}"]="${final_bam}"
GROUP_OF["${sample}"]="${group}"
KIND_OF["${sample}"]="${kind}"
INPUT_OF["${sample}"]="${input_id}"
return
fi
check_file "${fq1}"
check_file "${fq2}"
log "开始 fastp ${sample}"
fastp \
-i "${fq1}" \
-I "${fq2}" \
-o "${trim_r1}" \
-O "${trim_r2}" \
--detect_adapter_for_pe \
-q 20 \
-l 20 \
-w "${THREADS}" \
-h "${fastp_html}" \
-j "${fastp_json}" \
> "${fastp_log}" 2>&1
check_file "${trim_r1}"
check_file "${trim_r2}"
log "开始 bowtie2 ${sample}"
bowtie2 \
--very-sensitive \
--no-mixed \
--no-discordant \
-X 2000 \
-x "${INDEX}" \
-1 "${trim_r1}" \
-2 "${trim_r2}" \
-p "${THREADS}" \
2> "${bowtie_log}" | \
samtools view \
-@ "${SORT_THREADS}" \
-b \
-f 2 \
-F 3852 \
-q "${MAPQ}" \
- | \
samtools sort \
-n \
-@ "${SORT_THREADS}" \
-m "${SORT_MEM}" \
-T "${TMP_DIR}/${sample}.namesort" \
-o "${namesort_bam}" \
-
rm -f "${trim_r1}" "${trim_r2}"
log "fixmate ${sample}"
samtools fixmate \
-m \
-@ "${SORT_THREADS}" \
"${namesort_bam}" \
"${fixmate_bam}"
rm -f "${namesort_bam}"
log "坐标排序 ${sample}"
samtools sort \
-@ "${SORT_THREADS}" \
-m "${SORT_MEM}" \
-T "${TMP_DIR}/${sample}.coordsort" \
-o "${coordsort_bam}" \
"${fixmate_bam}"
rm -f "${fixmate_bam}"
log "去重复 ${sample}"
samtools markdup \
-r \
-@ "${SORT_THREADS}" \
"${coordsort_bam}" \
"${rmdup_bam}"
rm -f "${coordsort_bam}"
if [[ "${REMOVE_CHRM}" -eq 1 ]]; then
log "去除 chrM ${sample}"
samtools index -@ "${SORT_THREADS}" "${rmdup_bam}"
samtools idxstats "${rmdup_bam}" | \
cut -f1 | \
grep -v '^chrM$' | \
grep -v '^\*$' > "${keep_chrom_txt}"
samtools view \
-@ "${SORT_THREADS}" \
-b \
"${rmdup_bam}" \
$(tr '\n' ' ' < "${keep_chrom_txt}") \
-o "${final_bam}"
samtools index -@ "${SORT_THREADS}" "${final_bam}"
rm -f "${rmdup_bam}" "${rmdup_bam}.bai" "${keep_chrom_txt}"
else
mv "${rmdup_bam}" "${final_bam}"
samtools index -@ "${SORT_THREADS}" "${final_bam}"
fi
log "QC 统计 ${sample}"
samtools flagstat -@ "${SORT_THREADS}" "${final_bam}" > "${flagstat_txt}"
samtools idxstats "${final_bam}" > "${idxstats_txt}"
if [[ "${MAKE_BIGWIG}" -eq 1 ]]; then
log "生成 bigWig ${sample}"
bamCoverage \
-b "${final_bam}" \
-o "${BW_DIR}/${sample}.CPM.bw" \
-p "${THREADS}" \
--binSize 50 \
--normalizeUsing CPM
fi
FINAL_BAM["${sample}"]="${final_bam}"
GROUP_OF["${sample}"]="${group}"
KIND_OF["${sample}"]="${kind}"
INPUT_OF["${sample}"]="${input_id}"
log "完成样本 ${sample}"
}
call_peaks_for_one_ip() {
local ip_sample="$1"
local input_sample="$2"
local ip_bam="${FINAL_BAM[${ip_sample}]}"
local input_bam="${FINAL_BAM[${input_sample}]}"
local outdir="${PEAK_DIR}/${ip_sample}"
local macs_log="${LOG_DIR}/${ip_sample}.macs2.log"
[[ -s "${ip_bam}" ]] || die "IP BAM 不存在 ${ip_bam}"
[[ -s "${input_bam}" ]] || die "Input BAM 不存在 ${input_bam}"
mkdir -p "${outdir}"
log "开始 call peak ${ip_sample} vs ${input_sample}"
if [[ "${PEAK_MODE}" == "broad" ]]; then
macs2 callpeak \
-t "${ip_bam}" \
-c "${input_bam}" \
-f BAMPE \
-g "${GENOME_SIZE}" \
-n "${ip_sample}" \
--outdir "${outdir}" \
--broad \
--broad-cutoff 0.1 \
--keep-dup all \
2> "${macs_log}"
else
macs2 callpeak \
-t "${ip_bam}" \
-c "${input_bam}" \
-f BAMPE \
-g "${GENOME_SIZE}" \
-n "${ip_sample}" \
--outdir "${outdir}" \
-q "${MACS_QVALUE}" \
--keep-dup all \
2> "${macs_log}"
fi
log "完成 peak calling ${ip_sample}"
}
log "开始读取样本表 ${SAMPLE_SHEET}"
check_file "${SAMPLE_SHEET}"
while IFS=',' read -r sample group kind input_id fq1 fq2; do
sample="$(trim_spaces "${sample}")"
group="$(trim_spaces "${group}")"
kind="$(trim_spaces "${kind}")"
input_id="$(trim_spaces "${input_id}")"
fq1="$(trim_spaces "${fq1}")"
fq2="$(trim_spaces "${fq2}")"
[[ -z "${sample}" ]] && continue
[[ "${sample}" =~ ^# ]] && continue
[[ "${sample}" == "sample_id" ]] && continue
process_one_sample "${sample}" "${group}" "${kind}" "${input_id}" "${fq1}" "${fq2}"
done < "${SAMPLE_SHEET}"
log "开始 peak calling"
while IFS=',' read -r sample group kind input_id fq1 fq2; do
sample="$(trim_spaces "${sample}")"
kind="$(trim_spaces "${kind}")"
input_id="$(trim_spaces "${input_id}")"
[[ -z "${sample}" ]] && continue
[[ "${sample}" =~ ^# ]] && continue
[[ "${sample}" == "sample_id" ]] && continue
if [[ "${kind}" == "ip" ]]; then
[[ "${input_id}" == "NA" || -z "${input_id}" ]] && die "IP 样本 ${sample} 没有对应 input_id"
[[ -n "${FINAL_BAM[${input_id}]:-}" ]] || die "IP 样本 ${sample} 对应的 input ${input_id} 未处理成功"
call_peaks_for_one_ip "${sample}" "${input_id}"
fi
done < "${SAMPLE_SHEET}"
log "写出最终 BAM 清单"
{
echo "sample_id,group,kind,input_id,final_bam"
while IFS=',' read -r sample group kind input_id fq1 fq2; do
sample="$(trim_spaces "${sample}")"
group="$(trim_spaces "${group}")"
kind="$(trim_spaces "${kind}")"
input_id="$(trim_spaces "${input_id}")"
[[ -z "${sample}" ]] && continue
[[ "${sample}" =~ ^# ]] && continue
[[ "${sample}" == "sample_id" ]] && continue
echo "${sample},${group},${kind},${input_id},${FINAL_BAM[${sample}]}"
done < "${SAMPLE_SHEET}"
} > "${PROJECT_DIR}/final_bam_manifest.csv"
log "全部完成"
log "结果目录"
log "FASTP ${FASTP_DIR}"
log "BAM ${BAM_DIR}"
log "QC ${QC_DIR}"
log "Peaks ${PEAK_DIR}"
log "bigWig ${BW_DIR}"
log "Logs ${LOG_DIR}"
chmod +x /home/data/t210424/chipseq/scripts/run_chipseq_mm10_trim_fastp.sh
后台运行+日志输出
chmod +x /home/data/t210424/chipseq/scripts/run_chipseq_mm10_trim_fastp.sh
cd /home/data/t210424/chipseq
nohup stdbuf -oL -eL bash /home/data/t210424/chipseq/scripts/run_chipseq_mm10_trim_fastp.sh \
/home/data/t210424/chipseq/samples.csv \
/home/data/t210424/chipseq/run_trimmed \
> /home/data/t210424/chipseq/run_trimmed.nohup.log 2>&1 &
echo $! > /home/data/t210424/chipseq/run_trimmed.pid
disown
进程查看
ps -fp $(cat /home/data/t210424/chipseq/run_trimmed.pid)
