染色体挂载

如果有已发表的染色体水平的参考基因组，且参考基因组组装质量尚可，也没有大的结构变异，则可以通过共线性进行挂载，比如泛基因组。项目地址：https://github.com/malonge/RagTag。RaGOO已不再更新，直接使用RagTag

# 安装
conda create -n ragtag -c conda-forge -c bioconda ragtag
conda activate ragtag

# 运行
ragtag.py \
scaffold \
-t 16 \
-o ragtag_output \
-u \ # 未挂载序列加标记 
--aligner nucmer \ # 比对软件 
ref.fasta draft.fasta

输出结果：

ragtag.scaffold.agp # 序列连接关系文件
ragtag.scaffold.asm.paf
ragtag.scaffold.asm.paf.log
ragtag.scaffold.confidence.txt
ragtag.scaffold.err
ragtag.scaffold.fasta # 结果文件
ragtag.scaffold.stats

一般测基因组大小的100X

对于组装指标一般，序列条数较多的情况，推荐使用 ALLHiC 进行 Hi-C 数据染色体挂载。可以指定染色体条数。可用于多倍体物种的挂载，需要提供近缘物种的参考基因组。项目地址：https://github.com/tangerzhang/ALLHiC

# 安装
git clone https://github.com/tangerzhang/ALLHiC
cd ALLHiC
chmod +x bin/*
chmod +x scripts/*
export PATH=/your/path/to/ALLHiC/scripts/:/your/path/to/ALLHiC/bin/:$PATH

可以基于 hic 数据的比对情况对基因组中可能存在的组装错误进行打断，但该操作会降低序列连续性。该步骤可选，一般先直接挂载，然后手动调整即可。如果发现问题太多，则打断后重新挂载。

ref=../data/genome.fa
R1=../data/HiC_R1.fastq.gz
R2=../data/HiC_R2.fastq.gz
threads=10

# index reference genome
ln -s ${ref} ./seq.fasta
bwa index seq.fasta
samtools faidx seq.fasta

# mapping
bwa mem -SP5M -t $threads seq.fasta $R1 $R2 \
| samtools view -hF 256 - \
| samtools sort -@ $threads -o sorted.bam -T tmp.ali
samtools index sorted.bam

# correct contig
ALLHiC_corrector \
-m sorted.bam \
-r seq.fasta \
-o corrected.fasta \ -t $threads

挂载染色体

ref=../data/genome.fa
R1=../data/HiC_R1.fastq.gz
R2=../data/HiC_R2.fastq.gz
threads=10
group_count=5
enzyme=GATC # HindIII: AAGCTT; MboI: GATC

# index reference genome
ln -s ${ref}  ./seq.HiCcorrected.fasta
bwa index seq.HiCcorrected.fasta
samtools faidx seq.HiCcorrected.fasta
bwa mem -SP5M -t $threads seq.HiCcorrected.fasta $R1 $R2 \
| samtools view -hF 256 - \
| samtools sort -@ $threads -o sample.bwa_mem.bam -T tmp.ali

# filter bam
samtools view -bq 40 sample.bwa_mem.bam | \
samtools view -bt seq.HiCcorrected.fasta.fai > sample.unique.bam
PreprocessSAMs.pl sample.unique.bam seq.HiCcorrected.fasta $enzyme

# partition
ALLHiC_partition \
-r seq.HiCcorrected.fasta \
-e $enzyme \
-k $group_count \
-b sample.unique.REduced.paired_only.bam

# optimize
rm cmd.list
for((K=1;K<=$group_count;K++));do echo "allhic optimize \
sample.unique.REduced.paired_only.counts_${enzyme}.${group_\ count}g${K}.txt \
sample.unique.REduced.paired_only.clm" >> cmd.list;done
ParaFly -c cmd.list -CPU $threads 

# build
ALLHiC_build  seq.HiCcorrected.fasta

# plot
samtools faidx groups.asm.fasta
cut -f1,2 groups.asm.fasta.fai|grep sample > chrn.list
ALLHiC_plot sample.bwa_mem.bam groups.agp chrn.list 50k pdf

自动挂载的染色体可能出现染色体划分不正确的情况，可以生成 juicer- box 格式文件，使用 juicerbox 进行手动纠错。

juicebox_scripts网址：https://github.com/phasegenomics/juicebox_scripts

# 基于read name 排序bam
samtools sort -n -@ 12 -o aligned.sort_name.bam  sample.bwa_mem.bam

# 生成juicer软件的merged_nodups.txt文件
matlock bam2 juicer aligned.sort_name.bam merged_nodups.txt

# agp 格式转3d-dna的assembly 格式
agp2assembly.py groups.agp groups.assembly

# 生成juicebox的二进制hic文件
bash 3d-dna/visualize/run-assembly-visualizer.sh \ 
-q 1 -p true groups.assembly merged_nodups.txt

网址:https://github.com/aidenlab/juicer/

# 手动安装
git clone https://github.com/theaidenlab/juicer.git
# <myJuicerDir> 为想要运行juicer的目录
cp juicer/CPU <myJuicerDir>/scripts
cd scripts/common
wget https://github.com/aidenlab/Juicebox/releases/download/\
  v2.20.00/juicer_tools.2.20.00.jar
ln -s juicer_tools.2.20.00.jar  juicer_tools.jar

使用 juicer 生成 3d-dna 输入文件

ref=../../data/genome.fa
R1=../../data/HiC_R1.fastq.gz
R2=../../data/HiC_R2.fastq.gz
enzyme=MboI  # or HindIII

# 准备基因组相关文件
# 包括bwa index及染色体长度文件
ln -s $ref genome.fa
singularity exec ../../../container/Assembly202306.sif bwa \
index genome.fa
singularity exec ../../../container/Assembly202306.sif seqkit fx2tab \
-l -n -i genome.fa > chr.size

# 准备酶切位点文件
generate_site_positions.py $enzyme genome genome.fa

# 准备juicer script目录
cp -r /opt/juicer/CPU/   ./scripts

# 准备测序数据，不要修改R*.fastq.gz后缀
mkdir fastq 
ln -s PATH/HiC_R1.fastq.gz fastq/HiC_R1.fastq.gz
ln -s PATH/HiC_R2.fastq.gz  fastq/HiC_R2.fastq.gz

# 运行juicer
bash \
./scripts/juicer.sh \
-y genome_${enzyme}.txt \ # 酶切位点文件 
-z genome.fa \ # 基因组文件
-s $enzyme \ # 酶
-p chr.size \ # 染色体长度文件
-D ./ \ # scripts 和 fastq文件所在目录 -t 20 \
--assembly # 生成3d-dna输入文件

网址:https://github.com/aidenlab/3d-dna

conda install -c conda-forge -c bioconda 3d-dna

运行

3d-dna/run-asm-pipeline.sh \
-r 1 \ # 纠错次数，0表示不纠错，默认为2 
genome.fa juicer/aligned/merged_nodups.txt

结果文件:

genome.FINAL.fasta :3d-dna 的挂载结果文件
genome.final.assembly :用于 juicebox 输入的 assembly 文件
genome.final.hic :用于 juicebox 输入的 hic 文件

如果组装结果连续性较好，只有几十条序列，可以不运行 3d-dna, 仅生成 juice-box 输入, 后续全部手动调图。

# 生成assembly文件
awk -f 3d-dna/utils/generate-assembly-file-from-fasta.awk \ 
genome.fa > original.assembly

# 生成hic
bash 3d-dna/visualize/run-assembly-visualizer.sh \ 
original.assembly juicer/aligned/merged_nodups.txt

标准化一般选择“Balanced”。需要掌握的操作：

位置调整：Shift+点击选中某区块，然后将鼠标箭头放在需要重新放置的2个区域的十字附近，变为指向焦点处的箭头
旋转：鼠标箭头放在该区域的右上角，箭头变为圆圈
打断：鼠标箭头放在想要打断的对角线位置，箭头变为剪刀
合并：鼠标箭头放在需要合并的2个区域的十字附近，箭头变为一个直角

手动调整后，使用下面程序生成最终的 fasta 文件和 hic 文件

run-asm-pipeline-post-review.sh \
-q 1 -p true --build-gapped-map --sort-output \
-r review.assembly genome.fa merged_nodups.txt

2023-11-02

染色体挂载

基于已有参考基因组

基于HIC

ALLHIC

3D-DNA

运行 juicer

运行3d-dna

juicerbox

本页目录