SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
BBMap (aligner for DNA/RNAseq) is now open-source and available for download. Brian Bushnell Bioinformatics 533 Yesterday 10:55 PM
BBMap for BitSeq dietmar13 Bioinformatics 1 04-30-2015 08:40 AM
Please help my BBMap cannot remove Illumina adapter TofuKaj Bioinformatics 4 04-28-2015 08:53 AM
BBMap Error Phage Hunter Bioinformatics 5 01-14-2015 04:34 AM
Introducing BBMap, a new short-read aligner for DNA and RNA Brian Bushnell Bioinformatics 24 07-07-2014 09:37 AM

Reply
 
Thread Tools
Old 04-04-2017, 05:58 AM   #141
Dinab
Junior Member
 
Location: Israel

Join Date: Apr 2017
Posts: 3
Default Quality distribution graph

A bit late, but hope you can help me:
I want to create an average quality distribution graph of multiple sequences from fastq file. Is it possible with reformat?
Also, if I use reformat.sh with maq=XX, does it simply calculates the avvarge quality based on base quality and checks the error rate (2% max)?

Thank you!
Dinab is offline   Reply With Quote
Old 04-04-2017, 04:09 PM   #142
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,637
Default

Quote:
Originally Posted by Dinab View Post
A bit late, but hope you can help me:
I want to create an average quality distribution graph of multiple sequences from fastq file. Is it possible with reformat?
Yes, the histogram section of the help covers that. Depending on how you want to plot it, you can use qhist, qchist, aqhist, or bqhist. The most generic is aqhist, which plots the number of reads per average quality score.

Quote:
Also, if I use reformat.sh with maq=XX, does it simply calculates the average quality based on base quality and checks the error rate (2% max)?
"maq" transforms the qualities into probability space, calculates the average error rate, and then transforms that back into the phred scale. So, a 2% maximum error rate would be approximately maq=17.
Brian Bushnell is offline   Reply With Quote
Old 04-09-2017, 05:11 AM   #143
Dinab
Junior Member
 
Location: Israel

Join Date: Apr 2017
Posts: 3
Talking

Perfect, Thank you!
Dinab is offline   Reply With Quote
Old 04-09-2017, 11:26 PM   #144
Dinab
Junior Member
 
Location: Israel

Join Date: Apr 2017
Posts: 3
Default Another question

Hi Brian,
Your previous response worked perfectly, thank you again.
I now have additional complexity to add: is it possible to get a plot of the mean quality of each sequence as variable of the length?
Dinab is offline   Reply With Quote
Old 04-10-2017, 08:51 AM   #145
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,637
Default

Quote:
Originally Posted by Dinab View Post
I now have additional complexity to add: is it possible to get a plot of the mean quality of each sequence as variable of the length?
No, I've never implemented something like that... the closest you could get with the existing package is to use a loop for each length with something like...

Code:
reformat.sh in=reads.fq minlen=147 maxlen=147 out=stdout.fq int=f | reformat.sh in=stdin.fq aqhist=aqhist_147.txt
Brian Bushnell is offline   Reply With Quote
Old 04-10-2017, 04:33 PM   #146
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,405
Default

Quote:
Originally Posted by Dinab View Post
is it possible to get a plot of the mean quality of each sequence as variable of the length?
FastQC gives you an average per sequence quality score.
GenoMax is offline   Reply With Quote
Old 05-07-2017, 10:25 PM   #147
blsfoxfox
Member
 
Location: auburn

Join Date: Jan 2013
Posts: 12
Default fastq format error

Hi Brian,

I found there is some format error with my fastq. By using reformat.sh in=in.fastq out=out.fastq, I get following error:

Input is being processed as unpaired
java.lang.AssertionError:
Error in in.fastq, line 4241767, with these 4 lines:
>m150430_235943_42146_c100804572550000001823173810081565_s1_p0/108988/0_3288 RQ=0.868
CTCTTTCCGGTAACACGCGCCTAGATACTAACAACCTACCTTATACTAATCAGCATGCCTAATACCAACA
ACCTACTTTATACCTAATAACAATGCCTTGATAACTAACAACATACCTTAATACCTAATACACGCTTAAT
ACCTAACAACAATACCTTAATACCTAACAACATACTTTAATACCTAACAACCTACCTTAATAACTTAATA

