sort a fasta file containing the Oxford Nanopore Technologies (ONT) header by sequencing start_time ascending

Question

I have basecalled ONT reads and converted them to multifasta. The multifasta contains the original ONT headers in this format:

3ebb56cd-3671-4abd-b1ac-0c759bt068d0 runid=eb6214851489c8e00eb0dcbd00d737f7dddxxxx read=1078 ch=422 start_time=2017-11-21T20:35:34Z

By checking the sequencing output as fasta, the sequences are not ordered by start_time.

>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z
TTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATG
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
TGGTATGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATG 
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
CTCGTTTCAGGTATGCTTCGTTTCAGTTGCCCTCAGGAGCAGATTGTTTATGGGTTACTTTCCTGTCGCTCTATCTTCAC 

I learned how to extract header IDs and time_start from ONT fastq files by using this answer.

I'd like to have a similar approach but with fasta.

Yet, I would like to create a new fasta file where the sequences are written in an ascending order as they were sequenced (start_time) as follows:

>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
TGGTATGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATG 
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
CTCGTTTCAGGTATGCTTCGTTTCAGTTGCCCTCAGGAGCAGATTGTTTATGGGTTACTTTCCTGTCGCTCTATCTTCAC
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z
TTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATG

Answer 1

Edit: After doing some reading, if it can be assumed that all the timestamps are ISO8601 and of the same representation of ISO861 (there are a few), then they can be sorted lexicographically, directly, without converting to epoch-time integers. So I have edited this answer to add an answer without conversion steps:

---

Here's one approach that uses the timestamps directly:

(1) Extract the start times from the headers and sort them to another file, in ascending order:

$ awk '(/^>/){ match($0, /start_time=(.*)/, a); print a[1]; }' my.fa | sort > starts_iso8601.sorted.txt

(2) Use that file of ordered ISO8601-formatted starts to write out FASTA records in your sort order:

$ while IFS='' read -r line || [[ -n "$line" ]]; do awk -v line=$line -v RS=">" '{ m = match($0,line); if (m) { printf(">%s", $0); } }' my.fa ; done < starts_iso8601.sorted.txt > my_ordered.fa

Quite a bit cleaner and faster.

You might separate the timestamp extraction and sort steps into two commands and inspect the intermediate result, just to ensure that you have valid strings. Or use the approach below to convert the timestamps, which will catch bad input.

---

Here's a second way that may add unnecessary work. As an advantage, however, it adds some error checking by converting the ISO8601-formatted timestamp.

If the string you pull out of the header is not a valid timestamp, the first conversion step (ISO8601-to-epoch) will fail and you will then be able to catch invalid input pretty quickly, which can add confidence that your output is correct.

(1) Extract the start times from the headers:

$ awk '(/^>/){ match($0, /start_time=(.*)/, a); print a[1]; }' my.fa > starts_iso8601.txt

(2) Convert the ISO8601-formatted times to time-since-epoch times:

$ while IFS='' read -r line || [[ -n "$line" ]]; do date -d"$line" +%s; done < starts_iso8601.txt > starts_epoch.txt

(3) Sort the epoch times — which are just integers — to another file, in ascending order:

$ sort -n starts_epoch.txt > starts_epoch.sorted.txt

(4) Convert sorted epoch starts back to ISO8601 format:

$ while IFS='' read -r line || [[ -n "$line" ]]; do date -d"@$line" -u +"%Y-%m-%dT%H:%M:%SZ"; done < starts_epoch.sorted.txt > starts_iso8601.sorted.txt

(5) Use that file of ordered starts to write out FASTA records in your sort order:

$ while IFS='' read -r line || [[ -n "$line" ]]; do awk -v line=$line -v RS=">" '{ m = match($0,line); if (m) { printf(">%s", $0); } }' my.fa ; done < starts_iso8601.sorted.txt > my_ordered.fa

Note: Use gawk, gdate, gsort in place of awk, date and sort if you're using OS X, Homebrew, and GNU core utilities.

I'm showing the steps separately to explain the general procedure.

Many or all of these steps can be glued together in a Unix pipeline to save the time, memory and disk space costs of writing intermediate files.

Also, I'm not doing much error checking here. You might look at cases where headers do not contain start_time, or times are not all in ISO8601, or other cases not anticipated here.

Answer 2

I've just made a small tweak to my fastx-sort.pl script to get this working. It works on both fasta and fastq sequences, and should even work on a random mixture of the two.

This script can now do a dictionary sort by the bits of an ID string that match an arbitrary search pattern, as well as the previous ID and length-based sorting:

fastx-sort.pl -p '<perl regex pattern>' <fastq file> > sorted.fastq

Here's an example with your fasta input, sorting on the nanopore time string:

$ cat seq.fa | grep start_time
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
$ ~/scripts/fastx-sort.pl -p 'start_time=(.*?Z)' seq.fa | grep start_time
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z

Here's a worked example with one of my fastq files for the nanopore time string:

