<p>You could use <code>numpy.random.randint</code> to sample integers from a uniform distribution, from 1000 to 15000 nt, extracting subsequences of that length from a linearized multi-FASTA file.</p>
<p>Here's some code:</p>
<pre><code>#!/usr/bin/env python

import sys
import numpy as np

low = int(sys.argv[1])
high = int(sys.argv[2])

def process_record(header, sequence, low, high):
    next = 0
    fragment = 1
    while next &lt; len(sequence):
        new_next = np.random.randint(low, high) + next  
        if new_next &gt; len(sequence) - high:
            new_next = len(sequence)
        sys.stdout.write('{}:{}\n{}\n'.format(header, fragment, sequence[next:new_next]))
        next = new_next
        fragment += 1

header = &quot;&quot;
sequence = &quot;&quot;
for line in sys.stdin:
    if line.startswith('&gt;'):
        if len(sequence) &gt; 0:
            process_record(header, sequence, low, high)
            sequence = &quot;&quot;
        header = line.rstrip()
    else:
        sequence += line.rstrip()

if len(sequence) &gt; 0:
    process_record(header, sequence, low, high)
</code></pre>
<p>To use it:</p>
<pre><code>$ ./13879.py 1000 15000 &lt; in.fa &gt; out.fa
</code></pre>
<p>Fragments will be disjoint (non-overlapping, contiguous) subsequences.</p>
<p>I add a check to determine if the last fragment of a FASTA record would be shorter than <code>low</code> (1kb) in length, and, if so, append that fragment to the final subsequence. It is possible for that last fragment to be up to 16 kb - 1 nt in length, given your parameters:</p>
<pre><code>if new_next &gt; len(sequence) - high:
    new_next = len(sequence)
</code></pre>
<p>You could change this logic to discard the last fragment (i.e. <code>break</code>), or simply allow fragments less than 1kb in size, depending on what you need.</p>
<p><strong>EDIT</strong></p>
<p>Here's a slight modification that should allow to you preserve fragment sizes within your bounds.</p>
<p>A cartoon may help illustrate the problem. Let's set up a FASTA record, which has been broken up into <code>n</code> fragments:</p>
<pre><code>   1       2           n-2      n-1      n
|-----|---------|...|-------|---------|-----|
</code></pre>
<p>We may get unlucky and sample a random integer for fragment <code>n-1</code> (<code>frag(n-1)</code>) which makes the length of <code>frag(n)</code> less than <code>low</code> (1 kb in your example).</p>
<p>In my first suggested solution, <code>frag(n-1)</code> and <code>frag(n)</code> are appended to one another. This is why you get some fragments longer than <code>high</code> (15 kb).</p>
<p>Another option is to discard <code>frag(n)</code>, which is of length less than <code>low</code> (1 kb). Here's some code:</p>
<pre><code>#!/usr/bin/env python

import sys
import numpy as np
import timeit

low = int(sys.argv[1])
high = int(sys.argv[2])
in_fn = sys.argv[3]

&quot;&quot;&quot;
Globals
&quot;&quot;&quot;
fragment_lengths = []

&quot;&quot;&quot;
Process a record
&quot;&quot;&quot;
def process_record(header, sequence, low, high):
    next = 0
    last = 0
    fragment = 1
    fragment_sequence = &quot;&quot;
    while last &lt; len(sequence) - high:
        next = np.random.randint(low, high + 1) + last
        fragment_sequence += sequence[last:next]
        # Comment the sys.stdout.write and uncomment the fragment_lengths line, if measuring fragment length distributions 
        sys.stdout.write('{}:{}-{}:{}\n{}\n'.format(header, last, next, fragment, sequence[last:next]))
        #fragment_lengths.append(next - last)                                                                                                                                                                                                                                                                                                                                                                                          
        last = next
        fragment += 1
    &quot;&quot;&quot;
    Uncomment to compare input sequence with fragments
    &quot;&quot;&quot;
    #sys.stderr.write('In &gt; {}\nOut &gt; {}\n'.format(sequence, fragment_sequence))
    