at stream.FASTQ.quadToRead(FASTQ.java:779)
at stream.FASTQ.toReadList(FASTQ.java:710)
at stream.FastqReadInputStream.fillBuffer(FastqReadInputStream.java:110)
at stream.FastqReadInputStream.nextList(FastqReadInputStream.java:95)
at stream.ConcurrentGenericReadInputStream$ReadThread.readLists(ConcurrentGenericReadInputStream.java:656)
at stream.ConcurrentGenericReadInputStream$ReadThread.run(ConcurrentGenericReadInputStream.java:635)
Input: 1060400 reads 9378503777 bases
Output: 1060400 reads (100.00%) 9378503777 bases (100.00%)

Time: 75.193 seconds.
Reads Processed: 1060k 14.10k reads/sec
Bases Processed: 9378m 124.73m bases/sec
Exception in thread "main" java.lang.RuntimeException: ReformatReads terminated in an error state; the output may be corrupt.
at jgi.ReformatReads.process(ReformatReads.java:1130)
at jgi.ReformatReads.main(ReformatReads.java:45)

Could you please provide any suggestions on which steps shoud I follow to find the exact error of that read?

Thanks!
blsfoxfox is offline   Reply With Quote
Old 05-07-2017, 10:33 PM   #148
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,637
Default

The problem here is that you are using a fasta sequence named "in.fastq". BBTools is sensitive to filenames, and *.fastq will be processed as a fastq file. If the file has no extension it will usually look at the contents to try to figure out what it is, but when there is a known extension, it assumes it is correct.

So, just rename the input file to "in.fasta" or add the flag "extin=.fasta" to override filename. Although for most uses of PacBio data I do recommend that you go back and get the original fastq file and use that, because the quality scores are often useful.

Last edited by Brian Bushnell; 05-07-2017 at 10:35 PM.
Brian Bushnell is offline   Reply With Quote
Old 05-07-2017, 10:42 PM   #149
blsfoxfox
Member
 
Location: auburn

Join Date: Jan 2013
Posts: 12
Default

Hi Brian,

Thanks for your quick response!

Actually, I am doing that on purpose. I found there is format error in my fastq file while using another software, but I don't know where is it. By using reformat.sh with in=in.fastq out=out.fastq, bbmap will report where is that error. In this case, read >m150430_235943_42146_c100804572550000001823173810081565_s1_p0/108988/0_3288 RQ=0.868 causes a corruption of the output. However, I still don't know what's wrong and the only solution I can think of now is using filterbyname tool to exclude that read from my fastq.
Thus, I really appreciate if you could provide any suggestion on how could I identify and fix that format error.

Thanks again for your time,
blsfoxfox is offline   Reply With Quote
Old 05-08-2017, 06:39 AM   #150
blsfoxfox
Member
 
Location: auburn

Join Date: Jan 2013
Posts: 12
Default

Quote:
Originally Posted by Brian Bushnell View Post
The problem here is that you are using a fasta sequence named "in.fastq". BBTools is sensitive to filenames, and *.fastq will be processed as a fastq file. If the file has no extension it will usually look at the contents to try to figure out what it is, but when there is a known extension, it assumes it is correct.

So, just rename the input file to "in.fasta" or add the flag "extin=.fasta" to override filename. Although for most uses of PacBio data I do recommend that you go back and get the original fastq file and use that, because the quality scores are often useful.
Hi Brian,

Sorry for that I should double checked my input.fastq. After grep multiple lines after that read, I find some fasta sequences in that fastq file. As they are raw sequences, I don't have a 'clean' backup for it. So my problem now is to find a tool which could extract fasta sequences from fastq.

Thanks,
blsfoxfox is offline   Reply With Quote
Old 05-08-2017, 06:47 AM   #151
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,405
Default

Code:
reformat.sh in=seq.fq out=seq.fa
will convert fastq sequences to fasta. If your problem is malformed fastq records then you need to find and delete those records manually.
GenoMax is offline   Reply With Quote
Old 05-08-2017, 06:58 AM   #152
blsfoxfox
Member
 
Location: auburn

Join Date: Jan 2013
Posts: 12
Default

Quote:
Originally Posted by GenoMax View Post
Code:
reformat.sh in=seq.fq out=seq.fa
will convert fastq sequences to fasta. If your problem is malformed fastq records then you need to find and delete those records manually.
Hi,

