3
$\begingroup$

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.

Edit: added commands. I just want to know what regions are homologous across the six genomes (99/100% IDY). I don't know what I am supposed to merge to find this out. I was reading and I believe bedtools should be but I don't know. I am usually a wetlab but I try to do a bit of bioinformatics.

 # 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
$\endgroup$
1
  • 1
    $\begingroup$ Please edit your question and show us the exact commands you used. All of them. Including how you concatenated the files. We need to know exactly what each of the columns is and then we need to see or at least understand your expected output. You only say you want to use bedtools merge, but you don't tell us what you are merging or what you want to do. When should a line be merged? When it is 100% identical? More than N% identical? Any size limits? How much identity across genomes? Please explain what you need in more detail. $\endgroup$
    – terdon
    Commented Feb 7 at 10:03

1 Answer 1

2
$\begingroup$

According to show cords the format you are getting contains columns which you can use to obtain BED format with relevant information.
You need following columns:

TAGS_genome1 
[S1]    Start of the alignment region in the reference sequence.
[E1]    End of the alignment region in the reference sequence.
TAGS_genome2
[% IDY]

I would start with not merging all coordinates but getting bed files for each pairwise genome_a__vs__genome_b.bed etc. and start from there, taking a subset of 3 genomes first (say a, b,c). You can find the regions shared by all three bedtools intersect:

bedtools intersect -a a_vs_b.bed -b a_vs_c.bed

The issue is that you will have many to many relationships, since region_1 from b can have several similar regions in b, and each one of these (== regions in b) also can be similar to a number of regions in a. Depending on what you are looking for and how dis-similar the genomes are you may end up with reporting repetitive elements or if you filter hard missing i.e. exons.

$\endgroup$
1
  • 2
    $\begingroup$ Thanks @darked89. That's exactly what I needed, those 5 columns. I didn't consider that there are too many relationships. I will then continue to do comparisons. I think will give me the answer that I want. bedtools intersect -a a_vs_b.bed -b a_vs_c.bed > a_vs_b_vs_c.bed bedtools intersect -a a_vs_b_vs_c.bed vs a_vs_d_vs_e.bed > all_compared.bed $\endgroup$
    – LORL
    Commented Feb 9 at 2:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.