Question: Extracting last N bases of reads in FASTQ file
1
gravatar for lebrigand
2.4 years ago by
lebrigand10
lebrigand10 wrote:

Hello, I can't find an easy solution to extract the last N bases of all reads from a FASTQ file, is there an easy solution with prinseq-lite, cutadapt or Fastx-toolkit to do it ? Or an other tool. I know i can do it on my own with either a java or a perl script, but well it seems so obvious that one of those program can do it in a few seconds, but i do not find the way to do it. thanks a lot to save me time, and thanks for this so resourcefull forum. kevin

rna-seq sequence • 1.7k views
ADD COMMENTlink modified 6 months ago by Harumi10 • written 2.4 years ago by lebrigand10
1

You mean remove them?

ADD REPLYlink written 2.4 years ago by Gabriel R.2.5k
1

Do you want to remove last N bases from a fastq file or want to extract them into new file?

ADD REPLYlink written 2.4 years ago by venu5.6k
2

Answers below cover both possibilities.

ADD REPLYlink written 2.4 years ago by genomax55k
1

I was puzzled by the same thing.

ADD REPLYlink written 2.4 years ago by Gabriel R.2.5k

Hi all, thanks for the answers but it seems not working as i want. First all the sed command does not output fastq file so it is useless for me. Concerning the bbduk options, i did not try bbduk yet but it really seems interesting and maybe it can do the stuff, but this command line does not work. Let me explain better what i want, here is the first reads of my input.fastq file:

@NB500XXXXXX

CGCTCNAGAACAGGGTTTGTTAAGAGTAC

+

AAAAA#EEEEEEEEEEEEEEE/EEEEEEE

What i want as output is the last 4 bases in a fastq format file:

@NB500XXXXXX

GTAC

+

EEEE

NB: cat input.fastq | /share/apps/local/bbmap/bbduk.sh in=stdin.fastq forcetrimright=3 out=output.fq minlength=3

--> output.fastq contains the first 4 bases not the last ones

Note that reads can have differents length, what i want is really the last 4 bases of all reads in my fastq file in a fastq format (including QV) thanks for helping me ! kevin.

ADD REPLYlink modified 2.4 years ago • written 2.4 years ago by lebrigand10
1

If future use the "Add comment"/"Add reply" options against the respective answers/comments to provide additional information. Do not add a "new answer".

As for the BBMap solution I suggest that you use reformat.sh like so

reformat.sh in=input.fastq forcetrimleft=[seq_length - bases you want to keep] out=right.fastq

Replace [seq_length - bases you want to keep] this part with a real number.

Edit: I see that reads can have different length so this solution would not work.

ADD REPLYlink modified 2.4 years ago • written 2.4 years ago by genomax55k
1

It would help if the original question gave the exact details of what you want. The way it is worded makes it seem like all you wanted was the last four bases of each read simply printed out. Truncating fastq entries down to the last 4bp is a different problem.

ADD REPLYlink written 2.4 years ago by pld4.7k
4
gravatar for pld
2.4 years ago by
pld4.7k
United States
pld4.7k wrote:
sed -n '2~4p' myfile.fq | grep -o '.\{5\}$'

This will print the last 5 characters. The sed command prints the sequences (every 4th line of the fastq, starting with the second), the grep grabs the last n characters.

ADD COMMENTlink written 2.4 years ago by pld4.7k
3

+1

You can also do it all in sed: sed -rn 's/.*(.{3})/\1/; 2~4p' myfile.fq for the last 3 characters, for example

ADD REPLYlink modified 2.4 years ago • written 2.4 years ago by george.ry1.1k
1

Note. This solution does not keep the output reads in fastq format.

ADD REPLYlink written 2.4 years ago by genomax55k
1

Note. This solution in current form does not keep the output reads in fastq format.

ADD REPLYlink written 2.4 years ago by genomax55k
1

OP didn't ask for a way to truncate fastq entries, just to pull the last n bp of each read out of a fastq file.