&quot;&quot;&quot;
Parse records
&quot;&quot;&quot;
def parse_records(in_fn, low, high):
    header = &quot;&quot;
    sequence = &quot;&quot;
    with open(in_fn, 'r') as ifh:
        for line in ifh:
            if line.startswith('&gt;'):
                if len(sequence) &gt; 0:
                    process_record(header, sequence, low, high)
                    sequence = &quot;&quot;
                header = line.rstrip()
            else:
                sequence += line.rstrip()
        if len(sequence) &gt; 0:
            process_record(header, sequence, low, high)

def count_elements(seq):
    hist = {}
    for i in seq:
        hist[i] = hist.get(i, 0) + 1
    return hist

def wrapper(func, *args, **kwargs):
    def wrapped():
        return func(*args, **kwargs)
    return wrapped

def main(in_fn, low, high):
    &quot;&quot;&quot;
    Comment below if testing the distribution of fragments
    &quot;&quot;&quot;
    parse_records(in_fn, low, high)
    &quot;&quot;&quot;
    Uncomment to test the distribution of fragments from 
    running fns on same input 100k times
    &quot;&quot;&quot;
    #wrapped = wrapper(parse_records, in_fn, low, high)
    #timeit.timeit(wrapped, number=100000)
    #print(count_elements(fragment_lengths))

if __name__ == &quot;__main__&quot;:
    main(in_fn, low, high)
</code></pre>
<p>If <code>frag(n)</code> is needed, and you need to preserve the uniform distribution (or whatever fragment length distribution you want to model), then <a href="https://en.wikipedia.org/wiki/Rejection_sampling" rel="nofollow noreferrer">rejection sampling</a> might be an option.</p>
<p>With rejection sampling, you keep fragmenting your input sequence, until one of the distributions of lengths gets you a set of fragments that entirely cover the input sequence.</p>
<p>If you need that, follow up and I'll take another look.</p>

CC BY-SA 4.0
2020-07-14T16:26:51.910Z
Thanks for your answer! Do I need to remove all the fasta headers? I tried it on a multifasta and ended up with a few sequences over 15 kb. The distribution looks like this: https://imgur.com/a/FFr8xNv
Stanley Ho
CC BY-SA 4.0
2020-07-14T20:10:38.237Z
No, you do not need to remove headers. You could end up with a few sequences over 15kb, if the last fragment cannot be sized between `low` and `high` (i.e. between 1 and 15kb). I'll propose a solution in the edit.
Alex Reynolds
CC BY-SA 4.0
2020-07-15T10:47:27.460Z
Thanks that worked perfectly!
Stanley Ho
CC BY-SA 4.0
2020-07-15T10:57:53.693Z
Hi @StanleyHo please do sign the answer as "asscepted". I'll be making use of this code too, I like the args and kwargs
M__
CC BY-SA 4.0
2020-07-13T19:33:37.747Z
2020-07-15T18:58:53.347Z
2020-07-15T18:58:53.347Z
Alex Reynolds
<p>I've previously written <a href="https://gitlab.com/gringer/bioinfscripts/-/blob/master/normalise_seqlengths.pl" rel="nofollow noreferrer">my own script</a> to fragment DNA sequences (input as either fastq or fasta files) into identical-length fragments, with a given overlap. I used this because there was a particular assembler that demanded equal-length sequences, and I was working with nanopore reads, but it's since been useful for other purposes as well.</p>
<p>Because it is most concerned about making lengths identical, there may be some forced overlap when setting the overlap to zero. The overlap will only be exactly zero when the sequence length is an exact multiple of the fragment size. Here's an example use:</p>
<h2>zero overlap</h2>
<pre><code>normalise_seqlengths.pl -fraglength 25 -overlap 0 oneRead.fq
@A00187:180:HLGNGDSXX:3:1101:8486:1031#0 [150 bp, 6 fragments, 0 overlap] 1:N:0:CGCAAGAA+ACTCCTAA
ATTGCGCAATGACAGCCAGGGGAGT
+
FFFFFFFFFFFFFFFFF:,FFFFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#1:25-50  1:N:0:CGCAAGAA+ACTCCTAA
GGGAGCTGTCTGGGGTAACCTGGCT
+
:FFFF:F:FFFFF::FFFFF,::FF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#2:50-75  1:N:0:CGCAAGAA+ACTCCTAA
CTAGGACCAGGTGGAAAAAGCACCT
+
F:FFFFFFFFFFFFFFFFFFFFFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#3:75-100  1:N:0:CGCAAGAA+ACTCCTAA
TGGCTGGCCCAGCCCTGTTCTTCTC
+
:FF,FFFFF,:,FFF::F:FFF,:F
@A00187:180:HLGNGDSXX:3:1101:8486:1031#4:100-125  1:N:0:CGCAAGAA+ACTCCTAA
CTGGGGAGGGAAGGGCAGGTCCCCT
+
F:FFFFFFFFF,FFFFF:FFF,F,,
@A00187:180:HLGNGDSXX:3:1101:8486:1031#5:125-150  1:N:0:CGCAAGAA+ACTCCTAA
AATGCCTGAACCACATGGTGGGCCA
+
FFF,F:FFFFFFF:F,,,FFF,F,,
</code></pre>
<h2>partial overlap</h2>
<pre><code>normalise_seqlengths.pl -fraglength 63 -overlap 0 oneRead.fq
@A00187:180:HLGNGDSXX:3:1101:8486:1031#0 [150 bp, 3 fragments, 20 overlap] 1:N:0:CGCAAGAA+ACTCCTAA
ATTGCGCAATGACAGCCAGGGGAGTGGGAGCTGTCTGGGGTAACCTGGCTCTAGGACCAGGTG
+
FFFFFFFFFFFFFFFFF:,FFFFFF:FFFF:F:FFFFF::FFFFF,::FFF:FFFFFFFFFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#1:43-106  1:N:0:CGCAAGAA+ACTCCTAA
CCTGGCTCTAGGACCAGGTGGAAAAAGCACCTTGGCTGGCCCAGCCCTGTTCTTCTCCTGGGG
+
FF,::FFF:FFFFFFFFFFFFFFFFFFFFFFF:FF,FFFFF,:,FFF::F:FFF,:FF:FFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#2:87-150  1:N:0:CGCAAGAA+ACTCCTAA
CCCTGTTCTTCTCCTGGGGAGGGAAGGGCAGGTCCCCTAATGCCTGAACCACATGGTGGGCCA
+
FFF::F:FFF,:FF:FFFFFFFFF,FFFFF:FFF,F,,FFF,F:FFFFFFF:F,,,FFF,F,,
</code></pre>

