染色体挂载

一.AllHiC挂载

1.二倍体染色体挂载

1.1 设置hic测序数据的存储路径

hicdir=hic测序得到的数据文件存放的位置
fq1=$hicdir/hic_1.fastq.gz
fq2=$hicdir/hic_2.fastq.gz

1.2 链接待挂载基因组文件

ln -s $hicdir/GCA_000222345.1_C24_2010-09-30.fa ./draft.asm.fasta

1.3 建立索引

bwa index draft.asm.fasta
#生成了draft.asm.fasta.amb、draft.asm.fasta.ann、draft.asm.fasta.bwt、draft.asm.fasta.pac、draft.asm.fasta.sa五个文件
samtools faidx draft.asm.fasta
#生成了draft.asm.fasta.fai一个文件

1.4 将HiC数据比对到基因组文件

bwa aln -t 12 draft.asm.fasta $fq1 > sample_R1.sai #生成sample_R1.sai一个文件
bwa aln -t 12 draft.asm.fasta $fq2 > sample_R2.sai #生成sample_R2.sai一个文件
bwa sampe draft.asm.fasta sample_R1.sai sample_R2.sai $fq1 $fq2 > sample.bwa_aln.sam #将sample_R1.sai、sample_R2.sai、draft.asm.fasta、$fq1和$fq2合并到一个sample.bwa_aln.sam文件

1.5 对比对结果进行过滤

PreprocessSAMs.pl sample.bwa_aln.sam draft.asm.fasta HINDIII #识别酶切位点,指明用的内切酶HINDIII(如果不知道需要咨询测序公司)
#依次生成sample.bwa_aln.bam、draft.asm.fasta.pos_of_AAGCTT.txt、draft.asm.fasta.near_AAGCTT.500.bed、sample.bwa_aln.REduced.bam、**sample.bwa_aln.REduced.paired_only.bam、sample.bwa_aln.REduced.paired_only.flagstat
filterBAM_forHiC.pl sample.bwa_aln.REduced.paired_only.bam sample.clean.sam #过滤生成sample.clean.sam文件
samtools view -bt draft.asm.fasta.fai sample.clean.sam > sample.clean.bam #将sample.clean.sam文件转换成sample.clean.bam**文件

1.6 将组装好的contig划分为不同的group

k=5 #设置染色体数目
ALLHiC_partition -b sample.clean.bam -r draft.asm.fasta -e HindIII -k $k
#生成文件sample.clean.counts_AAGCTT.txt、sample.clean.distribution.txt、sample.clean.clm、sample.clean.pairs.txt、**sample.clean.counts_AAGCTT.5g5.txt、sample.clean.counts_AAGCTT.5g4.txt、sample.clean.counts_AAGCTT.5g3.txt、sample.clean.counts_AAGCTT.5g2.txt、sample.clean.counts_AAGCTT.5g1.txt**、sample.clean.clusters.txt10个文件

1.7 提取CLM文件和限制性酶切位点的数目-Extract CLM file and counts of restriction sites #(HindIII: AAGCTT; MboI: GATC)

allhic extract sample.clean.bam draft.asm.fasta –RE AAGCTT #–RE 设置限制酶识别结合位点序列与HiC测序使用的限制酶相对应
#生成sample.clean.counts_AAGCTT.txt、sample.clean.distribution.txt、**sample.clean.clm**、sample.clean.pairs.txt

1.8 排序和方向优化

for K in seq 1 $k;
do
allhic optimize sample.clean.counts_AAGCTT.${k}g${K}.txt sample.clean.clm
done
#生成**sample.clean.counts_AAGCTT.5g5.tour、sample.clean.counts_AAGCTT.5g4.tour、sample.clean.counts_AAGCTT.5g3.tour、sample.clean.counts_AAGCTT.5g2.tour、sample.clean.counts_AAGCTT.5g1.tour**

1.9 得到最终结果

