日志===== 初始化环境 =====
===== 验证工具版本 =====
samtools 1
.22
.1
1
.2
.2
/public1/home/scb3201/anaconda3/envs/Hic/share/3d-dna/run-asm-pipeline
.sh --version
1
.2
.2
cooler, version 0
.10
.2
3
.2
.4
===== 开始Hi-C数据预处理 =====
-- 运行fastp质控 --
-- 创建BWA索引 --
-- 创建染色体大小
文件 --
-- 生成酶切位点
文件 --
===== 运行Juicer+3D-DNA流程 =====
-- 运行Juicer主流程 --
-- 运行3D-DNA组装 --
/public1/home/scb3201/anaconda3/envs/Hic/share/3d-dna/run-asm-pipeline
.sh /tmp/hic_6309618/juicer/aligned/merged_nodups
.txt /public1/home/scb3201/results/genome/Medaka/genome
.nextpolish
.medaka
.fasta --output /public1/home/scb3201/results/genome/Hic/3d_dna_assembly --threads 24 --rounds 3
version: 190716
-- 整理染色体级组装 --
错误: 3D-DNA输出
文件未生成!
“#!/bin/bash
#SBATCH -p amd_256
#SBATCH -N 1
#SBATCH -n 24
#SBATCH -t 48:00:00
#SBATCH --mem=250G
#SBATCH -J hic_assembly
#SBATCH -o /public1/home/scb3201/rawdata/Hic/logs/%j_hic
.out
#SBATCH -e /public1/home/scb3201/rawdata/Hic/logs/%j_hic
.err
##########################################################
# 1
. 环境配置与工具验证
##########################################################
echo "===== 初始化环境 ====="
# 确保日志目录存在
mkdir -p /public1/home/scb3201/rawdata/Hic/logs
mkdir -p /public1/home/scb3201/results/genome/Hic
# 激活conda环境
source ~/
.bashrc
conda activate Hic || {
echo "错误：激活Hic环境失败！"
exit 1
}
# 设置工具路径
export REF_
FASTA="/public1/home/scb3201/results/genome/Medaka/genome
.nextpolish
.medaka
.fasta"
export HIC_R1="/public1/home/scb3201/ONT_BGI_RNA-seq/not/HiC/DL_R1
.fq
.gz"
export HIC_R2="/public1/home/scb3201/ONT_BGI_RNA-seq/not/HiC/DL_R2
.fq
.gz"
export OUTPUT_DIR="/public1/home/scb3201/results/genome/Hic"
export TMP_DIR="/tmp/hic_${SLURM_JOB_ID}"
# 创建临时目录
mkdir -p ${TMP_DIR}
# 验证工具可用性
echo "===== 验证工具版本 ====="
tools=("bwa" "samtools" "juicer" "3d-dna" "yahs" "cooler" "hicstuff")
for tool in "${tools[@]}"; do
$tool --version | head -1 || echo "✗ $tool 不可用"
done
##########################################################
# 2
. Hi-C数据预处理
##########################################################
echo "===== 开始Hi-C数据预处理 ====="
# 2
.1 数据质控
echo "-- 运行fastp质控 --"
fastp -i ${HIC_R1} -I ${HIC_R2} \
-o ${OUTPUT_DIR}/clean_R1
.fq
.gz \
-O ${OUTPUT_DIR}/clean_R2
.fq
.gz \
--html ${OUTPUT_DIR}/fastp_report
.html \
--json ${OUTPUT_DIR}/fastp_report
.json \
-w 8 \
-j 4
# 2
.2 创建参考基因组索引
echo "-- 创建BWA索引 --"
bwa index -p ${TMP_DIR}/genome_index ${REF_
FASTA}
echo "-- 创建染色体大小
文件 --"
samtools faidx ${REF_
FASTA}
awk '{print $1 "\t" $2}' ${REF_
FASTA}
.fai > ${TMP_DIR}/chrom
.sizes
# 2
.3 生成酶切位点
文件
echo "-- 生成酶切位点
文件 --"
generate_site_positions
.py DpnII ${TMP_DIR}/genome_DpnII ${REF_
FASTA}
##########################################################
# 3
. Juicer + 3D-DNA染色体挂载
##########################################################
echo "===== 运行Juicer+3D-DNA流程 ====="
# 3
.1 Juicer预处理
echo "-- 运行Juicer主流程 --"
juicer
.sh \
-t 24 \
-g genome \
-d ${TMP_DIR}/juicer \
-s DpnII \
-z ${REF_
FASTA} \
-p ${TMP_DIR}/chrom
.sizes \
-y ${TMP_DIR}/genome_DpnII
.txt \
-D /public1/home/scb3201/anaconda3/envs/Hic/share/juicer/ \
1> ${OUTPUT_DIR}/juicer
.log 2>&1
# 3
.2 3D-DNA染色体挂载
echo "-- 运行3D-DNA组装 --"
3d-dna ${TMP_DIR}/juicer/aligned/merged_nodups
.txt \
${REF_
FASTA} \
--output ${OUTPUT_DIR}/3d_dna_assembly \
--threads 24 \
--rounds 3
# 3
.3 整理最终组装结果
echo "-- 整理染色体级组装 --"
ASSEMBLY_FILE="${OUTPUT_DIR}/3d_dna_assembly
.3d_dna
.asm
.fasta"
if [ -f "${ASSEMBLY_FILE}" ]; then
# 添加contig长度信息
seqkit fx2tab -n -l ${ASSEMBLY_FILE} > ${OUTPUT_DIR}/contig_lengths
.tsv
echo "染色体级组装完成: ${ASSEMBLY_FILE}"
else
echo "错误: 3D-DNA输出
文件未生成!"
exit 1
fi
##########################################################
# 4
. yahs + chromap备选流程
##########################################################
echo "===== 运行yahs+chromap备选流程 ====="
# 4
.1 chromap比对
echo "-- chromap比对 --"
chromap -t 24 -r ${REF_
FASTA} \
-1 ${OUTPUT_DIR}/clean_R1
.fq
.gz \
-2 ${OUTPUT_DIR}/clean_R2
.fq
.gz \
-o ${TMP_DIR}/aligned
.bam \
--barcode-format CR
# 4
.2 转换BAM为BED
echo "-- 转换BAM为BED --"
bedtools bamtobed -i ${TMP_DIR}/aligned
.bam > ${TMP_DIR}/aligned
.bed
# 4
.3 yahs组装
echo "-- yahs组装 --"
yahs ${TMP_DIR}/aligned
.bed ${REF_
FASTA} \
-o ${OUTPUT_DIR}/yahs_assembly \
-t 24
# 4
.4 整理yahs结果
if [ -f "${OUTPUT_DIR}/yahs_assembly
.yahs
.out
.bin" ]; then
cooler makebins \
${OUTPUT_DIR}/yahs_assembly
.yahs
.out
.bin \
${OUTPUT_DIR}/yahs_assembly
.cool \
-c2
echo "yahs组装完成: ${OUTPUT_DIR}/yahs_assembly
.cool"
fi
##########################################################
# 5
. Hi-C数据可视化
##########################################################
echo "===== 生成Hi-C可视化 ====="
# 5
.1 创建
.cool
文件
echo "-- 创建cooler矩阵 --"
cooler cload pairix \
${TMP_DIR}/chrom
.sizes:10000 \
<(pairtools parse -c ${TMP_DIR}/chrom
.sizes \
${TMP_DIR}/juicer/aligned/merged_nodups
.txt) \
${OUTPUT_DIR}/hic_matrix_10k
.cool
# 5
.2 平衡矩阵
echo "-- 矩阵平衡 --"
cooler balance ${OUTPUT_DIR}/hic_matrix_10k
.cool
# 5
.3 生成全基因组热图
echo "-- 全基因组热图 --"
hicPlotMatrix \
--matrix ${OUTPUT_DIR}/hic_matrix_10k
.cool \
--out ${OUTPUT_DIR}/whole_genome_heatmap
.png \
--log \
--dpi 300
# 5
.4 生成染色体交互热图
echo "-- 染色体交互热图 --"
for chrom in $(cut -f1 ${TMP_DIR}/chrom
.sizes | head -5); do
hicPlotTAD \
--matrix ${OUTPUT_DIR}/hic_matrix_10k
.cool \
--out ${OUTPUT_DIR}/chr${chrom}_tad
.png \
--chromosome ${chrom} \
--dpi 300
done
# 5
.5 生成CoolBox可视化
echo "-- 交互式可视化 --"
cat > ${OUTPUT_DIR}/generate_coolbox
.py << EOF
from coolbox
.api import *
clr = Cool("${OUTPUT_DIR}/hic_matrix_10k
.cool")
frame = XAxis() + clr + CoolValue()
frame
.show("${OUTPUT_DIR}/interactive_hic
.html")
EOF
python ${OUTPUT_DIR}/generate_coolbox
.py
##########################################################
# 6
. 结果整理与清理
##########################################################
echo "===== 整理最终结果 ====="
# 6
.1 组装评估
echo "-- 组装评估 --"
assembly-stats ${ASSEMBLY_FILE} > ${OUTPUT_DIR}/assembly_stats
.txt
# 6
.2 创建结果报告
cat > ${OUTPUT_DIR}/report
.html << EOF
<html>
<head><title>Hi-C组装报告</title></head>
<body>
<h1>Hi-C染色体挂载结果报告</h1>
<h2>组装统计</h2>
<pre>$(cat ${OUTPUT_DIR}/assembly_stats
.txt)</pre>
<h2>可视化结果</h2>
<p><strong>全基因组热图:</strong><br>
<
img src="whole_genome_heatmap
.png" width="80%"></p>
<h3>染色体TAD结构</h3>
$(for f in ${OUTPUT_DIR}/chr*_tad
.png; do
echo "<p>$(basename $f):<br><
img src=\"$(basename $f)\" width=\"60%\"></p>"
done)
<h2>交互式可视化</h2>
<p><a href="interactive_hic
.html">点击查看交互式Hi-C浏览器</a></p>
</body>
</html>
EOF
# 6
.3 清理临时
文件
echo "-- 清理临时
文件 --"
rm -rf ${TMP_DIR}
echo "===== Hi-C分析完成! 结果保存在: ${OUTPUT_DIR} ====="
”
这篇博客介绍了GATK工具中的BwaMemIndexImageCreator如何用于创建GATK与BWA所需的.reference.fasta.img文件。这个过程涉及到对FASTA格式的基因组参考序列进行操作，生成的.index.img文件对于生物信息学中的序列比对至关重要。
扫一扫
480
到【灌水乐园】发言
