Clean to remove contigs with only 3 or less nucleotides

[BenjaminSchwessinger]

Hi Jon,

It would be great if clean had an option to remove contigs that have less than 4 nucleotide types. Predict raises this as error and it would be great to simply have an option to clean this up beforehand.

Thanks!

[nextgenusfs]

Yes that's a good suggestion I will add to do list.

[hyphaltip]

I might add this to AAFTF too - it is a good point. @stajichlab/AAFTF

[calliphriend]

Hello all. Apologies in advance, I am very new to this. I am working on a polychaete genome in 1.8.16 (Augustus would not work in 1.8.3). When I run clean, I still get the "contains only 3 nucleotides, removing" error numerous times until the process stops running (slurm log file remains same size) without completing/producing the output fasta file. I am definitely concerned the genome is too large to run clean as suggested in other threads, but even with increasing minlength significantly the job doesn't finish. I can't run predict without using force, which I don't want to do with so many contigs that need to be removed. I would sincerely appreciate any direction, thank you.

[nextgenusfs]

Hi @calliphriend, no need to apologize!

The clean functionality was written to clean up assembly artifacts from haploid assemblies -- basically it tries to remove contigs if they are fully contained within a larger contig. This might not be a safe assumption for all organisms. You don't need to use it. But it would probably help to clean up some of the contigs. There are likely many other ways to do this.

How large is the genome? What is the ploidy of organism and ploidy of assembly?

[calliphriend]

Thank you kindly for your help and very quick response! This is excellent to know. I have been looking for alternatives to clean up the only 3 nucleotide contigs and BLASTN or DeconSeq seems promising.

The assembly size is 1,509,999,849 bp; the genome assembly haploid length is 366,115,238 bp. The ploidy of the organism is unknown unfortunately. Thank you again.

[nextgenusfs]

Another FYI, I'm not really sure how well funannotate will perform assemblies where ploidy is > 1. I mostly work on fungi (haploid) so haven't experimented with the best approaches to annotating ploidy > 1. One approach might be to annotate the best haplotigs (ie haploid assembly) and then lift those gene models over to the rest of the assembly if that is necessary.

Would be easy to write a short little python script to remove/scrubb contigs that are either short and/or are likely "bad", ie do not contain all 4 nucleotides.

[calliphriend]

For sure – thank you so much for these excellent suggestions! I'll write that script and see if it solves the issue. Many thanks again.

[nextgenusfs]

I quite like the pyfastx library, its fast and you could do a length and composition quickly like this:

