下載mummer缎岗,并添加到環(huán)境變量
wget https://gigenet.dl.sourceforge.net/project/mummer/mummer/3.23/MUMmer3.23.tar.gz
tar-xfMUMmer3.23.tar.gz
cdMUMmer3.23
make install
vim ~/.bashrc
export PATH=$PATH:/home/liuzhh/software/mummer-4.0.0rc1/bin
source ~/.bashrc
執(zhí)行兩個基因組間的比對,并篩選結(jié)果
nucmer?/home/liuzhh/Project/01.Analysis/7hicGenome/10.Genome_Trim/G07802.genome.fasta??/ifs1/liuzhh/01.ChiZhi/00.Denovo/0406_110samples/06.bwa_freebayes-polish/G07802/G07802.contig
只保留1對1的最佳聯(lián)配
delta-filter -1 out.delta>filter.delta
輸出容易展示的信息
show-coords -T -q -H filter.delta>coord.txt
將比對結(jié)果可視化
用mummerplot可視化
下載:https://jpegclub.org/support/files/jpegsrc.v8d1.tar.gz
./configure --prefix=/home/chensiyuan/local
make
make install
conda install gnuplot=5.4
/home/chensiyuan/mummer-4.0.0/bin/mummerplot? -png? /home/chensiyuan/13.mummer/G07802/out.delta? ? --large
用R可視化
# 確保已經(jīng)加載了所需的庫
library(ggplot2)
library(dplyr)
library(scales)
file <- "E:/"
nucmer <- read.delim(paste0(file, "coord.txt"),
? ? ? ? ? ? ? ? ? ? header = F,
? ? ? ? ? ? ? ? ? ? sep = "\t",
? ? ? ? ? ? ? ? ? ? stringsAsFactors = FALSE)
colnames(nucmer) <- c("ref_start", "ref_end", "qry_start", "qry_end","ref_len", "qry_len", "identiy", "ref_tag","qry_tag")
#chr1_nuc <-? subset(nucmer,? ref_tag == 'chr1' & qry_tag == 'chr1')
# 定義染色體的順序,從 Chr15 到 Chr1
chromosome_order <- rev(c("Chr1", "Chr2", "Chr3", "Chr4", "Chr5", "Chr6","Chr7", "Chr8", "Chr9", "Chr10", "Chr11", "Chr12", "Chr13","Chr14", "Chr15"))
# 將 nucmer 數(shù)據(jù)框中的 ref_tag 和 qry_tag 轉(zhuǎn)換為因子镰官,并設(shè)置水平順序
nucmer$ref_tag <- factor(nucmer$ref_tag, levels = c("Chr1", "Chr2", "Chr3", "Chr4", "Chr5", "Chr6","Chr7", "Chr8", "Chr9", "Chr10", "Chr11", "Chr12", "Chr13", "Chr14", "Chr15"))
nucmer$qry_tag <- factor(nucmer$qry_tag, levels = chromosome_order)
# 使用 ggplot 創(chuàng)建圖形
p=ggplot(nucmer, aes(x = ref_start, y = qry_start)) +
? geom_segment(aes(xend = ref_end, yend = qry_end), size = 1) +
? facet_grid(qry_tag ~ ref_tag, switch = "both", scales = "free") +
? theme_bw() +
? theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 4),
? ? ? ? axis.text.y = element_text(size = 4, angle = 0, vjust = 0.5),? # 確保 Y 軸標(biāo)簽垂直顯示
? ? ? ? strip.text.x = element_text(size = 6),
? ? ? ? strip.text.y = element_text(size = 5)) +
? xlab("HI-C Genome ") +
? ylab("NO HI-C Genome")
print(p)