Subsampling using 'head -n #"?

[kga1978]

Hi all,

I often subsample my fastq files by using the unix 'head' command, rather than randomly getting reads from random positions in my fastq file. My setup is as follows:

Casava output

file1.fastq.gz
file2.fastq.gz
.
.
file#.fastq.gz

I concatenate these files using the following:

gunzip *.gz -c | gzip -9 > file.fastq.gz

In this way, the first part of this file will contain the reads from the original 'file1.fastq.gz'. I then subsample 25.000 reads from this file using the following command:

head -n 100000 file.fastq.gz | <downstream analysis, e.g. blastn>

My worry is that by doing that, I will get some sort of bias in my analysis as I am only taking the 'head' of the first part of all my reads. The question is - is this a valid concern? I.e. are the reads in the first part of, say, 'file1.fastq.gz' somehow different than, say, the middle part of 'file4.fastq.gz'?

Thanks very much in advance

[jbrwn]

it's definitely not very rigorous. htseq has a fastq reader in it which could help you subsample randomly.
htseq tour
previous discussion on subsampling paired-end data:
http://seqanswers.com/forums/showthread.php?t=12070

[gringer]

I do this often (mostly to make sure that the output data format is correct) and often find that it's a really bad representation of the data due to strange biases, low mapping quality, and odd mappings with unequal distribution to name a few. SOLiD or Illumina output files have a pile of rubbish at the start of the file from bad sequencing reads (the edges of chips / cells seem to be more prone to error), which heavily influences the results gleaned from the first few reads.

If you've been already doing this, I'm somewhat surprised you haven't discovered this already.

