ALLHIC使用｜ HiC辅助基因组组装（三）

最新推荐文章于 2024-09-10 19:02:52 发布

生信技术

最新推荐文章于 2024-09-10 19:02:52 发布

阅读量4.2k

点赞数 8

分类专栏： HiC 基因组

本文链接：https://blog.youkuaiyun.com/m0_49960764/article/details/119460897

版权

基因组同时被 2 个专栏收录

8 篇文章

订阅专栏

HiC

4 篇文章

订阅专栏

安装

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

依赖软件

samtools v1.9+
bedtools
matplotlib v2.0+

写在前面

ALLHIC官网提供了很详尽的内容，以及完整的pipeline，所以这里我主要是用来理清楚其整体思路，记录一下。
建议使用软件务必参照官网

官网链接手册

整体流程

ALLHiC一共分为五步:pruning, partition, rescue, optimization, building

prune 步骤去除了等位基因之间的联系，因此同源染色体更易于单独分离。
partition 功能将修剪的bam文件作为输入，并根据Hi-C建议的链接对链接的contigs进行聚类，大概是沿着相同同源染色体在预设数量的分区中进行。
rescue 功能从原始未修剪的bam文件中搜索分区步骤中不涉及的contigs，并根据Hi-C信号密度将它们分配给特定的群集。
optimize 步骤采用每个分区，并优化所有contigs的顺序和方向。
build 步骤通过连接contigs来重建每个染色体

如下图所示：

ALLHiC ]

Explanation of Prune

同源四倍体基因组的示意图。四个同源染色体显示为不同的颜色（分别为蓝色，橙色，绿色和紫色）。染色体中的红色区域表示具有高度相似性的序列。
检测自身四倍体基因组中的Hi-C信号。黑色虚线表示折叠区域和未折叠区域contigs之间的Hi-C信号。粉色虚线表示单体型Hi-C链接，灰色虚线表示单体型Hi-C链接。在组装过程中，红色区域会因高度的序列相似性而崩溃；同时，如果其他区域之间存在大量差异，则会将它们分为不同的contigs。由于塌陷区域与来自不同单倍型的contigs在物理上相关，因此将在塌陷区域与所有其他未塌陷的contigs之间检测到Hi-C信号。
传统的Hi-C脚手架方法将检测来自不同单倍型和折叠区域的contigs中的信号，并将所有序列聚在一起。
修剪Hi-C信号：1-去除等位基因区域之间的信号；2-仅在折叠区域和未折叠contigs之间保留最强的信号。
基于修剪的Hi-C信息进行分区。理想情况下，根据修剪结果将contigs分为不同的组。

Prune

运行ALLHIC

数据准备

Contig 水平的基因组组装结果
Hi-C测序数据

将 Hi-C reads Map 到基因组草图

bwa index and samtools faidx 对基因组草图建立索引

bwa index -a bwtsw draft.asm.fasta  
samtools faidx draft.asm.fasta

Hi-C reads Aligning 到基因组草图上

bwa aln -t 24 draft.asm.fasta reads_R1.fastq.gz > sample_R1.sai  
bwa aln -t 24 draft.asm.fasta reads_R2.fastq.gz > sample_R2.sai  
bwa sampe draft.asm.fasta sample_R1.sai sample_R2.sai reads_R1.fastq.gz reads_R2.fastq.gz > sample.bwa_aln.sam

过滤SAM文件

PreprocessSAMs.pl sample.bwa_aln.sam draft.asm.fasta MBOI
(*)filterBAM_forHiC.pl sample.bwa_aln.REduced.paired_only.bam sample.clean.sam
samtools view -bt draft.asm.fasta.fai sample.clean.sam > sample.clean.bam

*Tip: skip this step if you are using bwa mem for alignment

Prune(option)

Usage:  

ALLHiC_prune -i Allele.ctg.table -b sample.clean.bam -r draft.asm.fasta

在这里需要输入Allele.ctg.table文件，这个文件需要自己获取，参考官方参考链接;

ALLHiC 依靠等位基因重叠群表 (Allele.ctg.table) 来去除嘈杂的 Hi-C 信号

Allele.ctg.table 的格式：
前两列是染色体 ID 和参考基因组的位置。（如果您不使用参考基因组来生成 Allele.ctg.table，您可以将它们保留为 NA）
第 3 到第 N 列是我们确定的等位基因重叠群。修剪步骤将删除等位基因重叠群之间的 Hi-C 链接读数。

方法一：基于 BLAST 结果鉴定等位基因重叠群

将目标基因组中的 CDS Blast到相关参考的 CDS 文件
注意：请在运行 BLAST 之前修改 cds 名称。cds 名称应与 GFF3 中存在的基因名称相同

