SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Bowtie, an ultrafast, memory-efficient, open source short read aligner Ben Langmead Bioinformatics 513 05-14-2015 02:29 PM
Introducing BBMap, a new short-read aligner for DNA and RNA Brian Bushnell Bioinformatics 24 07-07-2014 09:37 AM
Miso's open source joyce kang Bioinformatics 1 01-25-2012 06:25 AM
Targeted resequencing - open source stanford_genome_tech Genomic Resequencing 3 09-27-2011 03:27 PM
EKOPath 4 going open source dnusol Bioinformatics 0 06-15-2011 01:10 AM

Reply
 
Thread Tools
Old 03-15-2018, 09:03 PM   #601
darthsequencer
Member
 
Location: california

Join Date: Feb 2012
Posts: 34
Default

Quote:
Originally Posted by GenoMax View Post
Have you tried ambig=all?
Yeah - I have but from what I can tell it maxes out at 6 reports unless maxsites is increased. I am wondering if there is a setting that avoids maxsites and reports everything above some score.
darthsequencer is offline   Reply With Quote
Old 03-16-2018, 09:41 AM   #602
SNPsaurus
Registered Vendor
 
Location: Eugene, OR

Join Date: May 2013
Posts: 450
Default

Quote:
Originally Posted by darthsequencer View Post
Yeah - I have but from what I can tell it maxes out at 6 reports unless maxsites is increased. I am wondering if there is a setting that avoids maxsites and reports everything above some score.
You might try the align2.BBMapPacBioSkimmer mapper with ambig=all and expectedsites= some high number. Not exactly what you were asking but I get dozens of hits with that.
__________________
Providing nextRAD genotyping and PacBio sequencing services. http://snpsaurus.com
SNPsaurus is offline   Reply With Quote
Old 03-21-2018, 11:26 AM   #603
SNPsaurus
Registered Vendor
 
Location: Eugene, OR

Join Date: May 2013
Posts: 450
Default

I have a question about the bbmap summary output. Here is an example:

Code:
Reads Used:             1061591 (151284618 bases)

Mapping:                15.250 seconds.
Reads/sec:              69611.57
kBases/sec:             9920.17


Read 1 data:            pct reads       num reads       pct bases          num bases

mapped:                  76.2432%          809391        77.6086%          117409909
unambiguous:             75.3885%          800318        76.8084%          116199316
ambiguous:                0.8547%            9073         0.8002%            1210593
low-Q discards:           0.0000%               0         0.0000%                  0

perfect best site:       11.4049%          121073         9.7282%           14717263
semiperfect site:        58.9222%          625513        60.0540%           90852423

Match Rate:                   NA               NA        99.1429%          116413381
Error Rate:              20.6634%          183594         0.2527%             296692
Sub Rate:                19.9279%          177059         0.2369%             278137
Del Rate:                 0.6244%            5548         0.0084%               9819
Ins Rate:                 0.6302%            5599         0.0074%               8736
N Rate:                  72.0298%          639985         0.6044%             709655
What does an N rate of 72% mean in this context? The reference has 0 bases that are N. There are 20,000 out of the 1 million reads that contain an N. The N Rate percent bases of 0.6% seems high compared to 20,000 reads with 1 or 2 Ns, but closer to what I think I see than the 72% N for reads.
__________________
Providing nextRAD genotyping and PacBio sequencing services. http://snpsaurus.com
SNPsaurus is offline   Reply With Quote
Old 03-21-2018, 03:13 PM   #604
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