If you can afford the time and want something more representative but keeping within quick command-line stuff, use 'shuf -n' instead of 'head -n'. You could also try dumping the first million or so reads (or more), then take the next million or so (but there'll still be some bias there as well).

[kga1978]

Thanks for the replies - that's very helpful.

I have only been running this over the last couple of days (but we have loads where I need to do it) and didn't notice anything majorly wrong, other than some weird QC scores for a couple of files (gringer, I guess that is what you have seen as well).

We have a script to subsample correctly, so I will probably just use that - the main problem being that it needs to read the entire file before being able to output a subsample. In that sense 'head' was a lot faster.

gringer, as for shuf, I don't think that will work? Reason being that each read spans 2 (fasta) or 4 (fastq) lines? Otherwise, this is exactly what I would be looking for - given a file of n lines, randomly take a read from line k (so 4 lines, always starting at a position that can be divided by 4). But I haven't been able to find a script or command that will do that.

[maubp]

Originally Posted by kga1978

I then subsample 25.000 reads from this file using the following command:

head -n 100000 file.fastq.gz | <downstream analysis, e.g. blastn>

I assume that is a typo, since you must unzip the gripped FASTQ before taking the first 1000000 lines.

Originally Posted by kga1978

IMy worry is that by doing that, I will get some sort of bias in my analysis as I am only taking the 'head' of the first part of all my reads. The question is - is this a valid concern?

Yes, since the read order produced by Illumina reflects the slide positions, so the first and last reads will be from the edges and generally of poorer quality (so I've read).

[gringer]

Originally Posted by kga1978

gringer, as for shuf, I don't think that will work? Reason being that each read spans 2 (fasta) or 4 (fastq) lines?

You're correct... it won't work because of the serial nature of the files. Yet another case of me not thinking enough before posting.

The only ways I can think of at the moment are fairly complex piped shell commands, or small programs.

I guess one "simple" way to do it would be to seek to a random file position, then read forward until you find a start of sequence marker, but that's a lot of work for too little benefit. FWIW, doing a strict random seek on a fastq file is not possible because both @ and + can appear in the quality strings, and you can get multi-line fastq files.

Edit: I've just thought of another way that would work for a single pass (i.e. could work for piped input), assuming you knew the approximate number of records in the file. You decide in advance what records will be selected (i.e. generate 1M numbers between 1 and 100M), sort the list, then run through the file and print only the records that appear in the list.

[kga1978]

I assume that is a typo

Sorry, yes, that was a typo - there's a gunzip in front of that.

So, a talented student wrote up a quick script in python for me and it works really well. It appears that it is about as fast as writing to disk, so overall what I was looking for. It takes the following argument:

python subsampler.py <filename> <number of sequences to extract>

The script is as follows:

#!/usr/bin/env python
import random
import sys

#name of the input file (fasta or fastq)
#assumes input file is standard fasta/fastq format
fileName = sys.argv[1]
#number of sequences to subsample
numSeq = int(sys.argv[2])
increment = 0

#if it's a fasta file
if (fileName.find(".fasta") != -1):
  increment = 2
#else if it's a fastq file
elif (fileName.find(".fastq") != -1):
  increment = 4
#quit if neither
else:
  sys.stdout.write("not a fasta/fastq file\n")
  sys.exit()

FILE = open(fileName, 'r')
reads = list()
index = 0
#put the entire file in a list
for line in FILE:
  if(index % increment == 0):
    reads.append(line)
  else:
    reads[(index/increment)] += line
  index += 1
FILE.close()

#figure out total number of reads, error if applicable
totalReads = index/increment
if(totalReads < numSeq):
  sys.stdout.write("You only have "+str(totalReads)+" reads!\n")
  sys.exit()

#shuffle the reads!
random.shuffle(reads)
#output the reads to stdout
for i in range(0, numSeq):
  sys.stdout.write(reads[i])

[gringer]

Originally Posted by kga1978

So, a talented student wrote up a quick script in python for me and it works really well. It appears that it is about as fast as writing to disk, so overall what I was looking for.

This code is fine, but bear in mind that while this process will work on piped input (possibly with minor modifications to file opening), it requires your entire input file to be stored in memory (in the 'reads' list). If you don't care about piped input and can afford a second sweep through the file, you can reduce the memory footprint substantially by only using the first sweep to count records, then the second sweep to print out the randomly selected sequences.

Edit: Just noticed that it won't work with multi-line FASTA or FASTQ files, but for most NGS things you don't need to worry about that.

[ETHANol]

kga1978, thanks for posting the script. I was just working on a perl script to do the same, but being a beginner at scripting it was going to take me an annoying amount of time. So thanks and now I can do something more productive today!!!

[kga1978]

Hey gringer,

Thanks for your comments - very helpful and I'll make sure to pass them on. For now though, I don't really mind holding everything in memory - should make it faster and memory isn't a problem (for now....).

And actually, I won't work on stdin at the moment as it requires the .fasta or .fastq extension to work (this could be easily fixed by specifying file-type though...).

[kga1978]

Just to follow up on the first part of my question, I now ran a couple of subsample experiments and compared them to the entire dataset as well as to 'head -n' and 'tail -n'. As noted by other posters, there are definitely biases when doing 'head' and 'tail'.

I subsampled 1 million reads from a total of 50 million and ran FastQC (to calculate avg. Q-score) and aligned the reads to hg19 using BWA. This is what I got:

[kga1978]

Also, following gringer's suggestion, here's a less memory-hungry version of the downsampling script:

#!/usr/bin/env python
import random
import sys

random.seed()
#name of the input file (fasta or fastq)
#assumes input file is standard fasta/fastq format
fileName = sys.argv[1]
#number of sequences to subsample
numSeq = int(sys.argv[2])
increment = 0

#if it's a fasta file
if (fileName.find(".fasta") != -1):
  increment = 2
#else if it's a fastq file
elif (fileName.find(".fastq") != -1):
  increment = 4
#quit if neither
else:
  sys.stdout.write("not a fasta/fastq file\n")
  sys.exit()

FILE = open(fileName, 'r')
totalReads = list()
index = 0
for line in FILE:
  if(index % increment == 0):
    totalReads.append(index/increment)
  index += 1
FILE.close()
if(len(totalReads) < numSeq):
  sys.stdout.write("You only have "+str(len(totalReads))+" reads!\n")
  sys.exit()

ttl = len(totalReads)
random.shuffle(totalReads)
totalReads = totalReads[0: numSeq]
totalReads.sort()

FILE = open(fileName, 'r')
listIndex = 0

for i in range(0, ttl):
  curRead = ""
  for j in range(0, increment):
    curRead += FILE.readline()
  if (i == totalReads[listIndex]):
    sys.stdout.write(curRead)
    listIndex += 1
    if(listIndex == numSeq):
      break
FILE.close()

[Adjuvant]

Nobody asked for it, but I wanted it so maybe somebody else can use it, too.

Here's a rewrite of the python script that can now take two paired files as input and maintain pairing in the subsets. It can also still take just a single-end file as input. It's not pretty, but it does the job.

Changes: Now have to enter the number of reads desired first, then the input file name(s). Instead of going to stdout, the reads are output to out_1.fastq (or .fasta) and out_2.fastq (or .fasta). It would be relatively simple to put in an output prefix parameter, but this was good enough for my purposes.

#!/usr/bin/env python
import random
import sys

random.seed()
#number of sequences to subsample
numSeq = int(sys.argv[1])
#name of the input file (fasta or fastq)
#assumes input file is standard fasta/fastq format
fileName1 = sys.argv[2]
fileName2 = sys.argv[3]
increment = 0

#if it's a fasta file
if (fileName1.find(".fasta") != -1):
  increment = 2
#else if it's a fastq file
elif (fileName1.find(".fastq") != -1):
  increment = 4
#quit if neither
else:
  sys.stdout.write("not a fasta/fastq file\n")
  sys.exit()

FILE1 = open(fileName1, 'r')
totalReads1 = list()
index = 0
for line in FILE1:
  if(index % increment == 0):
    totalReads1.append(index/increment)
  index += 1
FILE1.close()
if(len(totalReads1) < numSeq):
  sys.stdout.write("You only have "+str(len(totalReads))+" reads!\n")
  sys.exit()
  
if (fileName2):
  FILE2 = open(fileName2, 'r')
  totalReads2 = list()
  index = 0
  for line in FILE2:
    if (index % increment == 0):
      totalReads2.append(index/increment)
    index += 1
  FILE2.close()
  if (len(totalReads1) != len(totalReads2)):
    sys.stdout.write("read counts do not match\n")
    sys.exit()

ttl = len(totalReads1)
random.shuffle(totalReads1)
totalReads1 = totalReads1[0: numSeq]
totalReads1.sort()

FILE1 = open(fileName1, 'r')
if (fileName2):
  FILE2 = open(fileName2, 'r')
listIndex = 0

if (increment == 4):
  OUTFILE1 = open('out_1.fastq', 'w')
  if (fileName2):
    OUTFILE2 = open('out_2.fastq', 'w')
else:
  OUTFILE1 = open('out_1.fasta', 'w')
  if (fileName2):
    OUTFILE2 = open('out_2.fasta', 'w')

for i in range(0, ttl):
  curRead1 = ""
  curRead2 = ""
  for j in range(0, increment):
    curRead1 += FILE1.readline()
    if (fileName2):
      curRead2 += FILE2.readline()
  if (i == totalReads1[listIndex]):
    OUTFILE1.write(curRead1)
    if (fileName2):
      OUTFILE2.write(curRead2)
    listIndex += 1
    if(listIndex == numSeq):
      break
FILE1.close()
if (fileName2):
  FILE2.close()
OUTFILE1.close()
if (fileName2):
  OUTFILE2.close()

[silver_steve]

Originally Posted by kga1978

Sorry, yes, that was a typo - there's a gunzip in front of that.

So, a talented student wrote up a quick script in python for me and it works really well. It appears that it is about as fast as writing to disk, so overall what I was looking for. It takes the following argument:

python subsampler.py <filename> <number of sequences to extract>

The script is as follows:

#!/usr/bin/env python
import random
import sys

#name of the input file (fasta or fastq)
#assumes input file is standard fasta/fastq format
fileName = sys.argv[1]
#number of sequences to subsample
numSeq = int(sys.argv[2])
increment = 0

#if it's a fasta file
if (fileName.find(".fasta") != -1):
  increment = 2
#else if it's a fastq file
elif (fileName.find(".fastq") != -1):
  increment = 4
#quit if neither
else:
  sys.stdout.write("not a fasta/fastq file\n")
  sys.exit()

FILE = open(fileName, 'r')
reads = list()
index = 0
#put the entire file in a list
for line in FILE:
  if(index % increment == 0):
    reads.append(line)
  else:
    reads[(index/increment)] += line
  index += 1
FILE.close()

#figure out total number of reads, error if applicable
totalReads = index/increment
if(totalReads < numSeq):
  sys.stdout.write("You only have "+str(totalReads)+" reads!\n")
  sys.exit()

#shuffle the reads!
random.shuffle(reads)
#output the reads to stdout
for i in range(0, numSeq):
  sys.stdout.write(reads[i])

Thank you kga1978. I ran this script successfully, yet the output was 1000 lines, not 1000 random sequences. I believe this is because my file is formatted like:
>File_Name
GACGAACGCTGGCGGCGTGCCTAATACATGCAAGTCGAGCGAGTCTGCCTTGAAGATCGG
AGTGCTTGCACTCTGTGAAACAAGATACAGGCTAGCGGCGGACGGGTGAGTAACACGTGG...
i.e., the sequence is not one continuous string, but rather a set of lines.
Do you know how to modify the code or modify my files to work with subsampler.py?

[gringer]

Edit: Just noticed that it won't work with multi-line FASTA or FASTQ files, but for most NGS things you don't need to worry about that.

So we have a situation that the script wasn't designed for. If you have fastx-toolkit handy, you can convert multi-line fasta files into single-line fasta files using the fasta-formatter tool, and then run through the script. Changing this script to work in a memory-efficient manner for multi-line fasta and fastq files would increase the complexity a reasonable amount.

A rough pseudo-code for this is as follows:

First pass -- counting records

1. Look for magic record start sequences, if found, increment the record counter
2. If a fasta file, read through and discard up to the next record start, return to the first step
3. If a fastq file, read up to the next '+' to work out sequence length
4. read the next line and discard
5. read as many quality values as were in the sequence, return to the first step

Second pass -- printing sequences

1. Generate n random numbers between 1 and the number of records
2. iterate through the sequences as before, but print the entire sequence record if the current record number matches to the random list

[swbarnes2]

I think a simpler way to do it...just get all the reads from a particular set of tiles (I don't know how that translates on HiSeq machines) Since the first few and last few tiles are bad, pick some tiles in the middle, by grepping the read names, which have teh cluster cooridnates embedded in them. Or zgrep, if the fastqs are already zipped. That should be random with respect to where the read align to the genome, with better quality than picking the top and bottom tiles with head or foot.

But yes, this won't work on multi-line fastqs.

[jmun7616]

Hey, just notices this thread. This Python solution may be a bit late but it will work on both multi-line fasta and fastq files. Set up is such that all files in same dir will be subsampled, but this can be changed. Also three subsamplings taken, but you can change this easily enough or just use one of them.

import random
import glob
inFiles = glob.glob('*.fasta') + glob.glob('*.fastq')
outFiles = []
num = int(raw_input("Enter number of random sequences to select:\n"))

for i in range(len(inFiles)):

        for k in range(3):
                fNames = []
                fSeqs = []
                same = 0
                outFiles.append(file(str(inFiles[i])+'_'+'Rand_'+str(num)+'-'+str(k+1)+'.fa', 'wt'))
                for line in open(inFiles[i]):
                        if ((line[0] == '>') or (line[0] == '@')):
                                fNames.append(line)
                                same = 0
                        else:
                                if (same != 1):
                                        fSeqs.append(line)
                                        same = 1
                                else:
                                        fSeqs[(len(fSeqs)-1)] += line
                curr = (len(outFiles)-1)
                for j in range(num):
                        a = random.randint(0, (len(fNames)-1))
                        outFiles[curr].write(fNames.pop(a))
                        outFiles[curr].write(fSeqs.pop(a))
raw_input("Done.")

[gringer]

This Python solution may be a bit late but it will work on both multi-line fasta and fastq files.

A few warnings about this Python code:

• ASCII quality scores in fastq files can contain arbitrary characters like '>' and '@'. The only way to protect against this for multi-line fastq files is to keep a record of the sequence length.
• This code reads the entire set of reads into memory, then outputs random sequences -- this will not work on fasta/fastq files that are larger than the memory limit. While this is the only way to do it for streamed input, the code here uses file names, so can do a two-pass operation and reopen a file (see my previous post).

[Richard Finney]

A simple script extracting random fastq records using reservoir sampling ...

# usage: program inputfastq numsubsamples outputsubampfqfile
# example: ./job.fqreservoir SRR451888.fq 1000 subsamp.fq

cat $1 | \
awk '{if ((p%4)<3) printf("%s\t",$0); else printf("%s\n",$0); p++;}' | \
awk -v k=$2 'BEGIN{srand(systime() + PROCINFO["pid"]);}{y=x++<k?x-1:int(rand()*x);if(y<k)a[y]=$0}END{for(z in a)print a[z]}' | \
tr '\t' '\n' > $3

See previous discussion here: http://seqanswers.com/forums/showthread.php?t=16845

Note the acutal line #3 is a hack of lh3's script WITH seeding, i.e. BEGIN{srand(systime() + PROCINFO["pid"]);}
Remove the seed code to get same results every time you run it. Otherwise leave it in to get different results.

This program joins the 4 fastq lines together into one line , spearating them with a tab character.
These combined lines are subsampled.
At the end, the "tr" command seperates the selected combined lines back into 4 lines.
This script used the UNIX/LINUX/BSD standard utilities "tr", "cat", and "awk".

[gringer]

This program joins the 4 fastq lines together into one line , spearating them with a tab character.

That particular code won't work with multi-line fasta/fastq files (the linked subsampler.py script has the same issue, as well as pre-reading all the lines to determine the total number of lines, and some modulus code that may influence probabilities in a bad way), but the reservoir sampling method used would provide a way to process files from standard input in a single pass, and use multi-line files from standard input. However, you need to be a bit careful with the selection procedure to ensure a uniform chance of selection for every record. After a bit of hunting, I see some reasonable pseudo-code here (linked from a stackoverflow thread):

http://blogs.msdn.com/b/spt/archive/...-sampling.aspx

The modified pseudo-code for a fasta/fastq file using reservoir sampling would be as follows:

1. Look for magic record start sequences. If found, increment the record counter and clear 'current record', otherwise print out the output array and finish
2. If a fasta file, read through and store up to the next record start as the 'current record', then go to step 6
3. If a fastq file, read up to the next '+' to work out sequence length (store in 'current record')
4. read the next line and store in 'current record'
5. read as many quality values as were in the sequence, store in 'current record'
6. if record counter is less than desired number of records, add current record to end of output array, return to step 1
7. generate a random number between 0 and the total number of records read
8. if this number is less than the number of desired records, replace the array element indexed by that number by the current record, return to step 1

The selection procedure here is (I think) the same as in the awk script that Richard Finney wrote / copied, which is (I think) the same as in the msdn pseudo-code.