$ blastn -query rice.cds -db Bd.cds -out rice_vs_Sb.blast.out -evalue 0.001 -outfmt 6 -num_threads 4 -num_alignments 1

移除同一性 < 60% 和覆盖率 < 80% 的blast hits

blastn_parse.pl -i rice_vs_Sb.blast.out -o Erice_vs_Sb.blast.out -q rice.cds-b 1 -c 0.6 -d 0.8

根据 BLAST 结果对等位基因进行分类

classify.pl -i Eblast.out -p 2 -r Sbicolor_313_v3.1.gene -g rice.gff3

方法二：基于 GMAP 的方法

来生成 Allele.ctg.table，它不需要的目标基因组的注释。
这个方法适合于大多数人，因为如果是De novo组装，一般做Hi-C辅助基因组组装时肯定没有进行基因预测没有gff以及cds、pep序列，所以就需要进行下面的操作。

详细命令请看以下链接

运行 gmap 得到一个 gff3 文件

gmap_build -D . -d DB target.genome
gmap -D . -d DB -t 12 -f 2 -n $N reference.cds.fasta > gmap.gff3

注意：

target.genome 是多倍体基因组组装的 contig 水平
$N 是你的目标基因组的倍性，例如 $N=4 如果它是四倍体
reference.cds.fasta 是编码序列二倍体基因组，可参考得到等位基因表

生成 allelic.ctg.table

perl gmap2AlleleTable.pl referenece.gff3

gmap2Alleletable.pl的脚本链接点击此处

gmap2Alleletable.pl内容如下

#!/usr/bin/perl -w

die "Usage: perl $0 ref.gff3\n" if(!defined ($ARGV[0]));
my $refGFF = $ARGV[0];
open(IN, "grep 'gene' gmap.gff3 |") or die"";
while(<IN>){
	chomp;
	my @data = split(/\s+/,$_);
	my $gene = $1 if(/Name=([^;\n]*)/);
	$infordb{$gene} .= $data[0]."	";
	}
close IN;

open(OUT, "> Allele.ctg.table") or die"";
open(IN, "awk '\$3==\"gene\"' $refGFF |") or die"";
while(<IN>){
	chomp;
	my @data = split(/\s+/,$_);
	my $gene = $1 if(/Name=(\S+)/);
	   $gene =~ s/;.*//g;
	next if(!exists($infordb{$gene}));
	my @tdb = split(/\s+/,$infordb{$gene});
	my %tmpdb = ();
	map {$tmpdb{$_}++} @tdb;
	print OUT "$data[0]	$data[3]	";
	map {print OUT "$_	"} keys %tmpdb;
	print OUT "\n";
	}
close IN;
close OUT;

Partition

ps :如果跳过了第四步，那么可以直接用第三步分析的结果sample.clean.bam

Usage:  

ALLHiC_partition -b prunning.bam -r draft.asm.fasta -e AAGCTT -k 16  

Parameters: 
      -h : help and usage.
      -b : prunned bam (optional, default prunning.bam)
      -r : draft.sam.fasta
      -e : 酶切位点(HindIII: AAGCTT; MboI: GATC)
      -k : 分组数量 根据HiC信号将不同的contig进行分组
      -m : minimum number of restriction sites (default, 25)

Rescue

#如果做了第四步，会生成 -c prunning.clusters.txt 以及 -i prunning.counts_AAGCTT.txt
ALLHiC_rescue -b sample.clean.bam -r draft.asm.fasta -c prunning.clusters.txt -i prunning.counts_AAGCTT.txt

Optimize

# 生成.clm文件
allhic extract sample.clean.bam draft.asm.fasta --RE AAGCTT  
# 优化
for i in group*.txt; do
    allhic optimize $i sample.clean.clm
done

Build

ALLHiC_build draft.asm.fasta

plot

ALLHiC_plot sample.clean.bam groups.agp chrn.list 500k pdf

作者问题回复

修剪步骤中折叠区域和嵌合重叠群之间有什么区别？
在多倍体基因组中，一些同源区域（即等位基因序列）高度相似。这些序列经常被组装成一个重叠群，因为组装者不能分离等位基因。这种区域是折叠区域。
另一方面，一些重叠群包含来自不同单倍型或非同源染色体的序列，可以将其视为chimeric contigs。ALLHiC 旨在最大限度地减少collapsed contigs的负面影响，我们的模拟数据显示 ALLHiC 能够tolerate ~20% 的collapsed contigs；然而，只有约 5% 的chimeric contigs.。
如何确定哪些染色体组是同源的？