Thanks for your response!

My problem is that my file is mixed of fasta and fastq, and I don't know how to identify those fasta sequences and extract them.
blsfoxfox is offline   Reply With Quote
Old 05-08-2017, 07:02 AM   #153
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,405
Default

How did that happen?

Unless you have a defined reason I would suggest not trusting that file.
GenoMax is offline   Reply With Quote
Old 05-30-2017, 01:29 AM   #154
skatrinli
Junior Member
 
Location: turkey

Join Date: May 2017
Posts: 1
Default

Hey Brian,

Thanks for this great tool but I could not properly download it. I installed the latest version 37.25 but when i try to test the installation with the command
$ (installation directory)/stats.sh in=(installation directory)/resources/phix174_ill.ref.fa.gz
it comes this error:
Exception in thread "main" java.lang.RuntimeException: Unknown parameter Downloads/bbmap/resources/phix174_ill.ref.fa.gz
at jgi.AssemblyStats2.<init>(AssemblyStats2.java:166)
at jgi.AssemblyStats2.main(AssemblyStats2.java:39)

What did I do wrong?
skatrinli is offline   Reply With Quote
Old 05-30-2017, 04:26 AM   #155
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,405
Default

There is no "installation" needed for BBMap. Just uncompress and run (as long as you have Java available for your OS). Are you using Java 1.7 or 1.8?
@Brian no longer tests against Java 1.6 (which is what you may be using) if I recall.

I get the following when I run stats.sh.

Code:
$ stats.sh in=/path_to/bbmap/resources/phix174_ill.ref.fa.gz 
A       C       G       T       N       IUPAC   Other   GC      GC_stdev
0.2399  0.2144  0.2326  0.3130  0.0000  0.0000  0.0000  0.4471  0.0000

Main genome scaffold total:             1
Main genome contig total:               1
Main genome scaffold sequence total:    0.005 MB
Main genome contig sequence total:      0.005 MB        0.000% gap
Main genome scaffold N/L50:             1/5.386 KB
Main genome contig N/L50:               1/5.386 KB
Main genome scaffold N/L90:             1/5.386 KB
Main genome contig N/L90:               1/5.386 KB
Max scaffold length:                    5.386 KB
Max contig length:                      5.386 KB
Number of scaffolds > 50 KB:            0
% main genome in scaffolds > 50 KB:     0.00%


Minimum         Number          Number          Total           Total           Scaffold
Scaffold        of              of              Scaffold        Contig          Contig  
Length          Scaffolds       Contigs         Length          Length          Coverage
--------        --------------  --------------  --------------  --------------  --------
    All                      1               1           5,386           5,386   100.00%
   5 KB                      1               1           5,386           5,386   100.00%
GenoMax is offline   Reply With Quote
Old 05-30-2017, 09:02 AM   #156
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,637
Default

Quote:
Originally Posted by skatrinli View Post
Hey Brian,

Thanks for this great tool but I could not properly download it. I installed the latest version 37.25 but when i try to test the installation with the command
$ (installation directory)/stats.sh in=(installation directory)/resources/phix174_ill.ref.fa.gz
it comes this error:
Exception in thread "main" java.lang.RuntimeException: Unknown parameter Downloads/bbmap/resources/phix174_ill.ref.fa.gz
at jgi.AssemblyStats2.<init>(AssemblyStats2.java:166)
at jgi.AssemblyStats2.main(AssemblyStats2.java:39)

What did I do wrong?
You can't have spaces in the filenames without specific countermeasures like quotes. For example:

Code:
stats.sh in=foo bar.fa
Exception in thread "main" java.lang.RuntimeException: Unknown parameter bar.fa
        at jgi.AssemblyStats2.<init>(AssemblyStats2.java:166)
        at jgi.AssemblyStats2.main(AssemblyStats2.java:39)
That doesn't work...
Code:
stats.sh in="foo bar.fa"
Exception in thread "main" java.lang.RuntimeException: Unknown parameter bar.fa
        at jgi.AssemblyStats2.<init>(AssemblyStats2.java:166)
        at jgi.AssemblyStats2.main(AssemblyStats2.java:39)
That doesn't work either.

