Extracting specific sections from FASTQ determined by positions

Question

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.

Answer 1

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 2

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 <region>

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