mummer进行基因组间的比对

2024-09-08  本文已影响0人  小陈生信日记

下载mummer,并添加到环境变量

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

执行两个基因组间的比对,并筛选结果

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的最佳联配

delta-filter -1 out.delta>filter.delta

输出容易展示的信息

show-coords -T -q -H filter.delta>coord.txt

将比对结果可视化

用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可视化

# 确保已经加载了所需的库

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 数据框中的 ref_tag 和 qry_tag 转换为因子,并设置水平顺序

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 创建图形

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 轴标签垂直显示

        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)

上一篇 下一篇

猜你喜欢

热点阅读