Introducing Reformat, a fast read format converter

[darthsequencer]

Percent identity filter?

Hi is there a way to filter sam/bam filters by percent identity?

[Brian Bushnell]

Hi Brian,

I'm using the Mira assembler to assemble Illumina MiSeq reads I got from NCBI/SRA. I use "fastq-dump --split-files -F xxxx.sra" to extract the reads and get two files. Each set of paired end reads have exactly the same name, and mira needs the /1 and /2 added onto the reads and gives an error if it detects reads with the same name.

Thanks for your quick reply,
Rod

Rod,

You can now use the flags "addslash=t slashspace=f" together to accomplish that.

As an unrelated note, reformat now supports the "stoptag" flag, so it can process a sam file and add the stop coordinate of the read to the optional tags, prefixed by "YS:i:".

[rodf]

Hi Brian,

I'm using the Mira assembler to assemble Illumina MiSeq reads I got from NCBI/SRA. I use "fastq-dump --split-files -F xxxx.sra" to extract the reads and get two files. Each set of paired end reads have exactly the same name, and mira needs the /1 and /2 added onto the reads and gives an error if it detects reads with the same name.

Thanks for your quick reply,
Rod

[Brian Bushnell]

Hi Rod,

There is currently no way to fix that. I specifically made it that way to replicate the name structure of real Illumina reads. For paired reads, many tools or formats (such as sam) require both to have exactly matching names, excluding anything after the first whitespace (such as a 1 or 2). Most aligners, therefore, trim everything after the first whitespace, though with BBMap you can disable this with the "keepnames" flag.

What are you doing that requires no space? I could add an option for that, but I'm interested in why.

[rodf]

addslash=t problem

Hi Brian,

I tried to use reformat to add /1 and /2 to paired read names and a space was added between the name and the slash and this does not work for the assembler.

i.e.
want name/1 but get name /1

how to fix this?

Thanks,
Rod

[Brian Bushnell]

Yep, that's the correct approach. The "reads=-1" flag is not necessary (that's the default), and "int=f" also gets forced automatically when you have dual input files.

[jpummil]

Hey Brian!

So, to subsample a set of PE reads to reduce overall file size (creating quick running data set for a workshop), this would suffice?

reformat.sh in1=x1.fq in2=x2.fq out1=y1.fq out2=y2.fq reads=-1 samplerate=0.1 int=f

It would: Keep parings intact, give me 1/10 the data overall, ensure no interleaving (though I expect assigning the pairings at the beginning would do this as well)

I already did a quick quality trimming with:
reformat.sh in1=x1.fastq in2=x2.fastq out1=y1.fastq out2=y2.fastq outsingle=singletons.fq qtrim
=rl trimq=10 minlength=50

Edit: Seems to have worked!

-rw-rw-r-- 1 jpummil jpummil 2.0G Feb 24 13:06 L001_R1_001_Qt.fastq
-rw-rw-r-- 1 jpummil jpummil 2.0G Feb 24 13:06 L001_R2_001_Qt.fastq

-rw-rw-r-- 1 jpummil jpummil 199M Feb 24 13:43 L001_R1_001_Sub.fastq
-rw-rw-r-- 1 jpummil jpummil 199M Feb 24 13:43 L001_R2_001_Sub.fastq

And it ran just fine in SPAdes.

[Brian Bushnell]

Hi Brian,

Does the subsampling tool randomly subsample reads? The program I am using now (Geneious) only takes the first specified % of reads and does not randomize.

Thank you,
Marisa

Yes, it randomly subsamples. You can set the RNG seed with "sampleseed=NUMBER" if you want deterministic random sampling; by default, the seed is random.

Reformat can also give the first X reads with the "reads=X" flag, or combine the reads=X and samplerate=Y to subsample a fraction of the first X reads, etc.

[Marisa_Miller]

Hi Brian,

Does the subsampling tool randomly subsample reads? The program I am using now (Geneious) only takes the first specified % of reads and does not randomize.

Thank you,
Marisa

[Brian Bushnell]

This was caused by the bases and qualities having different lengths. Either the fna and qual file do not go together, or their order is different, or one of the files is misformatted.

"reads.fna" should already be in fasta format, though, which is what you specified as the output format. If you want fastq, the output filename needs to end with "fastq".

I suggest you post the first sequence in the reads.fna and reads.qual so we can see what's going on.

[Nanu]

When I use the command reformat.sh in bbtools package I am getting the following error::
java -ea -Xmx200m -cp /home/himanshu/Downloads/me2/bbmap/current/ jgi.ReformatReads -in=reads.fna qfin=reads.qual out=reads.fasta
Executing jgi.ReformatReads [-in=reads.fna, qfin=reads.qual, out=reads.fasta]

