Right now I have six genomes that I want to compare and identify homologous regions in the genomes. I have run nucmer, show-coords and obtained the output files. An example is shown below with Genome 1 vs Genome 2. More files go through Genome 1 vs Genome 3 until all are compared.
I was trying to use bedtools merge after concatenating the multiple output files from show-coords command to find the homologous regions. However, I am not quite sure on the command to use. Any help is appreciated.
# Loop through all pairs of genomes
for i in {1..6}; do
for j in $(seq $((i + 1)) 6); do
genome1="genome${i}.fasta"
genome2="genome${j}.fasta"
output_prefix="output_${i}_${j}"
# Run NUCmer
nucmer -p "${output_prefix}" "$genome1" "$genome2"
# Show coordinates
show-coords -r -l -c -T "${output_prefix}.delta" > "${output_prefix}_coords.txt"
done
done
After using the above I then concatenate the files with
cat output_*_coords.txt > combined_coords.txt
This gives me something that looks like this below.
[S1] [E1] [S2] [E2] [LEN 1] [LEN 2] [% IDY] [LEN R] [LEN Q] [COV R] [COV Q] [TAGS]
1 14347 1 14347 14347 14347 99.02 5064019 5004926 0.28 0.29 Genome1 Genome2
14474 32175 14350 32051 17702 17702 99.33 5064019 5004926 0.35 0.35 Genome1 Genome2
16541 22042 4698474 4692975 5502 5500 99.64 5064019 5004926 0.11 0.11 Genome1 Genome2
16548 22042 3728461 3722968 5495 5494 99.84 5064019 5004926 0.11 0.11 Genome1 Genome2
16583 22041 3904238 3898781 5459 5458 99.65 5064019 5004926 0.11 0.11 Genome1 Genome2
18727 22097 4426876 4423506 3371 3371 99.73 5064019 5004926 0.07 0.07 Genome1 Genome2
