Max length reads from fastq

Hello everyone,

I want to extract reads with maximum length from a fastq file, I've opted for these commands and I'm not sure if my syntax is correct.

zcat file.fastq.gz | awk 'BEGIN {max_len=0} {if(NR%4==2) {len=length($0); if(len>max_len) {max_len=len; max_seq=$0; header=prev_header}}; prev_header=$0} END {print header; print max_seq; getline; getline; print; getline; print}' | gzip > reads_max_length.fastq.gz

Also, when I align with minimap2, I don't get any results. Maybe I did something wrong? Thanks

fastq

• 2,447 views • link

updated 2.4 years ago by Michael Dondrup

• written 2.4 years ago by marco.barr •

1 answer

This code gets you the first maximum length entry, but only half of it. You should test your code

$ cat test.fastq
@Areadlength 4
ACGT
+
!!!!
@Areadlength 3
ACG
+
!!!
@Areadlength 5
ACGTt
+
!!!!!
@Anotherreadlength 5
NNNNN
+
!!!!!



$ cat test.fastq | awk 'BEGIN {max_len=0} {if(NR%4==2) {len=length($0); if(len>max_len) {max_len=len; max_seq=$0; header=prev_header}}; prev_header=$0} END {print header; print max_seq; getline; getline; print; getline; print}'

@Areadlength 5
ACGTt

END {print header; print max_seq; getline; getline; print; getline; print}: At the END is at the end of the file, so getline is empty.

If you want the complete record, you need to save the + and quality line as well. If you want all records of maximum length you likely have to do two passes through your file.

• 0 views • link

written 2.4 years ago by Michael Dondrup

You're right, so you think a modification to my code could be this:

awk 'BEGIN {max_len=0; delete max_indices} {if(NR%4==2) {len=length($0); if(len>max_len) {max_len=len; delete max_indices; max_indices[NR]}}}; {if(length($0)==max_len && (NR-2) in max_indices) {max_indices[NR-2]}}; END {for(i in max_indices) {start=i*4-3; end=i*4; for(j=start; j<=end; j++) {print $(j)}}}' test.fastq > output.fastq

In the first step, find the maximum length and store the indices of the rows corresponding to this length. In the second step, use these indices to print all the records that have the maximum length. Do you think it's correct? Thank you very much

• 0 views • link

written 2.4 years ago by marco.barr •

The output from my test case is this, so no.

37: 40:

Honestly, I am not a big fan of awk. I am assuming it is possible but awk tends to become unreadable with increasing complexity.

• 0 views • link

written 2.4 years ago by Michael Dondrup

Me too.. What can I do instead of awk?

• 0 views • link

written 2.4 years ago by marco.barr •

I think seqkit should be safe. In a first round determine max size, in a second call extract all sequences of that length. If you MUST automate it:

M=$(seqkit stats -T test.fastq | tail -n1 | cut -f 8); seqkit seq -m $M test.fastq # edit: sorry, missing semicolon

@Areadlength 5
ACGTt
+
!!!!!
@Anotherreadlength 5
NNNNN
+
!!!!!

Unlike an awk hack, this is more fool-proof (not necessarily faster).

• 0 views • link

written 2.4 years ago by Michael Dondrup

Thanks! runnig this commands give me an error: Error: invalid argument "/home/PI4KA_ATG_50R_RP.fastq" for "-m, --min-len" flag: strconv.ParseInt: parsing "/home/PI4KA_ATG_50R_RP.fastq": invalid syntaxe

Maybe Seqkit seq interprets the values for the minimum length limit as integers, but instead, it's getting the FASTQ file path. To solve this issue, do I need to provide an integer for the minimum length limit?"

• 0 views • link

written 2.4 years ago by marco.barr •

M=$(seqkit stats -T test.fastq | tail -n1 | cut -f 8); seqkit seq -m $M test.fastq -o result.fastq.gz

seqkit seq -m $(seqkit stats -T test.fastq | tail -n1 | cut -f 8) test.fastq -o result.fastq.gz

• 0 views • link

written 2.4 years ago by shenwei356

Thanks! It worked

• 0 views • link

written 2.4 years ago by marco.barr •