Chromosome-Separated Analysis Results

[CJHtea]

Dear @bxuecarnegie 😊,
Thank you for developing such a powerful tool​ for gene cluster detection!
I have a question about the results. To save runtime, I separated the genome.fa by chromosome for analysis (-dnaf "chr.fa").
However, the output curiously includes gene cluster information from other chromosomes.
Why might this be the case?
Looking forward to your response.
Best,
Jiahao Chen

[CJHtea]

On the other hand, when I tried to split the chromosomes for processing, I found that those assembled to a "no gaps" level would report errors, causing the subsequent steps to fail to run properly:

[bxuecarnegie]

Hi,

For your first question, I'm assuming the gene cluster information from other chromosomes were retrieved from your -glof file.
For the second question, as I'm not the main author of this script, I don't know if it were written to run with no gaps level assembled genome. You can test the "awk" commands listed in the log file and check the output, if the awk scripts were able to run correctly but produces empty results then it probably won't run with this pipeline.

Best,
Bo

[CJHtea]

Dear @bxuecarnegie,
Thanks for your quick response!
I have manually executed each command for debugging and confirmed that if the chromosome contains no Ns, the result will not be output correctly. Therefore, I modified the get_positions_of_gap.awk script. The following script can now correctly count chromosomes without gaps:

{
  # Takes a FASTA file with header lines (>) followed by a single-line
  # sequence per chromosome/contig.
  # Input: FASTA from stdin or file
  # Output: contig, type (N or A), start, length

    # Non-header lines: process sequence
    if ( $0 !~ /^>/ ) {

                # Normalize bases:
                #   - All A/T/C/G (upper/lower) and '*' -> 'A'
                #   - All n/N -> 'N'
                seq = $0
                gsub(/[aAtTcCgG*]/, "A", seq)
                gsub(/[nN]/, "N", seq)

                len = length(seq)

                # Arrays to store N segments and A segments for this sequence
                n_cnt = 0
                a_cnt = 0

                # Single pass over the sequence to find contiguous N and A runs
                i = 1
                while (i <= len) {
                        base = substr(seq, i, 1)

                        if (base == "N") {
                                # Scan one contiguous N run
                                start = i
                                while (i <= len && substr(seq, i, 1) == "N") {
                                        i++
                                }
                                seglen = i - start

                                # Store this N run
                                n_cnt++
                                n_start[n_cnt] = start
                                n_len[n_cnt]   = seglen

                        } else if (base == "A") {
                                # Scan one contiguous A run (non-N bases)
                                start = i
                                while (i <= len && substr(seq, i, 1) == "A") {
                                        i++
                                }
                                seglen = i - start

                                # Store this A run
                                a_cnt++
                                a_start[a_cnt] = start
                                a_len[a_cnt]   = seglen

                        } else {
                                # Should not occur after normalization; skip defensively
                                i++
                        }
                }

                # First output all N segments (same order as original N-pass)
                for (k = 1; k <= n_cnt; k++) {
                        print vheader[1], "N", n_start[k], n_len[k]
                }

                # Then output all A segments (same order as original A-pass)
                for (k = 1; k <= a_cnt; k++) {
                        print vheader[1], "A", a_start[k], a_len[k]
                }

    } else {
                # Header line: store contig/chromosome ID in vheader[1]
                split($0, vheader, " ")
                gsub(/>/, "", vheader[1])
    }
}

I will continue to monitor whether the modified script can output results smoothly and provide you with feedback.
Best,
Jiahao Chen

[CJHtea]

Another observation: the following is the log file from one of my runs:

My questions are:

1. Does 6587 represent the total number of gene clusters? But my output results only have 361 rows.
2. The log indicates that it seems to be writing out the results, but the software run has ended without outputting any error or warning messages. Why is this? Although the same phenomenon occurs with the test data, the output results match the expected results.

[CJHtea]

Dear @bxuecarnegie, Thanks for your quick response! I have manually executed each command for debugging and confirmed that if the chromosome contains no Ns, the result will not be output correctly. Therefore, I modified the get_positions_of_gap.awk script. The following script can now correctly count chromosomes without gaps:

{
  # Takes a FASTA file with header lines (>) followed by a single-line
  # sequence per chromosome/contig.
  # Input: FASTA from stdin or file
  # Output: contig, type (N or A), start, length

    # Non-header lines: process sequence
    if ( $0 !~ /^>/ ) {

                # Normalize bases:
                #   - All A/T/C/G (upper/lower) and '*' -> 'A'
                #   - All n/N -> 'N'
                seq = $0
                gsub(/[aAtTcCgG*]/, "A", seq)
                gsub(/[nN]/, "N", seq)

                len = length(seq)

                # Arrays to store N segments and A segments for this sequence
                n_cnt = 0
                a_cnt = 0

                # Single pass over the sequence to find contiguous N and A runs
                i = 1
                while (i <= len) {
                        base = substr(seq, i, 1)

                        if (base == "N") {
                                # Scan one contiguous N run
                                start = i
                                while (i <= len && substr(seq, i, 1) == "N") {
                                        i++
                                }
                                seglen = i - start

                                # Store this N run
                                n_cnt++
                                n_start[n_cnt] = start
                                n_len[n_cnt]   = seglen

                        } else if (base == "A") {
                                # Scan one contiguous A run (non-N bases)
                                start = i
                                while (i <= len && substr(seq, i, 1) == "A") {
                                        i++
                                }
                                seglen = i - start

                                # Store this A run
                                a_cnt++
                                a_start[a_cnt] = start
                                a_len[a_cnt]   = seglen

                        } else {
                                # Should not occur after normalization; skip defensively
                                i++
                        }
                }

                # First output all N segments (same order as original N-pass)
                for (k = 1; k <= n_cnt; k++) {
                        print vheader[1], "N", n_start[k], n_len[k]
                }

                # Then output all A segments (same order as original A-pass)
                for (k = 1; k <= a_cnt; k++) {
                        print vheader[1], "A", a_start[k], a_len[k]
                }

    } else {
                # Header line: store contig/chromosome ID in vheader[1]
                split($0, vheader, " ")
                gsub(/>/, "", vheader[1])
    }
}

I will continue to monitor whether the modified script can output results smoothly and provide you with feedback. Best, Jiahao Chen

I tested using a chromosome without gaps, and it output normally.😊

[bxuecarnegie]

Thank you for the update and testing the script. It produces the same result with the example input so I have updated the repo. I have added you to the acknowledgement section of the README.

Best,
Bo

[CJHtea]

Thank you for adopting my suggestion and including me in the acknowledgment.
I hope this will benefit more scholars focusing on plant gene clusters.
Meanwhile, when counting "N" and "A" in larger genomes, this step runs very slowly, likely because the original pipeline uses an awk-based strategy. Without altering the original workflow, I recommend:

1. Split the genome by chromosomes.
2. Manually perform "N" and "A" counting for each chromosome according to the pipeline steps.
3. Summarize and name the counting results according to the required input format, then place them in the folder where the pipeline originally outputs summarized results.
4. When performing formal analysis, create ghost genome.fa, genome.fa_temp1, genome.fa_temp2, genome.fa_temp3, and genome.fa_temp4 files.

I believe this approach will save a significant amount of time.

Best,
Jiahao Chen

[CJHtea]

If I’m only interested in terpene-related gene clusters, how should I set the -md?

[bxuecarnegie]

-md {'Terpenoids'} should limit the analyzed metabolic domain to terpenoids.