0

I have to retrieve a specific information with known length from my fastq file and add it in another position.
For example given as input the following fastq file:

@SRR5394526.1 1 length=150  
CGATGTTAAATCAACGATAACTACACCG    
+SRR5394526.1 1 length=150  
AA<AFJFJJJJJJJJJJAJJJJJJJJJF    

I would like as output:

@SRR5394526.1.CGATGT 1 length=150    
TAAATCAACGATAACTACACCG    
+SRR5394526.1.CGATGT 1 length=150    
FJJJJJJJJJJAJJJJJJJJJF  

As you can notice the first 6 nucleotides were removed from either the sequence in the second line but also the sequence in the 4th line and appended after the first number 1 in the 1st and the 3rd lines. I have millions of chunk of this size (4 rows) within the file and this is just an example.

I have already found how to add/append info in a file: sed 's/myinfo/&,/4' and how to remove info in a file sed -e '423s!//!!; 424s!printf!//&!' but it's not enough. Any idea is much appreciated.

3

Using awk:

awk '(FNR-1) % 2 == 0 { name=$1; chr=$2; len=$3; next }
     (FNR-2) % 4 == 0 { seq=substr($0,1,6) }
                      { print name "." seq, chr, len
                        print substr($0,7) }' file.fastq >newfile.fastq

This awk program has three blocks.

  1. The first block triggers on every second row (the sequence and the quality data header lines), starting on the first one. It saves the three bits of information on that row into three variables. It then skips to the next line of input immediately.

  2. The second block extracts the 6 first characters from the sequence line into seq, but only for ever fourth row starting at row 2 (the sequence lines only).

  3. The last block runs only on rows not handled by the first block (every sequence or quality data line) and constructs the output.

To use this on gzip-compressed files (or bgzip-compressed files, commonly used in bioinformatics projects), use

zcat file.fastq.gz | awk '...' | bgzip -c >newfile.gz

To use a variable as the value used for cutting, consider

awk -v n=6 '(FNR-1) % 2 == 0 { name=$1; chr=$2; len=$3; next }
            (FNR-2) % 4 == 0 { seq=substr($0,1,n) }
                             { print name "." seq, chr, len
                               print substr($0,n+1) }'

Where the -v n=6 controls the length of the cut.

You may also put the actual awk code (everything inside the single quotes) into its own script file and use it as

awk -v n=6 -f script.awk file.fastq
0

data in fastq file, use gnu sed on 4 by 4 lines,

$ sed -nE ' N;N;N;s/(.+\.1)(\s.+\n)(.{6})(\w+)\s*(\n.+\.1)(.+\n).{6}(\w+)/\1.\3\2\4\5.\3\6\7/p' fastq

@SRR5394526.1.CGATGT 1 length=150
TAAATCAACGATAACTACACCG
+SRR5394526.1.CGATGT 1 length=150
FJJJJJJJJJJAJJJJJJJJJF

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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