ALLHiC_build draft.asm.fasta
##生成文件
groups.asm.fasta挂载好的的染色体和一些没有挂载到染色体的contig
groups.agp挂载的日志文件,记录了挂载的contig顺序,在相邻contig之间可能加入100个N碱基相连接.

2.0 绘制heatmap图

samtools faidx groups.asm.fasta #建立索引 生成groups.asm.fasta.fai文件
cut -f 1-2 groups.asm.fasta.fai | grep “sample.clean.counts” > groups.chrn.list #获得染色体长度–切割groups.asm.fasta.fai文件中的1-2列,并搜索含有字符串“sample.clean.counts”的行,然后重定向输出到groups.chrn.list文件中.
ALLHiC_plot -b sample.clean.bam -a groups.agp -l groups.chrn.list -m 500k -o groups.hic.heatmap
#-m 500k基因组越大数值越大,每隔500k统计区域的互作信号强度.
#生成一个文件夹groups.hic.heatmap

2.1 可以使用近源物种比较挂载效果
参考QUAST中方法

#参考这里：https://github.com/tpoorten/dotPlotly/tree/master/example

ln -s $参考基因组路径 ref_genome.fa
mkdir 1.ali_ref
cd 1.ali_ref
minimap2 -t 10 -x asm5 ../ref_genome.fa ../groups.asm.fasta > overlaps.paf
Rscript $scriptsdir/pafCoordsDotPlotly.R -i overlaps.paf -o overlaps
-l -s -t -k $k -q 10000

二.3dna和Juicebox手动调整挂载结果

如果挂载结果不理想,使用juicebox手动纠错重新挂载

---

mkdir 2.juicebox
cd 2.juicebox

1.生成juicebox需要的.hic文件 group.hic、group.assembly

使用工具samtools、pairtools、matlock、3d-dna(conda安装)、juicebox(图形化界面)

#groups.agp 转换成 groups.assembly 文件

python3 $scriptsdir/agp2assembly.py ../groups.agp groups.assembly
#groups.hic 生成比较麻烦，

#Identify valid Hi-C pair

cut -f 1,2 ../draft.asm.fasta.fai > draft.asm.sizes
#生成所有contigid对应长度的文件draft.asm.sizes
samtools view -h ../sample.clean.bam | pairtools parse -c draft.asm.sizes -o parsed.pairsam.gz –cmd-out gzip
#生成文件parsed.pairsam.gz做为下一步的输入
pairtools sort –memory 10G -o sorted.pairsam.gz parsed.pairsam.gz –cmd-out gzip –cmd-in gunzip
#生成文件sorted.pairsam.gz做为下一步的输入
pairtools dedup –mark-dups -o deduped.pairsam.gz sorted.pairsam.gz –cmd-out gzip –cmd-in gunzip
#生成文件deduped.pairsam.gz做为下一步的输入
pairtools select ‘(pair_type == “UU”) or (pair_type == “UR”) or (pair_type == “RU”)’ -o filtered.pairsam.gz deduped.pairsam.gz –cmd-out gzip –cmd-in gunzip
#生成文件filtered.pairsam.gz做为下一步的输入
pairtools split –output-sam sample.filtered.bam filtered.pairsam.gz –cmd-out gzip –cmd-in gunzip
#最终生成文件sample.filtered.bam

#将sample.filtered.bam转换为long format.
#Convert .bam file to long format used by juicer (merged_nodups.txt)

matlock bam2 juicer sample.filtered.bam out.links.txt
#生成文件out.links.txt

#排序一下

sort -k2,2 -k6,6 out.links.txt > out.sorted.links.txt
#生成文件out.sorted.links.txt

#生成.hic 文件

run-assembly-visualizer.sh -p false groups.assembly out.sorted.links.txt
#生成文件groups.hic
#run-assembly-visualizer.sh 位于$miniforge3/envs/assmeble/share/3d-dna/visualize/run-assembly-visualizer.sh,可能需要给文件加执行权限.

