<p>For the specific regions you mentioned, you could do this in <code>awk</code>. It won't be fast for a large file, but it will work:</p>

<pre><code><span class="math-container">$ awk 'NR%4==2 || NR%4==0{$0=substr($0,1,1)substr($0,3,1)substr($</span>0,5,3)}1 ' file.fq 
@id1
ACGGT
+
13567
@id2
TCGGC
+
13567
</code></pre>

<p><code>$0</code> is the current line. The script above will replace the line with three substrings of itself, one starting at position <code>1</code> with a length of one, one starting at position 3 with a length of one and the last starting at plsition 5 with a length of 3. The <code>NR%4==2 || NR%4==0</code> ensures that this substitution only occurs on lines whose line number modulo 4 is either 2 or 0, so every 2nd and 4th line of the file.</p>

<p><strong>IMPORTANT</strong>: this assumes you <em>always</em> have only one line of sequence in each section. The fastq format allows multiple lines of sequence (and quality scores). So your data might not conform to this. However, for regular short-read data, it should be fine. </p>

<hr>

<p>And here's a slightly more sophisticated version in Perl, with some rudimentary error checking and which can take regions in the format you asked for (<code>1,3,5-7</code>):</p>

<pre class="lang-perl prettyprint-override"><code>#!/usr/bin/env perl

if (<span class="math-container">$#ARGV != 1) {
  die "Need exactly 2 arguments.\n";
}
my $</span>file = $ARGV[0];

my @regions = split(/,/,<span class="math-container">$ARGV[1]);
my (@startPositions, @lengths);
for my $</span>arg (0..<span class="math-container">$#regions) {
  if ($regions[$arg] =~ /^\d+$/) {
    push @startPositions, $regions[$arg];
  }
  elsif ($regions[$arg] =~ /^(\d+)-(\d+)$/) {
    for my $i ($1..$2) {
        push @startPositions, $</span>i;
    }
  }
}

open(my <span class="math-container">$inputFileHandle, '&lt;', $</span>file) or
    die "Failed to open <span class="math-container">$file for reading: $</span>!\n";

while (&lt;<span class="math-container">$inputFileHandle&gt;) {
  if ($. % 4 == 2 || $. % 4 == 0) {
    my $line = "";
    foreach my $start (@startPositions) {
      $line .= substr($_, $start-1, 1)
    }
    print "$line\n"
  }
  else {
    print
  }
  ## If we find a line that consists only of '+' but its line number
  ## modulo 4 isn't 3, that means that there was at least one multi-line
  ## sequence that we cannot handle correctly.
  if (/^\+$/ &amp;&amp; $. % 4 != 3) {
    die "Line $</span>. is just a '+', multi-line sequence found in file.\n"
  }
}
</code></pre>

<p>Save that as <code>foo.pl</code> and then run like this:</p>

<pre><code>perl foo.pl file.fq 1,3,5-7
</code></pre>