ADD REPLYlink written 2.4 years ago by pld4.7k
4
gravatar for John
2.4 years ago by
John12k
Germany
John12k wrote:

I don't have time to check it on a FASTQ file, but in python:

with open('input.fastq','r') as f:
    toggle = False
    for line in f:
        if toggle: print line[-5:-1]
        else:      print line[ 0:-1]
        toggle = not toggle

Change the path of input.fastq, save the code as dangerous.py, and run python dangerous.py > new.fastq

Works on reads of different lengths, and outputs the format you were looking for.

ADD COMMENTlink modified 2.4 years ago • written 2.4 years ago by John12k
2

Confirmed to work as advertised :-)

ADD REPLYlink written 2.4 years ago by genomax55k
2
gravatar for michael.ante
2.4 years ago by
michael.ante2.6k
Austria/Vienna
michael.ante2.6k wrote:

Hi Kevin,

give bbduk.sh from the bbmap suite a try. It will be something like:

cat myfile.fq | bbduk.sh in=stdin.fq forcetrimright=$N  out=stdout.fq minlength=$N

AFAIK, bbduk is 0-based; so if you want to have the last 5 nt, your N has to be 4.

Cheers,

Michael

ADD COMMENTlink written 2.4 years ago by michael.ante2.6k

I'm sorry. In my toy example it seemed to work.

If you want to make it complicated: Use reformat.sh to make the reverse complement, get the first N reads, and change it back with another reverse complement:

cat myfile.fq | reformat.sh in=stdin.fq  out=stdout.fq rcomp=t | bbduk.sh in=stdin.fq forcetrimright=$N  out=stdout.fq minlength=$N | reformat.sh in=stdin.fq  out=stdout.fq rcomp=t

Anyway, I'd go for John's Python solution.

ADD REPLYlink modified 2.4 years ago • written 2.4 years ago by michael.ante2.6k
1
gravatar for shenwei356
18 months ago by
shenwei3564.1k
China
shenwei3564.1k wrote:

I'm wondering why the No.1 user Biostar modified some threads and brought them to homepage these days.

But seqkit can do this easily.

N=100
gzip -d -c read_1.fq.gz | seqkit subseq -r -$N:-1 | gzip -c > read_1.t.fq.gz

The definition of region is 1-based and with some custom design.

 1-based index    1 2 3 4 5 6 7 8 9 10
negative index    0-9-8-7-6-5-4-3-2-1
           seq    A C G T N a c g t n
           1:1    A
           2:4      C G T
         -4:-2                c g t
         -4:-1                c g t n
         -1:-1                      n
          2:-2      C G T N a c g t
          1:-1    A C G T N a c g t n
          1:12    A C G T N a c g t n
        -12:-1    A C G T N a c g t n

Example:

$ echo -e "@seq\nACGTNacgtn\n+\nGGGGGCCCCC"                        
@seq
ACGTNacgtn
+
GGGGGCCCCC

first 6 bases:

$ echo -e "@seq\nACGTNacgtn\n+\nGGGGGCCCCC" | seqkit subseq -r 1:6  
@seq
ACGTNa
+
GGGGGC

last 6 bases:

$ echo -e "@seq\nACGTNacgtn\n+\nGGGGGCCCCC" | seqkit subseq -r -6:-1
@seq
Nacgtn
+
GCCCCC
ADD COMMENTlink modified 18 months ago • written 18 months ago by shenwei3564.1k
1

I'm wondering why the No.1 user Biostar modified some threads and brought them to homepage these days.

That is a programmed feature in Biostars that "bumps" old threads (I am sure there is some logic that @Istvan can comment on) to keep the main page periodically populated.It provides a chance for (relatively) new users like you to offer fresh takes on answers/comments.

ADD REPLYlink written 18 months ago by genomax55k
2

Since it's hard to find out what the robot modified. I tend to accept your guessing.

Yes, although I registered 4 years ago, I just started to comment few months ago :P

ADD REPLYlink written 18 months ago by shenwei3564.1k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 757 users visited in the last hour