13

fast5 is a variant of HDF5 the native format in which raw data from Oxford Nanopore MinION are provided. You can easily extract the reads in fast5 format into a standard fastq format, using for example poretools.

Say I have aligned these reads in fastq format to an external reference genome, resulting in a SAM file. Say I have then taken a subset of the SAM file, according to the bitwise flag, to include only the reads that map to the reference. With the read ID, I can then grep them out from the file containing the reads in fastq format, generating a subset file in fastq format containing only the IDs that have mapped to the reference.

Now my question is, can we subset reads from the fast5 archive according to the list of mapping reads as taken from the file with reads in fastq format? This is for educational purposes, so that we have a smaller starting archive, and the fast5 -> fastq extraction takes less cpu time.

Michael Hall
  • 663
  • 4
  • 11
aechchiki
  • 2,676
  • 11
  • 34
  • 2
    This seems a nice question to play around if we could have sample data. Could you provide (or a link to) a fast5 and the list of mappings ? – llrs Aug 15 '17 at 16:18
  • 1
    There are some sample reads from my TEDx talk last year here, reads from runs I did during the AGTA conference last year here, and reads from a chimeric read investigation here. The ONT FAST5 "standard" change so frequently that it's difficult to find something that represents what's currently available. – gringer Aug 15 '17 at 18:09

2 Answers2

6

NOTICE: I have altered my answer slightly from the original as I have turned the original script into a pip installable program (with tests) and have updated the links and code snippets accordingly. The essence of the answer is still exactly the same.

This is something I have been meaning to get around to for a while, so thanks for the prompt.

I have created a python program called fast5seek to do what (I think) you're after.

As you have mentioned this is for educational purposes I have added a tonne of comments to the code too so I think you shouldn't have any issues following it.

The docs on the GitHub repo have all the info, but for those reading along at home

pip3 install fast5seek
fast5seek -i /path/to/fast5s -r in.fastq in.bam in.sam -o out.txt

What it does is read in <in.fastq|in.bam|in.sam> and extract the read id from each header. It then goes through all the fast5 files under /path/to/fast5s and checks whether their read id is in the set of read ids from <in.fastq|in.bam|in.sam>. If it is, the path to the file is written to it's own line in out.txt.

If no output (-o) is given, it will write the output to stdout.

So if you wanted to pipe these paths into another program, you could do something like

mkdir subset_dir/
fast5seek -i /path/to/fast5s/ -r in.fastq | xargs cp -t subset_dir/

The above example would copy the fast5 files that are found in your fastq/BAM/SAM to subset_dir/.

Michael Hall
  • 663
  • 4
  • 11
  • 1
    As a caution, the run id should be checked (if that information is available). There's a small possibility that two FASTQ files from different runs will have the same read id. – gringer Aug 16 '17 at 11:40
  • 1
    I have now added this check for fastq files. Although, run_id is not contained within BAM/SAM files (at least not for the ones I have) so this check will only occur on fastq. – Michael Hall Aug 17 '17 at 00:56
3

I do this very frequently, using the read file name to identify FAST5 files associated with particular reads. If a FASTQ record includes the channel number and read number (and preferably the runID as well), then I use that information to find the associated FAST5 read.

If reads are called with Albacore, then the sequencing_summary.txt file has additional information. The first column is the exact name of the file that was called, the second column is the read ID of the FASTQ sequence, and the third column is the run ID associated with that read ID. This is much easier to work with, but does require calling reads with Albacore (which, admittedly, is what seems to be producing the best results at the moment).

gringer
  • 14,012
  • 5
  • 23
  • 79
  • yes this is what I did, basically: 1. extract fast5 to fastq; 2. align the fastq to the reference; 3. get the reads ID that map from the SAM of alignment; 4. grep the mapping reads ID from the fastq; 5. extract the name of the fast5 from the fasta header of each read (last field); 6. find the corresponding fast5 in the original fast5 directory. But I was wondering if I did not miss any software that could do this on one run. – aechchiki Aug 16 '17 at 05:38
  • 1
    Added a note about albacore, which has a direct link between FASTQ read IDs and FAST5 file names – gringer Aug 16 '17 at 07:39
  • 1
    In light of your above comment @AminaEchchiki I have added in compatibility for BAM and SAM files too. – Michael Hall Aug 16 '17 at 08:25