0

I'm using the following chain of commands to sort a FASTQ file:

cat reads.fq | paste - - - - | sort -k1,1 -S 3G | tr '\t' '\n' > sorted_reads.fq

The FASTQ file is separated into groups for four line first of which is the id we want to sort on.

But in the output I see this pattern:

@ERR894725.10000000/1
AGAAAAAGAAAAATTGTTTCCTGTTATATCCATTTCCTTTAATATAGTTTACAAATTGTGCATTTCAACAGCAGCACTCTGTCCATATGTCTAGCAGACTCCTTAACACACTGTGCTACAACTTCT
+
<B<<B<<<B<<<<<<<B<B<BB<BB<<<<<BB<BB<BB<B<B0<<<<BBB<<BBB<B<B<BB<<B<BB<B<BB0BB<B<B<<<<<<BB<B<BB<BB<B<BBBBBBB0B<B<BBB<B0B<BB<<07<
@ERR894725.100000002/1
AGATAGAGTCTTGCTCTGTCACCCAGCCTGGAGTGCAACGGTGCTATCTCTACTAACTGCAACCTCCGCCCCCCAGGTTCAAGCGATTCTCCTGCCTCAGCCTCCTGAGCAGCTGGGACTACAGTG
+
<<<<<BB<<<B<<BB<B<B<B<BBBBBBB<BBBB<BBB<7BB<BB<<<B<B<<B<B<B<BBB<<B<B7BBBBBBBBBBB<BBBB'B<B<BB<B<BBBBBBBBBBBBB<BBBBBBBBBB<BB<<<<<
@ERR894725.100000002/2
TGTACAGAAAGTATCTGTTTTATTAATTCAACACTGTAAACATTTGCCATGTCTTAAAGTTCAAGTCTCTGTAACCTCCTGAGGTCAGGAGTTTGACATCAGCCTGACCAACATGGTGAAACCCTG
+
<<<<<B<BB<<<<<<B<<<<<<<<<<<<<B<<B<B<<<<<<B<<<<B<B<<<<B<<<<<<<<<<<B<B<B<B<B<BB<BB<BBBB<BBBBBBBB<B<B<<BBBBB<B<BBB<BBBBBBBBB<<<<<
@ERR894725.10000000/2
CTACATAATTTCCCTTACTGGACTAGCTTTTTTTCGGTGCTATTTTAAAATATGTTTAATTCTCCAACTGCTTAGAAGTCTTTTTAGATATTTTGCTGGGTGCAATAGCACACTTTATATCTATAT
+
7B<<B<<<<<B<<<B<<<B<<B<B<<BB<<<<<<<7<B<BB<<<<<<B<<<<<<B<<<<<<<<<<BB<B0<B0<0BBBB<BBBBB<<<<<<BBB<BB<<BB<BBB<BBB0<B<BBBBBBBBB<<7<
@ERR894725.100000003/1
AAAGAGAGCCCGCATTGCCAAGTCAATCCTAAGCCAAAAGAACAAAGCTGGAGGCATCACGCTACCTGACTTCAAATTATACTACAAGGCTACAGTAACCAAAACAGCATGGTACTGGTACCAAAA
+
<<<<B<B<<<<7BB<<<BBBBBB<BB<<BB<BBBBBBBBBBB<BBBBBB<BBBBBB<<B<7B<<<BB<B<BB<BBB<B<<<<B<<BB<BBBB<BBB<B<BBBBB<BBBBBBBBB<BBBBB<B<<<B
@ERR894725.100000003/2
AAATTTTCTCCCATTCTGTAGGTTGCCTGTTCACTCTGATGGTAGTTTCTTTTGCTGTGCAGAAGCTCTTTAGTTTAATTAGAACCCATTTGTCAATTTTGGCTTTTGTTGCCATTGCTTTTGGTG
+
<<<<<<<<B<<<B<<<B<<<<<<<<B<B<<<<B<B<B<B<<<<<B<<<<B<<<<BB<<<BB<<<BBB<B<B<BBBB<B<B<BBB<BBB<BB<B<B<<BBB<BBB<BBBBB<B<0BBBBBBBB<<70