$ cat called_all.fastq | grep start_time | head -n 3
@3a151247-6720-4071-945e-7a92c1856e45 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=9896 ch=246 start_time=2017-11-27T00:58:42Z
@52e439d6-f14a-4af2-adba-132d240fa967 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=8634 ch=511 start_time=2017-11-27T00:58:41Z
@82e3dd50-c395-4ade-ab6b-0eed530d3c3f runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=8312 ch=299 start_time=2017-11-27T00:58:39Z
$ ~/scripts/fastx-sort.pl -p 'start_time=(.*?Z)' called_all.fastq | grep start_time | head -n 10
@195ad363-1d5a-40bf-abed-820c8acce8fc runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=9 ch=144 start_time=2017-11-26T22:42:32Z
@5df2fc6e-53c8-489c-8002-c97f151104bd runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=11 ch=129 start_time=2017-11-26T22:42:37Z
@23940fe6-0eef-418f-9156-d0a26df9bd80 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=22 ch=330 start_time=2017-11-26T22:42:39Z

This script now does a stable sort, so you can do things like sorting by channel, then by read number:

~/scripts/fastx-sort.pl -n -p 'read=([0-9]+)' called_all.fastq | ~/scripts/fastx-sort.pl -n -p 'ch=([0-9]+)' | grep 'ch='
...
@876327ec-ccb0-44ad-a9fb-e008744d6c15 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=5974 ch=66 start_time=2017-11-27T00:04:52Z
@e3ec5d92-ee8a-43d6-82ed-e359e0539b75 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=6117 ch=66 start_time=2017-11-27T00:06:55Z
@410fe169-c8d1-4686-b92f-cbeb0602b2f0 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=6280 ch=66 start_time=2017-11-27T00:09:26Z
@d27c2827-1746-4827-9fec-3075e08d3d64 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=6375 ch=66 start_time=2017-11-27T00:10:44Z
@43079baa-8591-4498-bbfd-5f5ad677c011 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=32 ch=69 start_time=2017-11-26T22:43:32Z
@2bd9b317-438d-4e45-958b-7f77df1dd51a runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=170 ch=69 start_time=2017-11-26T22:47:32Z
@e3599944-849a-492c-ad6d-8cc319b7c5d7 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=207 ch=69 start_time=2017-11-26T22:49:11Z
@8b9d80d0-b28d-46bb-a4e9-6bba95a79ff4 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=320 ch=71 start_time=2017-11-26T22:50:33Z
@b964f164-c36c-4f8a-8a1a-1ad196c996a3 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=780 ch=71 start_time=2017-11-26T22:59:14Z
@92497c10-a46a-4a31-bc18-cfa334491ca9 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=937 ch=71 start_time=2017-11-26T23:02:32Z
@2a715984-6e14-4f31-bd47-d36ee5cbc3e8 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=2271 ch=71 start_time=2017-11-26T23:24:11Z
@ff2a59ae-4909-4f3e-bea2-b695d521bf5e runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=2441 ch=71 start_time=2017-11-26T23:27:01Z
@f106792d-31a4-44f2-a4fe-83f35a42a573 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=2559 ch=71 start_time=2017-11-26T23:28:54Z
@8de1d1ca-c6f7-43f5-bbd6-3c79187c1bbe runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=2930 ch=71 start_time=2017-11-26T23:34:33Z
...

Additional notes:

• In order to do the sort, the script currently stores the entirety of the input in memory.

Answer 3

My go-to tools for any type of fasta file manipulation are these two simple awk scripts:

• FastaToTbl will convert fasta to "tbl" formato: the ID, a tab and then the sequence, all on the same line.

  #!/usr/bin/awk -f
  {
          if (substr($1,1,1)==">")
      		if (NR>1)
                      	printf "\n%s\t", substr($0,2,length($0)-1)
      		else 
      			printf "%s\t", substr($0,2,length($0)-1)
              else 
                      printf "%s", $0
  }END{printf "\n"}
  

• TblToFasta will convert tbl back to fasta:

  #!/usr/bin/awk -f
  {
    sequence=$NF
  
    ls = length(sequence)
    is = 1
    fld  = 1
    while (fld < NF)
    {
       if (fld == 1){printf ">"}
       printf "%s " , $fld
  
       if (fld == NF-1)
        {
          printf "\n"
        }
        fld = fld+1
    }
    while (is <= ls)
    {
      printf "%s\n", substr(sequence,is,60)
      is=is+60
    }
  }
  

Armed with those and a bit of Perl, you can do:

$ FastaToTbl file.fa | 
    perl -lne '/start_time=(\S+)/; 
               $d=`date -d $1 +%s`; 
               $k{$d}=$_; 
               END{
                    print "$k{$_}" for sort(keys(%k))
               }' | TblToFasta
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z 
TGGTATGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTAT
GTGGTATGTGGTATGTGGTATGTGGTATG
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z 
CTCGTTTCAGGTATGCTTCGTTTCAGTTGCCCTCAGGAGCAGATTGTTTATGGGTTACTT
TCCTGTCGCTCTATCTTCAC
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z 
TTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGG
TATGTTGGTATGTTGGTATGTTGGTATG

Explanation

The perl script will take the tbl format sequences (so everything on one line), and identify the start_time: /start_time=(\S+)/. The parentheses capture the matched pattern so the start time will now be $1. Then, we run the system date command giving it the start time as the date and getting the result in seconds since the epoch (%s). note that this assumes GNU date (the default on Linux systems) which can handle the -d option to specify a date. The next step is using $d, the date in seconds, as a key for a hash whose values are the input lines: $k{$d}=$_. Then, after the whole file has been processed, we print each entry of the hash in sorted order (print "$k{$_}" for sort(keys(%k))). Finally, TblToFasta changes back to fasta format.