This means 72 percent of the reads mapped with an "N" symbol in the match string, an internal data structure similar to a cigar string. The "N" symbol denotes either an N in the read or an N in the reference, which can include sequences that go off the end of scaffolds (the ends are internally padded with Ns). Additionally, degenerate IUPAC codes are counted as Ns.
Brian Bushnell is offline   Reply With Quote
Old 03-21-2018, 03:24 PM   #605
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Quote:
Originally Posted by kimng View Post
Hello, I've been using different parts of the BBMap suite and mixing them with other softwares for WGS of bacterial samples (primarily for the purpose of automated QC). I'm a fan of callvariants.sh as it works quickly and doesn't require me to sort .sam files before hand and it's simple to clearfilters to return all values. Unfortunately when using callvariants.sh with .sam files created from minimap2 (https://github.com/lh3/minimap2) using illumina nextseq reads. I receive the error seen below. I was curious if this was a minor bug that could be corrected?

I'm running callvariants.sh from bbmap v37.90 obtained from conda (https://anaconda.org/bioconda/bbmap). Appreciate any assistance given.

-Kim

Error:
java -ea -Xmx48779m -Xms48779m -cp /home/ssi.ad/kimn/.conda/envs/env_kim/opt/bbmap-37.90/current/ var2.CallVariants in=aln.sam ref=contigs.fasta vcf=minimap.vcf ploidy=1 clearfilters
Executing var2.CallVariants [in=aln.sam, ref=contigs.fasta, vcf=minimap.vcf, ploidy=1, clearfilters]

Loading reference.
Time: 0.304 seconds.
Processing input files.
Exception in thread "Thread-3" Exception in thread "Thread-4" Exception in thread "Thread-5" java.lang.AssertionError: TODO: Encountered a read with 'M' in cigar string but no MD tag.
at stream.SamLine.toShortMatch(SamLine.java:1200)
at stream.SamLine.toRead(SamLine.java:1972)
at stream.SamLine.toRead(SamLine.java:1832)
at stream.SamReadStreamer$ProcessThread.makeReads(SamReadStreamer.java:206)
at stream.SamReadStreamer$ProcessThread.run(SamReadStreamer.java:135)
java.lang.AssertionError: TODO: Encountered a read with 'M' in cigar string but no MD tag.
at stream.SamLine.toShortMatch(SamLine.java:1200)
at stream.SamLine.toRead(SamLine.java:1972)
at stream.SamLine.toRead(SamLine.java:1832)
at stream.SamReadStreamer$ProcessThread.makeReads(SamReadStreamer.java:206)
at stream.SamReadStreamer$ProcessThread.run(SamReadStreamer.java:135)
java.lang.AssertionError: TODO: Encountered a read with 'M' in cigar string but no MD tag.
at stream.SamLine.toShortMatch(SamLine.java:1200)
at stream.SamLine.toRead(SamLine.java:1972)
at stream.SamLine.toRead(SamLine.java:1832)
at stream.SamReadStreamer$ProcessThread.makeReads(SamReadStreamer.java:206)
at stream.SamReadStreamer$ProcessThread.run(SamReadStreamer.java:135)
The problem here is that minimap uses old-style cigar strings (M symbol instead of = and X) and also does not produce MD tags. I've added the ability to handle reads in that situation and it will be released in v37.95 (soon).

Prior to that you could probably add MD tags with "samtools calmd".
Brian Bushnell is offline   Reply With Quote
Old 03-21-2018, 03:26 PM   #606
SNPsaurus
Registered Vendor
 
Location: Eugene, OR

Join Date: May 2013
Posts: 450
Default

The reference in this case is a list of RAD loci sequences, each at 150 bp. The reads are 150 bp and should match one of the RAD loci. Since there are no N (or degenerate nucleotides) in the reference, and only a small number of reads with an N, it must have to do with the padding. A trimmed read would not have N padding, would it? This sample was degraded DNA, so plenty of reads were trimmed to less than 150 nucleotides. I doubt that 72% of the loci from this sample had an insertion that made the read extend past the reference (generated from consensus of many samples). Odd, but I guess it is not critical for me to fully figure this out.
__________________
Providing nextRAD genotyping and PacBio sequencing services. http://snpsaurus.com
SNPsaurus is offline   Reply With Quote
Old 03-29-2018, 06:28 AM   #607
tabwalsh
Junior Member
 
Location: Vienna

Join Date: Mar 2018
Posts: 2
Default Bias towards forward strand in msa.sh

Hello Brian,

Thanks for the lovely tool. I've been using msa.sh with cutprimers.sh (as described here) to extract sequences between pairs of primers.

However, I got some unexpectedly large sequences as output, and when I took a look at the SAM alignments output by msa.sh, none of them are mapping to the reverse strand.

For example, with the following template sequence…
>template
CTCGAGGGGCCTAGACATTGCCCTCCAGAGAGAGCACCCAACACCCTCCAGGCTTGACCGGCCAGGGTGTCCCCTTGACACCCGTCAAGCCTCACGA
…it should be possible to use the following commands to extract the sequence between the given primers:
msa.sh in=template.fa primer=GGGGCCTAGACATTGCCCTCCA out=fwd.sam
msa.sh in=template.fa primer=GACACCCTGGCCGGTCAAGCCT out=rev.sam
cutprimers.sh in=template.fa sam1=fwd.sam sam2=rev.sam include=t out=amplicon.fa
This outputs the following (primer alignments shown in uppercase):
GGGGCCTAGACATTGCCCTCCAgagagagcacccaacaccctccaggcttgaccggccagggtgtccccttGACACCCGTCAAGCCT
The sequence GACACCCGTCAAGCCT is the best match for the reverse primer on the forward strand, but not the best match overall.

That is on the reverse strand, which would result in the following (primer alignments again in uppercase):
GGGGCCTAGACATTGCCCTCCAgagagagcacccaacaccctccAGGCTTGACCGGCCAGGGTGTC
For the latest version of BBMap (BBMap_37.95.tar.gz), it looks like this issue might be fixed by changing line 183 of file 'bbmap/current/jgi/FindPrimers.java' to:
score2=ss.score=ss.quickScore;

Last edited by tabwalsh; 03-29-2018 at 06:31 AM.
tabwalsh is offline   Reply With Quote
Old 03-29-2018, 08:02 AM   #608
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,783
Default

Quote:
Originally Posted by tabwalsh View Post
For the latest version of BBMap (BBMap_37.95.tar.gz), it looks like this issue might be fixed by changing line 183 of file 'bbmap/current/jgi/FindPrimers.java' to:
score2=ss.score=ss.quickScore;
I suggest that you create a ticket on BBMap SF site since you have specifically identified a fix.
GenoMax is offline   Reply With Quote
Old 03-29-2018, 11:25 PM   #609
tabwalsh
Junior Member
 
Location: Vienna

Join Date: Mar 2018
Posts: 2
Default

Thanks, I've submitted a ticket here: https://sourceforge.net/p/bbmap/tickets/7/
tabwalsh is offline   Reply With Quote
Old 04-04-2018, 08:29 AM   #610
bwubb
Member
 
Location: Philadelphia

Join Date: Jan 2012
Posts: 55
Default

Greetings,

First time bbmap user and Im having some odd compatibility issues with certain tools using a bbmap generated bam. I do not think I did anything particularly odd in my bam workflow:

Left trim and right trim using bbduk | bbmerge | bbmap | samtools addreplacerg

From this I was able to successfully make variant calls and generate some metrics using GATK/Picard. However bedtools was reporting zero coverage at all intervals and I am unable to view reads restricting to regions using samtools.

Code:
[[email protected] GRCh37]$ samtools view 4721.ready.bam 1
[main_samview] region "1" specifies an unknown reference name. Continue anyway.


[[email protected] GRCh37]$ samtools view -H 4721.ready.bam | head -5
@HD     VN:1.4  SO:coordinate
@SQ     SN:1 dna:chromosome chromosome:GRCh37:1:1:249250621:1   LN:249250621
@SQ     SN:2 dna:chromosome chromosome:GRCh37:2:1:243199373:1   LN:243199373
@SQ     SN:3 dna:chromosome chromosome:GRCh37:3:1:198022430:1   LN:198022430
@SQ     SN:4 dna:chromosome chromosome:GRCh37:4:1:191154276:1   LN:191154276
Are they SN incorrect? Nothing looks really unusual. I guess the first line of my fastq does look like

Code:
>1 dna:chromosome chromosome:GRCh37:1:1:249250621:1
If this is the issue can I "correct" the bbmap index or should I just modify the fasta? Thanks.

-bwubb
bwubb is offline   Reply With Quote
Old 04-06-2018, 03:23 PM   #611
luc
Senior Member
 
Location: US

Join Date: Dec 2010
Posts: 337
Default

Convenient way to avoid re-indexing references?

When working with big genomes one would really like to avoid re-indexing reference sequences. BBmap is always looking for the reference indices in the working directory and not in the directory where the reference resides. Is there any way to point BBMap to the index? When using "path=" it seems to automatically re-index i this location.
Thanks in advance.
luc is offline   Reply With Quote
Old 04-07-2018, 06:17 AM   #612
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,783
Default

Quote:
Originally Posted by bwubb View Post
Greetings,

First time bbmap user and Im having some odd compatibility issues with certain tools using a bbmap generated bam. I do not think I did anything particularly odd in my bam workflow:

Left trim and right trim using bbduk | bbmerge | bbmap | samtools addreplacerg

From this I was able to successfully make variant calls and generate some metrics using GATK/Picard. However bedtools was reporting zero coverage at all intervals and I am unable to view reads restricting to regions using samtools.

Code:
[[email protected] GRCh37]$ samtools view 4721.ready.bam 1
[main_samview] region "1" specifies an unknown reference name. Continue anyway.


[[email protected] GRCh37]$ samtools view -H 4721.ready.bam | head -5
@HD     VN:1.4  SO:coordinate
@SQ     SN:1 dna:chromosome chromosome:GRCh37:1:1:249250621:1   LN:249250621
@SQ     SN:2 dna:chromosome chromosome:GRCh37:2:1:243199373:1   LN:243199373
@SQ     SN:3 dna:chromosome chromosome:GRCh37:3:1:198022430:1   LN:198022430
@SQ     SN:4 dna:chromosome chromosome:GRCh37:4:1:191154276:1   LN:191154276
Are they SN incorrect? Nothing looks really unusual. I guess the first line of my fastq does look like

Code:
>1 dna:chromosome chromosome:GRCh37:1:1:249250621:1
If this is the issue can I "correct" the bbmap index or should I just modify the fasta? Thanks.

-bwubb
BBMap is one of the aligners that uses the full header present in your fasta file when creating the index and passes it along to alignment file. If there are spaces in the header name they are written to alignment. Some downstream programs have a problem with this.

You can use the option "trd=t" to truncate the fasta header names after the first space in the name. There is an option for reformat.sh that can do this after the fact for aligned data. I can look this up later if you don't find it.
GenoMax is offline   Reply With Quote
Old 04-07-2018, 06:20 AM   #613
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,783
Default

Quote:
Originally Posted by luc View Post
Convenient way to avoid re-indexing references?

When working with big genomes one would really like to avoid re-indexing reference sequences. BBmap is always looking for the reference indices in the working directory and not in the directory where the reference resides. Is there any way to point BBMap to the index? When using "path=" it seems to automatically re-index i this location.
Thanks in advance.
Index the genome independently by doing
Code:
bbmap.sh ref=your_genome.fa
This will produce the index in the local directory. There will be a top level "ref" directory with several things in it. Don't worry about the exact organization/names of files in it. This is the pre-made index.

When you want to do the alignments instead of "ref=" you pass "path=/path_to_dir_containing_ref_dir". That should use the pre-made index.

Last edited by GenoMax; 04-07-2018 at 09:37 AM.
GenoMax is offline   Reply With Quote
Old 04-07-2018, 10:17 AM   #614
kevin199011
Junior Member
 
Location: Austin

Join Date: May 2014
Posts: 4
Default

Hi Brian, can BBduk deal with trimming adapter sequence only but retain other bases? Forexample:
SEQUENCEADAPTERSEQUENCE; trimming adapter leaves: SEQUENCESEQUENCE?

Thank you!

Kevin
kevin199011 is offline   Reply With Quote
Old 04-09-2018, 09:16 AM   #615
Coyk
Member
 
Location: West Virginia

Join Date: Jun 2011
Posts: 14
Default basic java windows error help

Hello,
I've been trying to get any basic command to work on my windows 7 machine. I've installed bbmap in the C:\Users\me directory. My .fastq.gz files are in the bbmap folder. When I type the command
C:\Users\me>java -ea -Xmx1g -cp \bbmap\current\jgi.BBDukF in1=inputR1.fastq.gz in2=inputR2.fastq.gz out1=test1.fastq out2=test2.fastq ktrim=r mink=8 k=15

I get "could not find or load main class in1=inputR1.fastq.gz"

I've messed with the type of slash, putting the infiles in the current directory, and taking away parts of the directory names. I even truncated the argument after -cp to just jgi.BBDukF, which I know is wrong, and it still wants the main class of in1. I don't know a thing about java and when I read instructions about paths, classpaths, path environments, etc, I get confused quickly. When I type "echo %CLASSPATH% into the command prompt, I get CLASSPATH back so that's not set. I tried using Git Bash but when I typed "bbduk.sh" I got "command not found" so I'm stuck, but I would really like to use this tool on my windows machine. Anyone willing to give me very basic instructions on the exact commands that get bbduk to run on a windows machine?
Coyk is offline   Reply With Quote
Old 04-09-2018, 10:00 AM   #616
JenBarb
Member
 
Location: Bethesda, MD

Join Date: Oct 2010
Posts: 47
Default filterbyname.sh is not pulling out reads

Hello,
Any suggestions as to why filterbyname is not pulling out reads that I know are there.

Help??

Jen
JenBarb is offline   Reply With Quote
Old 04-09-2018, 11:41 AM   #617
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,783
Default

Quote:
Originally Posted by JenBarb View Post
Hello,
Any suggestions as to why filterbyname is not pulling out reads that I know are there.

Help??

Jen
Are you including the entire fastq/fasta header in your retrieval command?
GenoMax is offline   Reply With Quote
Old 04-09-2018, 11:54 AM   #618
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,783
Default

Quote:
Originally Posted by Coyk View Post
Hello,
I've been trying to get any basic command to work on my windows 7 machine. I've installed bbmap in the C:\Users\me directory. My .fastq.gz files are in the bbmap folder. When I type the command
C:\Users\me>java -ea -Xmx1g -cp \bbmap\current\jgi.BBDukF in1=inputR1.fastq.gz in2=inputR2.fastq.gz out1=test1.fastq out2=test2.fastq ktrim=r mink=8 k=15

I get "could not find or load main class in1=inputR1.fastq.gz"

I've messed with the type of slash, putting the infiles in the current directory, and taking away parts of the directory names. I even truncated the argument after -cp to just jgi.BBDukF, which I know is wrong, and it still wants the main class of in1. I don't know a thing about java and when I read instructions about paths, classpaths, path environments, etc, I get confused quickly. When I type "echo %CLASSPATH% into the command prompt, I get CLASSPATH back so that's not set. I tried using Git Bash but when I typed "bbduk.sh" I got "command not found" so I'm stuck, but I would really like to use this tool on my windows machine. Anyone willing to give me very basic instructions on the exact commands that get bbduk to run on a windows machine?
Try the following. Make sure you change /path/to/ to a real value on your computer.
Code:
java -ea -Xmx200m -cp /path/to/bbmap/current/ jgi.BBDukF in=reads.fq out=processed.fq
GenoMax is offline   Reply With Quote
Old 04-11-2018, 07:40 PM   #619
susanklein
Senior Member
 
Location: oceania

Join Date: Feb 2014
Posts: 115
Default bbmap SAM output question

Hi,

I don't quite understand the SAM output of bbmap.sh

here is my line:

bbmap.sh ref=ref.fasta ambig=all nodisk=t nzo=t mappedonly=t outu=unmapped.fasta scafstats=stats.map in=reads.fastq out=mapped.sam

I usually just use the stats file to get read counts per feature using a ref with features as a multifasta, but I need to annotate from a .gff3 file using a draft genome ref and the only way I can see to do that is from the sam (to get the mapping coordinate. However, the sam seems to have low quality mappings too (MAPQ=3)..

What is the threshold actually used for mapping quality? I can see that in the options, and how can I get the sam to only contain the mapped reads above that threshold?

Thanks,

S.
susanklein is offline   Reply With Quote
Old 04-12-2018, 03:38 AM   #620
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,783
Default

Since you are using ambig=all you are going to get multi-mappers with their MAPQ being set to a small value (3). You could post filter your sam file or set ambig=toss (or ambig=random/best to get only one alignment).
GenoMax is offline   Reply With Quote
Reply

Tags
bbmap, metagenomics, rna-seq aligners, short read alignment

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 01:36 AM.


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