raw read processing: /1, /2 not always added/specified properly

[bluegenes]

sharing issue arising from troubleshooting with @jeanzzhao

tl;dr: Grist initial read processing is (silently) incorrect for samples that use .1 or .2 in read names (we assume this is used for paired end information). I'm not sure what issues arise downstream of the incorrect processing (if any) -- will update as we assess.

In rule download_sra_wc, grist attempts to standardize read names to the /1 and /2 expected/required by many programs:

echo processing R1...
       seqtk seq -C {output.temp_r1} | \
           perl -ne 's/\.([12])$/\/$1/; print $_' | \
           gzip -c > {output.r1} &

However, read names in some files may contain .1 or .2 as part of the name, e.g. raw reads from ERR257715_1.fastq:

@ERR257715.1 9L6V3M1:316:D0JM0ACXX:5:1101:1393:2121 length=101
TCACACGGATCTTACGCTCTTTATTCATCATGTTCATGAGAAAGTTCCGGAGATTCTTCGGGTTAGTCTCTTCGATGGGACGTTTGAACGTGTACCAGAGC
+ERR257715.1 9L6V3M1:316:D0JM0ACXX:5:1101:1393:2121 length=101
@?@DDEFF<FFHFIIIIIGIIIAGIHGIIGGGFHG>DEFFDBDFHGABGG?FBHICEHGIIGACHFFFFDFCED2;@BB?BB;@CB:@CBB@CCCC>CAB9
@ERR257715.2 9L6V3M1:316:D0JM0ACXX:5:1101:1430:2152 length=101
TACCCGTCCACGGTCCATCCCATACTTTTTCCCAGAGTGATGGACGTTCACTAATTAAAAAAATATTTAAGTGCTCACAATTTACAGTCCATGCAAGTAAA
+ERR257715.2 9L6V3M1:316:D0JM0ACXX:5:1101:1430:2152 length=101
@CCFFFFFDHHHHCGHIHIJJEICGHGGIJIGHGGHGFEFEHGGIJGGFHJJJIJJGGHIIHB;CCDEEECA@CDA>CCCDDDCDDDCCCCCC?CCDCDCC
@ERR257715.3 9L6V3M1:316:D0JM0ACXX:5:1101:1335:2158 length=101
AGTGATGGACGTTCACTAATTAAAAAAATATTTAAGTGCTCACAATTTACAGTCCATGCAAGTAAAATGGTCGCACGTCGTACTCTATTTCCATCTACTG
+ERR257715.3 9L6V3M1:316:D0JM0ACXX:5:1101:1335:2158 length=101
@CCFDFFFFHAHHIIIIIGIIFFHIIJIGJHIIIJHIIIJGGIJIGJJDCFHHIJIIIIIIG@EHGIGEHHFFB<9@DDDDDDDDACCDEDDECDDDDE

Grist processing by the code above yields this result:

@ERR257715/1
TCACACGGATCTTACGCTCTTTATTCATCATGTTCATGAGAAAGTTCCGGAGATTCTTCGGGTTAGTCTCTTCGATGGGACGTTTGAACGTGTACCAGAGC
+
@?@DDEFF<FFHFIIIIIGIIIAGIHGIIGGGFHG>DEFFDBDFHGABGG?FBHICEHGIIGACHFFFFDFCED2;@BB?BB;@CB:@CBB@CCCC>CAB9
@ERR257715/2
TACCCGTCCACGGTCCATCCCATACTTTTTCCCAGAGTGATGGACGTTCACTAATTAAAAAAATATTTAAGTGCTCACAATTTACAGTCCATGCAAGTAAA
+
@CCFFFFFDHHHHCGHIHIJJEICGHGGIJIGHGGHGFEFEHGGIJGGFHJJJIJJGGHIIHB;CCDEEECA@CDA>CCCDDDCDDDCCCCCC?CCDCDCC
@ERR257715.3
AGTGATGGACGTTCACTAATTAAAAAAATATTTAAGTGCTCACAATTTACAGTCCATGCAAGTAAAATGGTCGCACGTCGTACTCTATTTCCATCTACTGC
+
@CCFDFFFFHAHHIIIIIGIIFFHIIJIGJHIIIJHIIIJGGIJIGJJDCFHHIJIIIIIIG@EHGIGEHHFFB<9@DDDDDDDDACCDEDDECDDDDEC#

I've reprocessed the files for this sample by just adding /1 or /2 rather than substituting for .1/.2

for infile in ERR257715_1.fastq ERR257715_2.fastq
  do
    output="${infile%.fastq}.fq.gz"
    suffix="/${infile##*_}"
    suffix="${suffix%.fastq}"
    seqtk seq -C "$infile" | awk -v s="$suffix" '{ if (NR % 2 == 1 && /^@/) { print $0 s } else { print $0 } }' | gzip -c > "$output"
  done

which produces the following for the _1 file:

@ERR257715.1/1
TCACACGGATCTTACGCTCTTTATTCATCATGTTCATGAGAAAGTTCCGGAGATTCTTCGGGTTAGTCTCTTCGATGGGACGTTTGAACGTGTACCAGAGC
+
@?@DDEFF<FFHFIIIIIGIIIAGIHGIIGGGFHG>DEFFDBDFHGABGG?FBHICEHGIIGACHFFFFDFCED2;@BB?BB;@CB:@CBB@CCCC>CAB9
@ERR257715.2/1
TACCCGTCCACGGTCCATCCCATACTTTTTCCCAGAGTGATGGACGTTCACTAATTAAAAAAATATTTAAGTGCTCACAATTTACAGTCCATGCAAGTAAA
+
@CCFFFFFDHHHHCGHIHIJJEICGHGGIJIGHGGHGFEFEHGGIJGGFHJJJIJJGGHIIHB;CCDEEECA@CDA>CCCDDDCDDDCCCCCC?CCDCDCC
@ERR257715.3/1
AGTGATGGACGTTCACTAATTAAAAAAATATTTAAGTGCTCACAATTTACAGTCCATGCAAGTAAAATGGTCGCACGTCGTACTCTATTTCCATCTACTGC
+
@CCFDFFFFHAHHIIIIIGIIFFHIIJIGJHIIIJHIIIJGGIJIGJJDCFHHIJIIIIIIG@EHGIGEHHFFB<9@DDDDDDDDACCDEDDECDDDDEC#

This is not ideal replacement code however, since it doesn't handle cases where /1,/2 already exist, and doesn't substitute .1/.2 if they are used for specifying pair information.

I'm not sure what issues arise downstream of the incorrect processing (if any) -- will update as we assess.