2.使用juicebox手动调整

将groups.assembly、groups.hic下载到本地电脑.
下载安装juicebox调整组装结果.

3.使用调整后的文件groups.review.assmebly重新组装染色体.

##得到调整后的基因组

python3 $scriptsdir/juicebox_assembly_converter.py -a groups.review.assembly -f ../draft.asm.fasta -p groups.review
#输入draft.asm.fastacontig序列,将groups.review.assembly 转换成新版本的基因组文件groups.review.fasta
#还会得到.break_report.txt、.bed、.agp三个文件

##只输出调整后基因组的contig顺序(不太重要)

python3 $scriptsdir/juicebox_assembly_converter.py -a groups.review.assembly -f draft.asm.fasta -p groups.review.contig -c
#得到文件groups.review.contig.fasta
#还会得到.break_report.txt、.bed、.agp三个文件

##比对一下

minimap2 -t 10 -x asm5 ../ref_genome.fa groups.review.fasta > review.overlaps.paf
Rscript $scriptsdir/pafCoordsDotPlotly.R -i review.overlaps.paf -o review.overlaps -l -s -t -k $k -q 10000

##结果满意之后，将最终结果移入单独文件final对染色体重新命名和绘图，如果前人有研究可根据共线性结果命名染色体，如果没有研究可以按照大小排序后命名

mkdir ../final
cp groups.review.fasta ../final/groups.final.fasta
cp groups.review.agp ../final/groups.final.agp
cd ../final/

4.命名染色体并重新绘图

grep “>” groups.final.fasta #将输出结果输入到下面的echo命令中
eg: echo " PGA_scaffold_1__291_contigs__length_25707888 LG01 PGA_scaffold_2__192_contigs__length_22020340 LG02 PGA_scaffold_3__196_contigs__length_17780816 LG03 PGA_scaffold_4__162_contigs__length_14381897 LG04 PGA_scaffold_5__157_contigs__length_13401995 LG05 PGA_scaffold_6__1_contigs__length_357 Unknown1 PGA_scaffold_7__1_contigs__length_407 Unknown2 PGA_scaffold_8__1_contigs__length_414 Unknown3 PGA_scaffold_9__739_contigs__length_4945170 Unknown4 " > chrname_change_map.txt

##制作染色体名字修改对应文件

echo “
$grep输出结果
“ > chrname_change_map.txt

##修改染色体名称

python3 $scriptsdir/rename_seqID.py -r chrname_change_map.txt -f groups.final.fasta -a groups.final.agp -p genome.final
#生成genome.final.agp、genome.final.fa两个文件

##重新绘制hic热图

samtools faidx genome.final.fa
cut -f 1-2 genome.final.fa.fai > genome.final.chrn.list
ALLHiC_plot -b ../sample.clean.bam -a genome.final.agp -l genome.final.chrn.list -m 500k -o genome.final.hic.heatmap

#挂载率计算

cat genome.final.agp|grep -v “#” |awk ‘BEGIN{chrsum=0;unsum=0}{if($1!=”Unknown” && $5==”W”){chrsum=chrsum+$8-$7+1};if($1==”Unknown” && $5==”W”){unsum=unsum+$8-$7+1}}END{r=chrsum/(chrsum+unsum);print r}’

#最后比较一下出图

minimap2 -t 10 -x asm5 ref_genome.fa genome.final.fa > genome.final.overlaps.paf
Rscript $scriptsdir/pafCoordsDotPlotly.R -i genome.final.overlaps.paf -o genome.final.overlaps -l -s -t -k $k -m 10000
Rscript $scriptsdir/pafCoordsDotPlotly.R -i genome.final.overlaps.paf -o genome.final.overlaps.ordered -l -s -t -k $k -m 10000 -r 1,2,3,4,5
#生成排序好的染色体的图