CC BY-SA 4.0
2020-07-14T04:54:58.570Z
2020-07-14T04:54:58.570Z
gringer
<p>I have bacterial chromosomes/plasmids that I have downloaded off refseq in a multi-fasta that I would like to fragment into non-overlapping fragments of between 1 kb and 15 kb.</p>
<p>I have tried googling for related tools but to no avail. Does anyone know how I would achieve this? Something very similar is done in this paper:
<a href="https://link.springer.com/article/10.1007/s40484-019-0187-4" rel="nofollow noreferrer">https://link.springer.com/article/10.1007/s40484-019-0187-4</a></p>
<p>Thanks in advance!</p>

<assembly><bacteria>
How to fragment genomes into non-overlapping sequences of differing sizes?
# How to fragment genomes into non-overlapping sequences of differing sizes?

(Asked by: username_0 on 2020-07-13)

I have bacterial chromosomes/plasmids that I have downloaded off refseq in a multi-fasta that I would like to fragment into non-overlapping fragments of between 1 kb and 15 kb.

I have tried googling for related tools but to no avail. Does anyone know how I would achieve this? Something very similar is done in this paper:
<https://link.springer.com/article/10.1007/s40484-019-0187-4>

Thanks in advance!

**Tags:** assembly, bacteria

---

## 2 answers

### Answer by username_1 on 2020-07-13. Score: 3 | Accepted answer

You could use numpy.random.randint to sample integers from a uniform distribution, from 1000 to 15000 nt, extracting subsequences of that length from a linearized multi-FASTA file.
Here's some code:
#!/usr/bin/env python

import sys
import numpy as np

low = int(sys.argv[1])
high = int(sys.argv[2])

