How do I determine read length from a user-inputted fastq file and then input that information into a tool's command line?

Question

Can you guide me on how could take an input fastq file, read the length from read and then feed it into the findtail's command?

I have suggested concerns about documentation and program defaults to my prof and asked him to switch to cutadapt instead of using findtail that has no documentation. But he refuses so I have no option but to proceed with what he wants.

I want the user to input a RNA-seq data fastq file. Since the first line of each file contains the sequence length information, such as described:

@SRR3003130.1 HWI-ST741:596:H5FMVBCXX:1:1101:1244:1952 length=51

How do I pick the length=51, the 51 from this and feed in into the following command's --seqlength parameter:

    ./findtail --input_file Rrp6_S38_all_R1_001_egal.fastq --seqlength 51 \
--endgap 5 --taillength 10 --identity 100 --ptype A --stype T \
--output_format fasta --output_type pr > pr.fasta

The tool's command that I'm using only requires us to enter --seqlength once. So I do not think that determining the length of each and every read in the fastq file is of any use. Just simply the length from the first line of first read would suffice.

If something is lacking from my question, please let me know how can I improve.

Answer 1

I think that you cannot rely on that information being in a fastq header. Note also that each read has its own header! This is not just the first line! For example when I google "sample fastq file" the first examples do not have it. Spec doesn't seem to say anything about it.

For a one-off on a specific dataset this might be adequate:

grep -Eo "length=([0-9]+)" test.fq | cut -d"=" -f2
# prints "75" for each occurrence, so if you trust the first line you could pass to head:
grep -Eo "length=([0-9]+)" test.fq | cut -d"=" -f2 | head -1
# or you could explicitly sort and grab the longest read len:
grep -Eo "length=([0-9]+)" test.fq | cut -d"=" -f2 | sort | tail -1

For a more general application I'd suggest using some existing/explicit tools for measuring read length like seqkit or awk, rather than trusting some information that may or may not be in the header.

I can't find any docs for findtail, depending on what your application is I'd suggest using cutadapt instead, which does not require a single sequence length input.

Answer 2

To get the length of the first sequence from a single-line-per-sequence fastq, the second line can be piped through wc -c, and stored in a shell variable:

seqLength=$( head -n 2 input_file_name.fastq | tail -n 1 | wc -c)

This variable can then be used for the program arguments using ${variable}:

seqLength=$(head -n 2 input_file_name.fastq | tail -n 1 | wc -c)
./findtail --input_file input_file_name.fastq --seqlength ${seqLength} \
    --endgap 5 --taillength 10 --identity 100 --ptype A --stype T \
    --output_format fasta --output_type pr > pr.fasta

If the input file is gzipped [I have no idea if findtail supports gzipped files], then the length can be determined by piping through zcat:

seqLength=$( zcat input_file_name.fastq.gz | head -n 2 | tail -n 1 | wc -c)

As mentioned by many comments, note that some FASTQ files won't have the same length for each read, and it's unclear how findtail will behave for such a file.