Can "polish" gaps"?

[desmodus1984]

Hi,
I am interested in improving a genome assembly, and found your paper: "GoldPolish-target: targeted long-read genome assembly polishing".
It was very hard to find this github, since it is not mentioned in the paper.
I did check the "quality" of the genome and I used Hapog - since my genome is pseudohaploid, and based on BUSCO and proteins with errors, it didn't improve much. We did a first assembly using just PacBio HiFi reads and then it was scaffolded using Hi-C. The k-mer analysis was shocking because, despite better contiguity, gaps were more and longer, but the QV value was exactly the same.
Thus, I am interested now in "correcting" or "filling" the gaps.
I got very frustrated because according to an analysis, gap-filling can mess up with assembly, and I tried Samba from MaSurCa and it failed.
Amazingly, while looking for software I would GoldPolish-target, and I would like to try it. I was able to get the locations in bed for the gaps, but I wanted to consult you beforehand, if it is fine to pressume that it will "fill" the gaps. For instance, I asked the authors of Hapog if it will fill the gaps, and the answer was negative.

Thus, I wanted to ask you before if my pressumption/goal is actually correct/possible.

Thank you.

[emilyyzhangg]

Hi,

Thank you for reaching out! If you have large gaps comprised of Ns, GoldPolish-target alone likely will not work. Instead, I would suggest running our scaffolder tool, ntLink (https://github.com/bcgsc/ntLink) with the gap-fill feature, to first fill the gaps. Then you could run GoldPolish-target with the BED file that you've created to polish the filled gaps. This may help to improve the quality of the genome.

Please let me know if you have any additional questions or concerns.

Best,

Emily

[desmodus1984]

Hi,
Thanks for the response.
May I ask what do you mean with "large gaps"?
I used gfastats to get gaps information and maybe the gaps aren't "large" as you might think. There are 33 gaps, with mean gap length of 90.67, largest gap of 100, and smallest of 23. I was thinking of extending the gap start/end perhaps like 500 or 1000 bp, and give this new bed location for GoldPolish-target to polish, since the reads are long (PacBio Hifi).
But if you think that ntLink is better I will use this instead.
Should I provide the fastq.gz files or fasta input? Is it beneficial to filter out reads shorter than 5 - 10k?
Thanks.

[emilyyzhangg]

Hi,

Unfortunately, the k-mer size that we use is smaller than the gap size that you've indicated so you'll want to use ntLink to gap-fill before polishing with GoldPolish-target. I'd recommend that you use the two tools in conjunction with each other. You should use fasta files.

Best,

Emily

[desmodus1984]

Hi,
I ran ntlink and I ran into an error.
My script is this:

# Init conda
source /data/common/apps/conda_init.sh
which conda

#activate conda environment
conda activate synteny

# Save starting time
start_time=$(date +%s)
echo "HL-gaps-ntlink: Job started: $(date)"

#Change directory
cd  /data/common/data/HLep/

#HL gap filling?

ntLink scaffold gap_fill t=64 target=final_assembly.fasta reads=Hlep-PacBio.fasta k=32 w=250


# Save ending time
end_time=$(date +%s)
echo "HL-ntlink-gaps: Job completed: $(date)"

and the error message is this:

indexlr 1.7.5: Using more than 5 threads does not scale, reverting to 5.
indexlr 1.7.5: Using more than 5 threads does not scale, reverting to 5.
indexlr 1.7.5: Using more than 5 threads does not scale, reverting to 5.
The minimum coverage of single-end contigs is inf.
The minimum coverage of merged contigs is inf.
Traceback (most recent call last):
  File "/data/common/apps/miniconda3/envs/synteny/bin/share/ntlink-1.3.11-1/bin/ntlink_patch_gaps.py", line 840, in <module>
    main()
  File "/data/common/apps/miniconda3/envs/synteny/bin/share/ntlink-1.3.11-1/bin/ntlink_patch_gaps.py", line 822, in main
    map_long_reads(pairs, sequences, args)
  File "/data/common/apps/miniconda3/envs/synteny/bin/share/ntlink-1.3.11-1/bin/ntlink_patch_gaps.py", line 428, in map_long_reads
    assert source_id == source and label == "source"
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AssertionError
make: *** [/data/common/apps/miniconda3/envs/synteny/bin/share/ntlink-1.3.11-1/ntLink:267: final_assembly.fasta.k32.w250.z1000.ntLink.scaffolds.gap_fill.fa] Error 1

Is there an error on the program or on the input files?
Is ran fine for a while but then it failed.

Thanks;

[emilyyzhangg]

Hi! It may be an issue with the naming of the sequences. Could give us an example of a couple of sequence headers in your assembly file? We also want to confirm that Hlep-PacBio.fasta is your raw reads file. Further, can you give us some examples of some of these read names?

[emilyyzhangg]

If you'd like, you could try renaming all your sequences with this command seqtk rename file.fa seq_ and see if it fixes the problem!

[desmodus1984]

Hi! It may be an issue with the naming of the sequences. Could give us an example of a couple of sequence headers in your assembly file? We also want to confirm that Hlep-PacBio.fasta is your raw reads file. Further, can you give us some examples of some of these read names?

Hi @emilyyzhangg

Thanks and sorry. I just found out about this reply. I don't know why I didn't receive a notification about it.
I checked the headers and they look like this:

m84050_231222_221154_s3/31526375/ccs
m84050_231222_221154_s3/15140379/ccs
m84050_231222_221154_s3/47187042/ccs
m84050_231222_221154_s3/102241236/ccs
m84050_231222_221154_s3/214437175/ccs
m84050_231222_221154_s3/99747610/ccs
m84050_231222_221154_s3/94836224/ccs
m84050_231222_221154_s3/237111186/ccs
m84050_231223_001121_s4/128126969/ccs
m84050_231222_221154_s3/22876723/ccs
m84050_231222_221154_s3/42470631/ccs

Should I rename it to something like Hlep-PacBio-{nr}?

Thanks