def process_record(header, sequence, low, high):
    next = 0
    fragment = 1
    while next < len(sequence):
        new_next = np.random.randint(low, high) + next  
        if new_next > len(sequence) - high:
            new_next = len(sequence)
        sys.stdout.write('{}:{}\n{}\n'.format(header, fragment, sequence[next:new_next]))
        next = new_next
        fragment += 1

header = ""
sequence = ""
for line in sys.stdin:
    if line.startswith('>'):
        if len(sequence) > 0:
            process_record(header, sequence, low, high)
            sequence = ""
        header = line.rstrip()
    else:
        sequence += line.rstrip()

if len(sequence) > 0:
    process_record(header, sequence, low, high)

To use it:
$ ./13879.py 1000 15000 < in.fa > out.fa

Fragments will be disjoint (non-overlapping, contiguous) subsequences.
I add a check to determine if the last fragment of a FASTA record would be shorter than low (1kb) in length, and, if so, append that fragment to the final subsequence. It is possible for that last fragment to be up to 16 kb - 1 nt in length, given your parameters:
if new_next > len(sequence) - high:
    new_next = len(sequence)

You could change this logic to discard the last fragment (i.e. break), or simply allow fragments less than 1kb in size, depending on what you need.
EDIT
Here's a slight modification that should allow to you preserve fragment sizes within your bounds.
A cartoon may help illustrate the problem. Let's set up a FASTA record, which has been broken up into n fragments:
   1       2           n-2      n-1      n
|-----|---------|...|-------|---------|-----|

We may get unlucky and sample a random integer for fragment n-1 (frag(n-1)) which makes the length of frag(n) less than low (1 kb in your example).
In my first suggested solution, frag(n-1) and frag(n) are appended to one another. This is why you get some fragments longer than high (15 kb).
Another option is to discard frag(n), which is of length less than low (1 kb). Here's some code:
#!/usr/bin/env python

import sys
import numpy as np
import timeit

low = int(sys.argv[1])
high = int(sys.argv[2])
in_fn = sys.argv[3]

"""
Globals
"""
fragment_lengths = []

"""
Process a record
"""
def process_record(header, sequence, low, high):
    next = 0
    last = 0
    fragment = 1
    fragment_sequence = ""
    while last < len(sequence) - high:
        next = np.random.randint(low, high + 1) + last
        fragment_sequence += sequence[last:next]
        # Comment the sys.stdout.write and uncomment the fragment_lengths line, if measuring fragment length distributions 
        sys.stdout.write('{}:{}-{}:{}\n{}\n'.format(header, last, next, fragment, sequence[last:next]))
        #fragment_lengths.append(next - last)                                                                                                                                                                                                                                                                                                                                                                                          
        last = next
        fragment += 1
    """
    Uncomment to compare input sequence with fragments
    """
    #sys.stderr.write('In > {}\nOut > {}\n'.format(sequence, fragment_sequence))
    
"""
Parse records
"""
def parse_records(in_fn, low, high):
    header = ""
    sequence = ""
    with open(in_fn, 'r') as ifh:
        for line in ifh:
            if line.startswith('>'):
                if len(sequence) > 0:
                    process_record(header, sequence, low, high)
                    sequence = ""
                header = line.rstrip()
            else:
                sequence += line.rstrip()
        if len(sequence) > 0:
            process_record(header, sequence, low, high)

def count_elements(seq):
    hist = {}
    for i in seq:
        hist[i] = hist.get(i, 0) + 1
    return hist

def wrapper(func, *args, **kwargs):
    def wrapped():
        return func(*args, **kwargs)
    return wrapped

def main(in_fn, low, high):
    """
    Comment below if testing the distribution of fragments
    """
    parse_records(in_fn, low, high)
    """
    Uncomment to test the distribution of fragments from 
    running fns on same input 100k times
    """
    #wrapped = wrapper(parse_records, in_fn, low, high)
    #timeit.timeit(wrapped, number=100000)
    #print(count_elements(fragment_lengths))

if __name__ == "__main__":
    main(in_fn, low, high)

If frag(n) is needed, and you need to preserve the uniform distribution (or whatever fragment length distribution you want to model), then rejection sampling might be an option.
With rejection sampling, you keep fragmenting your input sequence, until one of the distributions of lengths gets you a set of fragments that entirely cover the input sequence.
If you need that, follow up and I'll take another look.

> **Comments:**
> **username_0**: Thanks for your answer! Do I need to remove all the fasta headers? I tried it on a multifasta and ended up with a few sequences over 15 kb. The distribution looks like this: https://imgur.com/a/FFr8xNv  
> **username_1**: No, you do not need to remove headers. You could end up with a few sequences over 15kb, if the last fragment cannot be sized between `low` and `high` (i.e. between 1 and 15kb). I'll propose a solution in the edit.  
> **username_0**: Thanks that worked perfectly!  
> **username_2**: Hi @StanleyHo please do sign the answer as "asscepted". I'll be making use of this code too, I like the args and kwargs  

---

### Answer by username_3 on 2020-07-14. Score: 1

I've previously written my own script to fragment DNA sequences (input as either fastq or fasta files) into identical-length fragments, with a given overlap. I used this because there was a particular assembler that demanded equal-length sequences, and I was working with nanopore reads, but it's since been useful for other purposes as well.
Because it is most concerned about making lengths identical, there may be some forced overlap when setting the overlap to zero. The overlap will only be exactly zero when the sequence length is an exact multiple of the fragment size. Here's an example use:
zero overlap
normalise_seqlengths.pl -fraglength 25 -overlap 0 oneRead.fq
@A00187:180:HLGNGDSXX:3:1101:8486:1031#0 [150 bp, 6 fragments, 0 overlap] 1:N:0:CGCAAGAA+ACTCCTAA
ATTGCGCAATGACAGCCAGGGGAGT
+
FFFFFFFFFFFFFFFFF:,FFFFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#1:25-50  1:N:0:CGCAAGAA+ACTCCTAA
GGGAGCTGTCTGGGGTAACCTGGCT
+
:FFFF:F:FFFFF::FFFFF,::FF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#2:50-75  1:N:0:CGCAAGAA+ACTCCTAA
CTAGGACCAGGTGGAAAAAGCACCT
+
F:FFFFFFFFFFFFFFFFFFFFFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#3:75-100  1:N:0:CGCAAGAA+ACTCCTAA
TGGCTGGCCCAGCCCTGTTCTTCTC
+
:FF,FFFFF,:,FFF::F:FFF,:F
@A00187:180:HLGNGDSXX:3:1101:8486:1031#4:100-125  1:N:0:CGCAAGAA+ACTCCTAA
CTGGGGAGGGAAGGGCAGGTCCCCT
+
F:FFFFFFFFF,FFFFF:FFF,F,,
@A00187:180:HLGNGDSXX:3:1101:8486:1031#5:125-150  1:N:0:CGCAAGAA+ACTCCTAA
AATGCCTGAACCACATGGTGGGCCA
+
FFF,F:FFFFFFF:F,,,FFF,F,,

partial overlap
normalise_seqlengths.pl -fraglength 63 -overlap 0 oneRead.fq
@A00187:180:HLGNGDSXX:3:1101:8486:1031#0 [150 bp, 3 fragments, 20 overlap] 1:N:0:CGCAAGAA+ACTCCTAA
ATTGCGCAATGACAGCCAGGGGAGTGGGAGCTGTCTGGGGTAACCTGGCTCTAGGACCAGGTG
+
FFFFFFFFFFFFFFFFF:,FFFFFF:FFFF:F:FFFFF::FFFFF,::FFF:FFFFFFFFFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#1:43-106  1:N:0:CGCAAGAA+ACTCCTAA
CCTGGCTCTAGGACCAGGTGGAAAAAGCACCTTGGCTGGCCCAGCCCTGTTCTTCTCCTGGGG
+
FF,::FFF:FFFFFFFFFFFFFFFFFFFFFFF:FF,FFFFF,:,FFF::F:FFF,:FF:FFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#2:87-150  1:N:0:CGCAAGAA+ACTCCTAA
CCCTGTTCTTCTCCTGGGGAGGGAAGGGCAGGTCCCCTAATGCCTGAACCACATGGTGGGCCA
+
FFF::F:FFF,:FF:FFFFFFFFF,FFFFF:FFF,F,,FFF,F:FFFFFFF:F,,,FFF,F,,


---