Why is there a gap between the @ERR894725.10000000/1 group and the @ERR894725.10000000/2 one? The ASCII code for / is smaller than that of any number so I don't see why @ERR894725.100000002/1 takes precedence over @ERR894725.10000000/2.

I see the same pattern across the entire file. There is a gap between @ERR894725.10000001/1 and @ERR894725.10000001/1 filled by @ERR894725.100000012/1 and @ERR894725.100000012/1 and so on.

Any explanation for why this is happening and ways to avoid it is appreciated.

** Update: The actual data is quite large (1TB+) and can't be provided in it's entirety. I tried sorting the example above using the same commands again and it gives a correct result. I'll try resorting the already sorted complete output to see if it fixes things. Still weird why this has happened in first place.

| improve this question | | | | |
  • 2
    Hard to say why you're seeing the given output with those commands without seeing the input. – DopeGhoti Jun 11 '18 at 22:33
  • What happens if you remove the tr command? – ctrl-alt-delor Jun 11 '18 at 22:40
  • 2
    Try to post a minimal failing example. i.e. remove all the extra characters, and simplify your input file to contain only the parts that fail. It will be much easier to troubleshoot. – Sparhawk Jun 11 '18 at 23:33
  • 1
    Just a hunch: I wonder if your input file uses a different line terminator than the utilities expect. Do you get the same results if you use tr -s '\r' '\n' < reads.fq (to convert all linebreaks to Unix style) instead of cat reads.fq? – Gaultheria Jun 11 '18 at 23:53
  • 3
    If this is GNU sort, you can use --debug to see exactly what it's sorting on. – muru Jun 12 '18 at 1:30
3

Locale locale locale

You get unexpected results because your default locale is something other than C or C.UTF-8.

sort sorts in the order of byte values only in the C locale; in other locales sort sorts in the mysterious order defined by well-intentioned souls in the locale definition files.

Compare:

$ <<< $'@ERR894725.10000000/1\t1\n@ERR894725.100000002/1\t2\n@ERR894725.100000002/2\t3\n@ERR894725.10000000/2\t4' \
LC_ALL=en_US.UTF-8 sort -s -k1.1
@ERR894725.10000000/1   1
@ERR894725.100000002/1  2
@ERR894725.100000002/2  3
@ERR894725.10000000/2   4

$ <<< $'@ERR894725.10000000/1\t1\n@ERR894725.100000002/1\t2\n@ERR894725.100000002/2\t3\n@ERR894725.10000000/2\t4' \
LC_ALL=C.UTF-8 sort -s -k1.1
@ERR894725.10000000/1   1
@ERR894725.10000000/2   4
@ERR894725.100000002/1  2
@ERR894725.100000002/2  3

If you value consistency, always say LC_ALL=C sort (or LC_ALL=C.UTF-8 sort, or in general LC_ALL=<your-preferred-locale> sort) and never just sort.

The manual page has this to say:

*** WARNING *** The locale specified by the environment affects sort order. Set LC_ALL=C to get the traditional sort order that uses native byte values.

| improve this answer | | | | |
  • Thanks for your answer. This seems like it could be the source of the problem. I'll need to re-sort the file which might take a day or two and update the result. – DarthPaghius Jun 12 '18 at 19:48
  • So I ran LC_ALL=C cat reads.fq | paste - - - - | sort -k1,1 -S 3G | tr '\t' '\n' > sorted_reads.fq and the result is the same. Can it be due to setting the locale not immediately before calling sort? – DarthPaghius Jun 13 '18 at 18:22
  • LC_ALL=C sort. Or export LC_ALL=C and then < reads.fq paste - - - - | sort. But LC_ALL=C cat | paste | sort does not set LC_ALL for sort. (Explanation: LC_ALL is an environment variable. You can either set such a variable as a separate command and then it will take effect in all subsequent commands, or you can set it temporarily for just one command.) – AlexP Jun 13 '18 at 18:31
  • I see. I had assumed LC_ALL=C <command> is the equivalent of export LC_ALL=C; <command> . – DarthPaghius Jun 13 '18 at 18:37
  • It is, but it takes effect only for the <command>; and it must be a "simple command", so in LC_ALL=C cat | paste - - - - | sort it takes effect for cat but not for paste or sort. – AlexP Jun 13 '18 at 18:43

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