>>> import pyfastx
>>> fa = pyfastx.Fasta('GCA_030734395.1_ASM3073439v1_genomic.fna', uppercase=True)
>>> for contig in fa:
...     print(contig.name, len(contig), contig.composition)
... 
JAUZTG010000001.1 117374 {'A': 27159, 'C': 31240, 'G': 30389, 'T': 28586}
JAUZTG010000010.1 17123 {'A': 4082, 'C': 4605, 'G': 4366, 'T': 4070}
JAUZTG010000011.1 631397 {'A': 149375, 'C': 166937, 'G': 165844, 'T': 149241}
JAUZTG010000012.1 1356386 {'A': 321292, 'C': 357175, 'G': 357370, 'T': 320549}
JAUZTG010000013.1 369857 {'A': 87558, 'C': 98041, 'G': 98746, 'T': 85512}
JAUZTG010000014.1 323906 {'A': 75668, 'C': 87788, 'G': 84413, 'T': 76037}
JAUZTG010000015.1 711434 {'A': 169210, 'C': 188405, 'G': 184840, 'T': 168979}
JAUZTG010000016.1 2337667 {'A': 552855, 'C': 615380, 'G': 620778, 'T': 548654}
JAUZTG010000017.1 1183939 {'A': 279705, 'C': 307699, 'G': 315515, 'T': 281020}
JAUZTG010000018.1 15334 {'A': 3562, 'C': 4095, 'G': 4264, 'T': 3413}
JAUZTG010000019.1 255293 {'A': 60553, 'C': 67736, 'G': 66704, 'T': 60300}
JAUZTG010000002.1 171319 {'A': 39829, 'C': 45522, 'G': 45688, 'T': 40280}
JAUZTG010000020.1 2236409 {'A': 521766, 'C': 596155, 'G': 596289, 'T': 522199}
JAUZTG010000021.1 207546 {'A': 50416, 'C': 53556, 'G': 54198, 'T': 49376}
JAUZTG010000022.1 19207 {'A': 4840, 'C': 4608, 'G': 4803, 'T': 4956}
JAUZTG010000023.1 642669 {'A': 152032, 'C': 170557, 'G': 167617, 'T': 152463}
JAUZTG010000024.1 370210 {'A': 86522, 'C': 97253, 'G': 99196, 'T': 87239}
JAUZTG010000025.1 86159 {'A': 20657, 'C': 22134, 'G': 23278, 'T': 20090}
JAUZTG010000026.1 1082 {'A': 267, 'C': 259, 'G': 290, 'T': 266}
JAUZTG010000027.1 148686 {'A': 37995, 'C': 36448, 'G': 36822, 'T': 37421}
JAUZTG010000028.1 389535 {'A': 92603, 'C': 102837, 'G': 101755, 'T': 92340}
JAUZTG010000029.1 7198218 {'A': 1693551, 'C': 1900575, 'G': 1904200, 'T': 1699892}
JAUZTG010000003.1 1272425 {'A': 298686, 'C': 339195, 'G': 338119, 'T': 296425}
JAUZTG010000030.1 55508 {'A': 13099, 'C': 14443, 'G': 14471, 'T': 13495}
JAUZTG010000031.1 65210 {'A': 15098, 'C': 18454, 'G': 16620, 'T': 15038}
JAUZTG010000032.1 134135 {'A': 30961, 'C': 36127, 'G': 35937, 'T': 31110}
JAUZTG010000033.1 144693 {'A': 35814, 'C': 37012, 'G': 35907, 'T': 35960}
JAUZTG010000034.1 2106387 {'A': 494349, 'C': 557823, 'G': 559727, 'T': 494488}
JAUZTG010000035.1 109034 {'A': 25628, 'C': 28663, 'G': 28557, 'T': 26186}
JAUZTG010000036.1 851245 {'A': 202406, 'C': 225009, 'G': 223196, 'T': 200634}
JAUZTG010000037.1 55100 {'A': 13288, 'C': 14354, 'G': 14322, 'T': 13136}
JAUZTG010000038.1 10593 {'A': 2585, 'C': 2723, 'G': 2851, 'T': 2434}
JAUZTG010000039.1 634250 {'A': 151350, 'C': 167371, 'G': 164606, 'T': 150923}
JAUZTG010000004.1 462262 {'A': 107929, 'C': 122828, 'G': 122638, 'T': 108867}
JAUZTG010000040.1 56101 {'A': 13855, 'C': 14379, 'G': 14100, 'T': 13767}
JAUZTG010000041.1 2650198 {'A': 620314, 'C': 704832, 'G': 701284, 'T': 623768}
JAUZTG010000042.1 2899401 {'A': 680290, 'C': 769683, 'G': 770009, 'T': 679419}
JAUZTG010000043.1 31136 {'A': 7775, 'C': 7790, 'G': 7997, 'T': 7574}
JAUZTG010000044.1 1766052 {'A': 412071, 'C': 467727, 'G': 469286, 'T': 416968}
JAUZTG010000045.1 10721 {'A': 2687, 'C': 2841, 'G': 2921, 'T': 2272}
JAUZTG010000046.1 39834 {'A': 10483, 'C': 9454, 'G': 9538, 'T': 10359}
JAUZTG010000047.1 662718 {'A': 157206, 'C': 176761, 'G': 171759, 'T': 156992}
JAUZTG010000048.1 3045196 {'A': 713334, 'C': 811023, 'G': 807608, 'T': 713231}
JAUZTG010000049.1 13925 {'A': 3190, 'C': 3934, 'G': 3600, 'T': 3201}
JAUZTG010000005.1 29108 {'A': 6845, 'C': 7591, 'G': 7602, 'T': 7070}
JAUZTG010000050.1 247851 {'A': 59902, 'C': 64658, 'G': 64739, 'T': 58552}
JAUZTG010000051.1 449909 {'A': 105568, 'C': 117277, 'G': 118821, 'T': 108243}
JAUZTG010000006.1 5416 {'A': 1407, 'C': 1267, 'G': 1325, 'T': 1417}
JAUZTG010000007.1 527676 {'A': 123720, 'C': 139326, 'G': 141054, 'T': 123576}
JAUZTG010000008.1 470228 {'A': 109072, 'C': 125266, 'G': 126100, 'T': 109790}
JAUZTG010000009.1 28168 {'A': 7896, 'C': 6221, 'G': 6272, 'T': 7779}
>>> 

Trivial to then just remove contigs that are too short or have less than 4 nucleotides, ie len(contig.composition) in the above example. passing uppercase=True will force all uppercase, otherwise you can have both lowercase and uppercase nucleotides in the composition dictionary.

https://pyfastx.readthedocs.io/en/latest/usage.html#fasta