Input is being processed as unpaired
Exception in thread "Thread-1" java.lang.AssertionError
at stream.FastaQualReadInputStream3.makeRead(FastaQualReadInputStream3.java:257)
at stream.FastaQualReadInputStream3.toReadList(FastaQualReadInputStream3.java:147)
at stream.FastaQualReadInputStream3.toReads(FastaQualReadInputStream3.java:113)
at stream.FastaQualReadInputStream3.fillBuffer(FastaQualReadInputStream3.java:97)
at stream.FastaQualReadInputStream3.hasMore(FastaQualReadInputStream3.java:56)
at stream.ConcurrentGenericReadInputStream$ReadThread.readLists(ConcurrentGenericReadInputStream.java:745)
at stream.ConcurrentGenericReadInputStream$ReadThread.run(ConcurrentGenericReadInputStream.java:737)

Please help me

[dpryan]

Consider this the first feature request then

[Brian Bushnell]

No, the conversion from sam/bam is stateless and 'dumb' - each line in the sam file will generate a FASTQ read, so secondary and chimeric alignments would cause problems.

[dpryan]

Looks nice. Out of curiousity, can it handle chimeric (i.e., non-linear) alignments when converting from BAM to fastq? That's been a real weak point of a lot of other tools.

[Brian Bushnell]

Reformat is a member of the BBMap/BBTools package. It is a multipurpose tool designed for converting reads or other nucleotide data between different formats. It supports, and can inter-convert:

fastq
fasta
fasta+qual
sam
scarf (an old Illumina format)
bam (if samtools is installed)
gzip
zip
ascii-33 (sanger)
ascii-64 (old Illumina)
paired files
interleaved files

It is multithreaded and can process data at over 500 megabytes per second, and can accept streams from standard in and write to standard out, allowing it to be easily dropped into the middle of a pipeline for format conversion. Reformat autodetects formats based on file extensions and content, making it very easy to use; and the autodetection can be overridden, allowing flexibility for people who don't like to follow naming conventions, or out-of-spec fastq files with qualities values like -17 or 120.

The program has been gradually expanded, and can now perform various other functions. None of these will break pairing, if the input is paired.

Quality trimming (either or both ends)
Quality filtering
Fixed-length trimming
Generation of histograms (base composition, quality, etc)
Subsampling (to a fraction of input reads, or an exact number of reads or bases)
Changing fasta line-wrapping length
Reverse-complementing (all reads or only read 2)
Adding /1 and /2 suffix to read names
GC-content filtering
Length-filtering
Testing for corrupted interleaved files

Reformat is compatible with any platform that supports Java 1.7 or higher. It also has a bash shellscript for simpler invocation. Typical usage examples:

Reformat fastq into fasta:
reformat.sh in=x.fq out=y.fa

Interleave paired reads:
reformat.sh in1=x1.fq in2=x2.fq out=y.fq

Note - you can actually use a shortcut if paired read files have the same name with a 1 and a 2. This is equivalent to the above command:
reformat.sh in=x#.fq out=y.fq

De-interleave reads:
reformat.sh in=x.fq out1=y1.fq out2=y2.fq

Verify that interleaving appears correct, assuming Illumina namimg conventions:
reformat.sh in=x.fq vint

Convert ASCII-33 to ASCII-64:
reformat.sh in=x.fq out=y.fq qin=33 qout=64

Quality-trim paired reads to Q10 on the left and right ends and discard reads shorter than 50bp after trimming:
reformat.sh in1=x1.fq in2=x2.fq out1=y1.fq out2=y2.fq outsingle=singletons.fq qtrim=rl trimq=10 minlength=50

Subsample 10% of the first 20000 pairs in an interleaved file:
reformat.sh in=x.fq out=y.fq reads=20000 samplerate=0.1 int=t
(in this case "int=t" overrides interleaving autodetection, to ensure reads are treated as pairs)

Pipe in a gzipped sam file and pipe out fasta:
reformat.sh in=stdin.sam.gz out=stdout.fa

Reverse-complement reads:
reformat.sh in=x.fq out=y.fq rcomp

For reformatting a file with very long sequences, Reformat will need more memory; just add the additional flag "-Xmx2g". For example, to change the line-wrapping length on the human genome (which has individual sequences over 200Mbp long) to 70 characters:
reformat.sh -Xmx2g in=HG19.fa.gz out=HG19_wrapped.fa.gz fastawrap=70

For additional functions, please run the shellscript with no arguments, or just read it with a text editor. If you have any questions, please post them in this thread.

For people using a non-bash terminal, you may need to type "bash reformat.sh" instead of just "reformat.sh".
For users of Windows or other platforms that do not support bash shellscripts, replace "reformat.sh" with "java -ea -Xmx200m /path/to/bbmap/current/ jgi.ReformatReads"
for example,
java -ea -Xmx200m C:\bbmap\current\ jgi.ReformatReads in=x.fq out=y.fa

Reformat can be downloaded with BBTools here:

https://sourceforge.net/projects/bbmap/