Code:
stats.sh in="foo\ bar.fa"
A       C       G       T       N       IUPAC   Other   GC      GC_stdev
NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN     0.0000

Main genome scaffold total:             0
Main genome contig total:               0
Main genome scaffold sequence total:    0.000 MB
Main genome contig sequence total:      0.000 MB        NaN% gap
Main genome scaffold N/L50:             0/0
Main genome contig N/L50:               0/0
Main genome scaffold N/L90:             0/0
Main genome contig N/L90:               0/0
Max scaffold length:                    0
Max contig length:                      0
Number of scaffolds > 50 KB:            0
% main genome in scaffolds > 50 KB:     0.00%


Minimum         Number          Number          Total           Total           Scaffold
Scaffold        of              of              Scaffold        Contig          Contig
Length          Scaffolds       Contigs         Length          Length          Coverage
--------        --------------  --------------  --------------  --------------  --------
That does work (though I ran it on an empty file).

The exact way to deal with spaces is system-specific. In the normal Windows shell you can just use quotes; in Linux bash it looks like you need quotes and an escape character (backslash); in Windows under a Linux emulator I'm not entirely sure. The easiest thing to do is to put files in a path that does not have any spaces (so, not in My Documents, but in C:\data\ or something like that.)
Brian Bushnell is offline   Reply With Quote
Old 05-31-2017, 06:12 PM   #157
TomHarrop
Member
 
Location: New Zealand

Join Date: Jul 2014
Posts: 17
Default

Hi Brian,

I'm using reformat.sh to play around with some Nanopore reads. Is there any way to get the histograms (e.g. mhist, qhist, bhist) to track longer reads, like the `max` parameter in readlength.sh?

Thanks,

Tom
TomHarrop is offline   Reply With Quote
Old 06-01-2017, 04:13 PM   #158
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,637
Default

Quote:
Originally Posted by TomHarrop View Post
Hi Brian,

I'm using reformat.sh to play around with some Nanopore reads. Is there any way to get the histograms (e.g. mhist, qhist, bhist) to track longer reads, like the `max` parameter in readlength.sh?

Thanks,

Tom
Sorry, the max lengths are currently constants. I'll add support for changing them. Generally I didn't find them all that useful for variable-length reads since I kind of designed them to find position-related anomalies, but it's fairly easy to change.
Brian Bushnell is offline   Reply With Quote
Old 08-12-2017, 12:29 AM   #159
jazz710
Member
 
Location: Iowa

Join Date: Oct 2012
Posts: 40
Default

Re: reformat.sh

Is there now, or could there be in the future, be a way to specify multiple sets of paired read inputs (ie. different libraries) which could be randomly sampled and output to a single FASTQ file?

Ie) in=FASTQ_A_R1,FASTQ_B_R1,FASTQ_C_R1 in2=FASTQ_A_R2,FASTQ_B_R2,FASTQ_C_R2 out=Subset_R1.fastq out2=Subset_R2.fastq

Can do it via a previous cat command (A+B+C -> reformat) but with large files cat can be an I/O issue.

Best and Thanks,
Bob
jazz710 is offline   Reply With Quote
Old 08-12-2017, 07:47 AM   #160
SNPsaurus
Registered Vendor
 
Location: Eugene, OR

Join Date: May 2013
Posts: 414
Default

Quote:
Originally Posted by jazz710 View Post
Re: reformat.sh

Is there now, or could there be in the future, be a way to specify multiple sets of paired read inputs (ie. different libraries) which could be randomly sampled and output to a single FASTQ file?

Ie) in=FASTQ_A_R1,FASTQ_B_R1,FASTQ_C_R1 in2=FASTQ_A_R2,FASTQ_B_R2,FASTQ_C_R2 out=Subset_R1.fastq out2=Subset_R2.fastq

Can do it via a previous cat command (A+B+C -> reformat) but with large files cat can be an I/O issue.

Best and Thanks,
Bob
I take a subset from each file and then cat the subsets together. That would be faster than combining the inputs with cat, at least. If the subset size is a significant fraction of most of the files it would still be slow, but if you are just collecting a small fraction of reads it is fast.
__________________
Providing nextRAD genotyping services. http://snpsaurus.com
SNPsaurus is offline   Reply With Quote
Reply

Tags
bbmap

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 12:53 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO