Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-14T16:26:51.910Z
Answers.Text: 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
Answers.UserDisplayName: Stanley Ho
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-14T20:10:38.237Z
Answers.Text: 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.
Answers.UserDisplayName: Alex Reynolds
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-15T10:47:27.460Z
Answers.Text: Thanks that worked perfectly!
Answers.UserDisplayName: Stanley Ho
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-15T10:57:53.693Z
Answers.Text: Hi @StanleyHo please do sign the answer as "asscepted". I'll be making use of this code too, I like the args and kwargs
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-13T19:33:37.747Z
Answers.LastActivityDate: 2020-07-15T18:58:53.347Z
Answers.LastEditDate: 2020-07-15T18:58:53.347Z
Answers.UserDisplayName: Alex Reynolds
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-14T04:54:58.570Z
Answers.LastActivityDate: 2020-07-14T04:54:58.570Z
Answers.UserDisplayName: gringer
Body: <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>
  
Tags: <assembly><bacteria>
Title: How to fragment genomes into non-overlapping sequences of differing sizes?
ThreadText: # 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](https://en.wikipedia.org/wiki/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 @username_0 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](https://gitlab.com/username_3/bioinfscripts/-/blob/master/normalise_seqlengths.pl) 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,,
  
  ```
  
  
  ---
  
  