3
$\begingroup$

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
$\endgroup$
3
  • $\begingroup$ The answer you link to is about fastq, but you talk about fasta in your question and show a line that isn't either. You have received one answer which assumes your input is fasta and one which assumes your input is fastq. Please edit your question and include an example of your input file so we can know which is correct. Ideally, show us a few lines of the original file and then show us the output you would like to see from those lines. $\endgroup$
    – terdon
    Dec 6, 2017 at 9:32
  • $\begingroup$ Thank you @terdon. I have edited the question to clarify! $\endgroup$
    – andresito
    Dec 6, 2017 at 19:41
  • $\begingroup$ Perfect! Thank you, now it's perfectly clear! By the way, this site supports a specific version of markdown syntax for formatting, you can see the details here: formatting tools. $\endgroup$
    – terdon
    Dec 6, 2017 at 22:58

3 Answers 3

3
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Thank you, @Alex Reynolds. Everything is working fine, however, the last code (5) seem to be writing some sequences twice. I can notice this because the resulting file my_ordered.fa contains more sequences than the original file my.fa . Could you check on that, please? $\endgroup$
    – andresito
    Dec 6, 2017 at 22:01
  • $\begingroup$ Can you provide an example FASTA file which is having this problem? $\endgroup$ Dec 6, 2017 at 23:33
2
$\begingroup$

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.
$\endgroup$
2
  • $\begingroup$ Thank you @gringer. Can you modify your script to use fasta and write an ascending ordered fasta as output? $\endgroup$
    – andresito
    Dec 6, 2017 at 22:27
  • 1
    $\begingroup$ It works with both fasta and fastq (and even a mixture of the two). The sort order is ascending by default. $\endgroup$
    – gringer
    Dec 7, 2017 at 0:05
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.