CC BY-SA 4.0
2020-05-11T11:28:13.973Z
2020-05-11T12:20:28.103Z
2020-05-11T12:20:28.103Z
terdon
<p>This looks like a job for samtools fqidx (<a href="http://www.htslib.org/doc/samtools-fqidx.html" rel="nofollow noreferrer">http://www.htslib.org/doc/samtools-fqidx.html</a>), assuming that your file is not enormous... read the warning in the link above.</p>

<p>usage would be</p>

<p><code>samtools fqidx yourfile.fq &lt;region&gt;</code></p>

<p>where region is specified using the fastq record identifier and coordinates.</p>

CC BY-SA 4.0
2020-05-11T10:54:40.007Z
How would a region like what the OP wants be written? You can run `samtools fqidx file.fq id1` to get the data for `id1`, but how would you specify only the 1st, 3rd, 5th, 6th and 7th residue as requested by the OP? The docs don't show the format. And wouldn't this need each id in the fastq file to be mentioned individually? How would you use this for what the OP asked for?
terdon
CC BY-SA 4.0
2020-05-11T12:26:27.433Z
I second terdons, fqidx is certainly a useful tool, but it doesn't address the specific question here.
Sebastian Müller
CC BY-SA 4.0
2020-05-11T08:21:19.313Z
2020-05-11T08:21:19.313Z
Wouter De Coster
<p>Given a FASTQ file, I'd like to generate a new FASTQ file including only certain subsections of the original sequences specified by their positions. </p>

<p>For example, say I want to extract the nucleotides at position <strong>1, 3, 5-7</strong> from each sequence:</p>

<p>Input:</p>

<pre><code>@id1
AACCGGTCC
+
123456789
@id2
TACCGGCCC
+
123456789
</code></pre>

<p>Output:</p>

<pre><code>@id1
ACGGT
+
13567
@id2
TCGGC
+
13567
</code></pre>

<p>Note, the phred-scores are constructed such that they correspond to the position for easier reading.
Note2, the phred-score subsections should be included accordingly (i.e. the same positions)</p>

<p>Ideally, a command should work this way:</p>

<pre><code>subselect 1,3,5-7 input.fq &gt; output.fq
</code></pre>

<p>Having looked through several poplar software, I couldn't find anything working out of the box. I was thinking to use either awk or even cutadapt (systematically trimming from the 5<span class="math-container">`</span> and pasting) but maybe out of an box solution does already exist which evaded me. </p>

<fastq><awk>
Extracting specific sections from FASTQ determined by positions
# Extracting specific sections from FASTQ determined by positions

(Asked by: username_0 on 2020-05-11)

Given a FASTQ file, I'd like to generate a new FASTQ file including only certain subsections of the original sequences specified by their positions. 
For example, say I want to extract the nucleotides at position 1, 3, 5-7 from each sequence:
Input:
@id1
AACCGGTCC
+
123456789
@id2
TACCGGCCC
+
123456789

Output:
@id1
ACGGT
+
13567
@id2
TCGGC
+
13567

Note, the phred-scores are constructed such that they correspond to the position for easier reading.
Note2, the phred-score subsections should be included accordingly (i.e. the same positions)
Ideally, a command should work this way:
subselect 1,3,5-7 input.fq > output.fq

Having looked through several poplar software, I couldn't find anything working out of the box. I was thinking to use either awk or even cutadapt (systematically trimming from the 5` and pasting) but maybe out of an box solution does already exist which evaded me.

**Tags:** fastq, awk

---

## 2 answers

### Answer by username_1 on 2020-05-11. Score: 1 | Accepted answer

For the specific regions you mentioned, you could do this in awk. It won't be fast for a large file, but it will work:
$ awk 'NR%4==2 || NR%4==0{$0=substr($0,1,1)substr($0,3,1)substr($0,5,3)}1 ' file.fq 
@id1
ACGGT
+
13567
@id2
TCGGC
+
13567

$0 is the current line. The script above will replace the line with three substrings of itself, one starting at position 1 with a length of one, one starting at position 3 with a length of one and the last starting at plsition 5 with a length of 3. The NR%4==2 || NR%4==0 ensures that this substitution only occurs on lines whose line number modulo 4 is either 2 or 0, so every 2nd and 4th line of the file.
IMPORTANT: this assumes you always have only one line of sequence in each section. The fastq format allows multiple lines of sequence (and quality scores). So your data might not conform to this. However, for regular short-read data, it should be fine. 

And here's a slightly more sophisticated version in Perl, with some rudimentary error checking and which can take regions in the format you asked for (1,3,5-7):
#!/usr/bin/env perl

if ($#ARGV != 1) {
  die "Need exactly 2 arguments.\n";
}
my $file = $ARGV[0];

my @regions = split(/,/,$ARGV[1]);
my (@startPositions, @lengths);
for my $arg (0..$#regions) {
  if ($regions[$arg] =~ /^\d+$/) {
    push @startPositions, $regions[$arg];
  }
  elsif ($regions[$arg] =~ /^(\d+)-(\d+)$/) {
    for my $i ($1..$2) {
        push @startPositions, $i;
    }
  }
}

open(my $inputFileHandle, '<', $file) or
    die "Failed to open $file for reading: $!\n";

while (<$inputFileHandle>) {
  if ($. % 4 == 2 || $. % 4 == 0) {
    my $line = "";
    foreach my $start (@startPositions) {
      $line .= substr($_, $start-1, 1)
    }
    print "$line\n"
  }
  else {
    print
  }
  ## If we find a line that consists only of '+' but its line number
  ## modulo 4 isn't 3, that means that there was at least one multi-line
  ## sequence that we cannot handle correctly.
  if (/^\+$/ && $. % 4 != 3) {
    die "Line $. is just a '+', multi-line sequence found in file.\n"
  }
}

Save that as foo.pl and then run like this:
perl foo.pl file.fq 1,3,5-7


---

### Answer by username_2 on 2020-05-11. Score: -1

This looks like a job for samtools fqidx (<http://www.htslib.org/doc/samtools-fqidx.html>), assuming that your file is not enormous... read the warning in the link above.

usage would be

`samtools fqidx yourfile.fq`

where region is specified using the fastq record identifier and coordinates.

> **Comments:**
> **username_1**: How would a region like what the OP wants be written? You can run `samtools fqidx file.fq id1` to get the data for `id1`, but how would you specify only the 1st, 3rd, 5th, 6th and 7th residue as requested by the OP? The docs don't show the format. And wouldn't this need each id in the fastq file to be mentioned individually? How would you use this for what the OP asked for?  
> **username_0**: I second terdons, fqidx is certainly a useful tool, but it doesn't address the specific question here.  

---

