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 09-22-2014, 12:14 PM   #81
kbseah
Junior Member
 
Location: Germany

Join Date: Aug 2013
Posts: 6
Default

Thanks for the quick reply, Brian!
kbseah is offline   Reply With Quote
Old 09-22-2014, 03:55 PM   #82
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Quote:
Originally Posted by kbseah View Post
Thanks for the quick reply, Brian!
You're welcome!
Quote:
Originally Posted by HGV View Post
Hi Brian
I was looking at the hitstats files and I realized that the %unambiguousReads
can add up to more than 100%, and similarly the unambiguousReads count can be higher than the total input reads...
Fixed now, as of v33.46.
Brian Bushnell is offline   Reply With Quote
Old 09-24-2014, 09:32 AM   #83
andrej-gnip
Junior Member
 
Location: Copenhagen

Join Date: Aug 2014
Posts: 4
Default

Hi,

I'm developing a tool for analysis of sequence reads from viral genetic material, and mapping to reference viral genomes is part of the process. First version used bowtie2, but now I'm trying to make a new version with better user flexibility, more options etc. And I'd like to use bbmap, as it seems to be the best for this purpose. The user manual to bbmap seems more straightforward, and also people I talked to who used both mappers told me bbmap is generally better.

Now, I did some tests, e.g. I ran bbmap on made up sequences to see how it would perform. I created a fastq file with one read, and two fasta files as references. One fasta file had only one sequence similar to the read, while the other fasta file had 4 such sequences (one of which was the same as the one in the other file, and this one was most similar). Naturally, the read was always mapped to the sequence with highest similarity. The mapping quality was the same in the 2 cases. According to that, I can say that mapping quality is completely independent from other sequences in the reference, it only depends on the read and the particular reference sequence to which it was mapped. I'm not sure, though, so I wanted to ask whether this assumption is actually true. Thanks a lot.
andrej-gnip is offline   Reply With Quote
Old 09-24-2014, 09:59 AM   #84
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

The mapping quality is dependent on the other reference sequences, but only if they are within some threshold of similarity (roughly 6 edits) to the best site. If you copy the same reference sequence twice in the fasta file, the read will map ambiguously and get a score of 3 or less. If you add one or two edits to one copy, the read will map to the unedited one but get a reduced score. But if they differ by, say, 10 edits, then the best mapping location will not get any score penalty. The penalty is also influenced by the number of alternative sites; for example, if there are 5 sites that are each 2 edits worse than the best site, that will give a greater score penalty than if there is only one alternative site.
Brian Bushnell is offline   Reply With Quote
Old 10-06-2014, 07:47 AM   #85
kbseah
Junior Member
 
Location: Germany

Join Date: Aug 2013
Posts: 6
Default pileup.sh inconsistent with samtools pileup

Hi Brian,

I tried using the "fastaorf" function in pileup.sh, to look at the read depth for a bunch of ORFs predicted with Prodigal. The input is in Prodigal's output format as specified. However, I get per-orf coverage results (the depthSum field) that are inconsistent with the output from samtools mpileup.

Briefly: I produced a pileup file with samtools mpileup and for each orf, simply summed the read depth (4th column of the pileup file) for each position that falls within the orf.
I double-checked this by converting the Prodigal output to a BED feature table, and used bedtools multicov and the original BAM file to produce a per-feature read depth. This gives per-feature read depths which are not identical to what I got by summing the depths but roughly a multiple (i.e. plotting the depths per orf from both methods against each other gives an approximately linear relationship).
The output from pileup.sh, on the other hand, doesn't give anything close to a linear relationship.

Is there something different in how pileup.sh calculates the per-orf coverage?

Thanks a lot,
Brandon
kbseah is offline   Reply With Quote
Old 10-06-2014, 09:11 AM   #86
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Brandon,

I did introduce a bug recently when I added support for tracking only read start positions rather than total coverage, which manifested in some situations. It's fixed now and I just uploaded the fixed version (33.57). Would you mind downloading that and confirming whether it works correctly?

Thanks!
Brian Bushnell is offline   Reply With Quote
Old 10-07-2014, 12:10 AM   #87
kbseah
Junior Member
 
Location: Germany

Join Date: Aug 2013
Posts: 6
Default

Hi Brian,

Thanks for your reply. Unfortunately the output seems to be the same. Should I send you the output from pileup.sh vs. the output from bedtools so that you can see what I mean?

Best,
Brandon

Quote:
Originally Posted by Brian Bushnell View Post
Brandon,

I did introduce a bug recently when I added support for tracking only read start positions rather than total coverage, which manifested in some situations. It's fixed now and I just uploaded the fixed version (33.57). Would you mind downloading that and confirming whether it works correctly?

Thanks!
kbseah is offline   Reply With Quote
Old 10-07-2014, 08:36 AM   #88
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Quote:
Originally Posted by kbseah View Post
Hi Brian,

Thanks for your reply. Unfortunately the output seems to be the same. Should I send you the output from pileup.sh vs. the output from bedtools so that you can see what I mean?

Best,
Brandon
Yes, please do, as well as the command line and stdout/stderr messages.
Brian Bushnell is offline   Reply With Quote
Old 10-20-2014, 05:38 AM   #89
holmrenser
Junior Member
 
Location: Netherlands

Join Date: Jul 2013
Posts: 9
Default

Hi Brian,

I would be interested in knowing when you intend to publish BBmap in a paper, can you enlighten us?

Cheers
holmrenser is offline   Reply With Quote
Old 10-20-2014, 10:06 AM   #90
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

I hope to have it submitted by the end of this year.

-Brian
Brian Bushnell is offline   Reply With Quote
Old 10-22-2014, 10:24 AM   #91
dbrami
Member
 
Location: San Diego

Join Date: Sep 2008
Posts: 14
Default

Hi Brian,
Thank you for previous answer through email which I am posting:
"Daniel,

Unfortunately, there is not currently a master document, though I certainly want to make one. BBMap is documented in /docs/readme.txt while everything else is just documented in its shellscript. But I do have individual threads on SeqAnswers about the most commonly-used tools:

http://seqanswers.com/forums/showthread.php?t=41057

(first post contains links to the other threads)

For metagenomics, I consider bbmap.sh, bbsplit.sh, bbduk.sh, bbnorm.sh, pileup.sh, dedupe.sh, and stats.sh to be the most directly relevant.

As for your specific questions:

Is there any documentation on the mapping output summary? Ie: How do numbers in mapped pairs vs read 1 mapped differ?

To be considered proper pairs, by default, the reads must map to opposite strands, on the same scaffold, with an insert size of under 32kbp. Read 1 will typically have a higher mapping rate than the pairing rate, because sometimes the reads map discordantly, and sometimes only one read maps. This is especially true in a metagenome which may have very short scaffolds.

Can we output the mapping results to .fast.gz directly?

All tools in the BBMap package can input and output gzipped files, fasta or fastq. It detects this based on filename, so a gzipped fastq file should be named "*.fq.gz" or "*.fastq.gz".

Is it possible to bin out mappings where both reads map vs either read maps vs unmapped through command line without bitwise comparison of SAM flags?

Yes. BBMap supports a "po" or "pairedonly" flag; you can for example do this:

bbmap.sh in1=read1.fq in2=read2.fq po=t outm=mapped.fq outu=unmapped.fq

...which will split the reads into proper pairs with both mapped (outm stream) and pairs where either one or both didn't map, or they did not map in the proper orientation (outu stream). With "po=f" (pairedonly=false), outm will accept anything where at least one read mapped, and outu will only get pairs in which neither mapped."
I am trying to find a value in mapping summary to would correspond to bowtie2's '% overall alignment rate'

Also, which BBMap parameters should I tweak if I want to accurately want to reproduce bowtie2 results obtained using "-k 1 --fast-local" ?
dbrami is offline   Reply With Quote
Old 10-22-2014, 10:48 AM   #92
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Daniel,

Quote:
Originally Posted by dbrami View Post
I am trying to find a value in mapping summary to would correspond to bowtie2's '% overall alignment rate'

Also, which BBMap parameters should I tweak if I want to accurately want to reproduce bowtie2 results obtained using "-k 1 --fast-local" ?
For paired reads, the overall alignment rate is ((percent mapped for read 1)+(percent mapped for read 2))/2. Or in other words half of the sum of the two numbers in red below.

Code:
Pairing data:           pct reads       num reads       pct bases          num bases

mated pairs:             99.5206%            3114        99.5206%             934200
bad pairs:                0.3196%              10         0.3196%               3000
insert size avg:          321.31


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

mapped:                  99.9361%            3127        99.9361%             469050
unambiguous:             98.6897%            3088        98.6897%             463200
ambiguous:                1.2464%              39         1.2464%               5850
low-Q discards:           0.0000%               0         0.0000%                  0

perfect best site:        1.9175%              60         1.9175%               9000
semiperfect site:         1.9175%              60         1.9175%               9000
rescued:                  0.0000%               0

Match Rate:                   NA               NA        62.0582%             441661
Error Rate:              96.5462%            3019        37.6163%             267711
Sub Rate:                87.4320%            2734         2.3412%              16662
Del Rate:                42.8846%            1341        34.0933%             242638
Ins Rate:                49.8881%            1560         1.1818%               8411
N Rate:                  49.4084%            1545         0.3254%               2316


Read 2 data:            pct reads       num reads       pct bases          num bases

mapped:                  99.9041%            3126        99.9041%             468900
unambiguous:             98.6897%            3088        98.6897%             463200
ambiguous:                1.2144%              38         1.2144%               5700
low-Q discards:           0.0000%               0         0.0000%                  0

perfect best site:        1.6619%              52         1.6619%               7800
semiperfect site:         1.6619%              52         1.6619%               7800
rescued:                  0.0000%               0

Match Rate:                   NA               NA        61.2099%             441930
Error Rate:              96.7370%            3024        38.4685%             277739
Sub Rate:                88.2278%            2758         2.2755%              16429
Del Rate:                43.9859%            1375        35.0546%             253091
Ins Rate:                49.1683%            1537         1.1384%               8219
N Rate:                  50.7358%            1586         0.3216%               2322
Flags equivalent to bowtie2's "-k 1 --fast-local" would be something like:

"fast local maxindel=10"
Brian Bushnell is offline   Reply With Quote
Old 10-22-2014, 02:38 PM   #93
dbrami
Member
 
Location: San Diego

Join Date: Sep 2008
Posts: 14
Default de-interleave fastq outtput

Can BBMap de-interleave the fastq output for reads 1 and reads 2?
Ie: if I Use "pairedonly=f outu=${B}.noplant.fastq.gz", can BBMap put the unmapped forward reads ${B}_R1.noplant.fastq.gz and the unmapped reverse reads in ${B}_R2.noplant.fastq.gz.
Or should I just de-interleave the reads as a post process?
dbrami is offline   Reply With Quote
Old 10-22-2014, 02:46 PM   #94
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Daniel,

You can set "outu1=" and "outu2=" or "outm1=" and "outm2=" for output in paired files.
Brian Bushnell is offline   Reply With Quote
Old 11-07-2014, 05:46 AM   #95
prashant_singh
Junior Member
 
Location: Farmington CT

Join Date: Nov 2014
Posts: 1
Default Issues with BBMap

Hi
I am a newbie to NGS data analysis- basically a bench scientist trying to do analysis.
I have sequenced human tissue samples on the PacBio RSII and am trying to align to human genome using the bbmap tool.

As mentioned in the forum, I am trying to use the MapPacBio.sh script for the alignment. Previously I was getting an error where it could not find the class align2.MapPacBio . I hard coded the path in the file and then

Here is the new error I am getting:

java -Djava.library.path=/opt/compsci/bbmap/33.73/jni/ -ea -Xmx51798m -cp BBMapPacBio build=1 overwrite=true minratio=0.40 fastareadlen=6000 ambiguous=best minscaf=100 startpad=10000 stoppad=10000 midpad=6000 in=/data/*****/PacBio/****/consensus_isoforms.fasta out=/data/***/pac88 ref=/data/shared/research_pipelines_reference_data/human/RNA/genome.fa
Error: Could not find or load main class build=1

Java version is "1.7.0_51"
OpenJDK Runtime Environment (rhel-2.4.4.1.el6_5-x86_64 u51-b02)
OpenJDK 64-Bit Server VM (build 24.45-b08, mixed mode)


Could you please help?

Thanks
Prashant
prashant_singh is offline   Reply With Quote
Old 11-07-2014, 06:09 AM   #96
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

@prashant_singh. I am hardly an expert on BBTools. Brian will probably log in soon -- he is on USA west coast time as far as I know -- and so will give a better answer. But I happened to start up a mapPacBio run just as your post came in. Comparing your Java to my Java reveals:

Yours ...
Code:
java 
  -Djava.library.path=/opt/compsci/bbmap/33.73/jni/ 
  -ea -Xmx51798m 
  -cp BBMapPacBio 
  build=1 
  overwrite=true 
  minratio=0.40 
  fastareadlen=6000 
  ambiguous=best 
  minscaf=100 
  startpad=10000 
  stoppad=10000 
  midpad=6000 
  in=/data/*****/PacBio/****/consensus_isoforms.fasta 
  out=/data/***/pac88  
  ref=/data/shared/research_pipelines_reference_data/human/RNA/genome.fa
Mine ...
Code:
java 
  -Djava.library.path=/group/bioinfo/apps/apps/BBMap-33.57/jni/ 
  -ea -Xmx126058m 
  -cp /group/bioinfo/apps/apps/BBMap-33.57/current/ 
  align2.BBMapPacBio 
  build=1 
  overwrite=true 
  minratio=0.40 
  fastareadlen=6000 
  ambiguous=best 
  minscaf=100 
  startpad=10000 
  stoppad=10000 
  midpad=6000 
  ref=kmer60-10.fa 
  in=all_spades_2000plus.fa 
  out=LR_vs_scaffolds.sam
The major different is the 'cp' where I have the complete path plus the line after the 'cp' which shows the program to run.

You mentioned:

Quote:
Previously I was getting an error where it could not find the class align2.MapPacBio . I hard coded the path in the file and then
....
Error: Could not find or load main class build=1
I suspect that you mis-coded the path and probably erased the part which states the program to run. I suggest starting anew and if you really need to hard-code the path then be careful in doing so.
westerman is offline   Reply With Quote
Old 11-07-2014, 08:38 AM   #97
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

That's correct, and yes, I'm on West Coast time The command line should probably be:

java -ea -Xmx51798m -cp /opt/compsci/bbmap/33.73/current/ align2.BBMapPacBio build=1 overwrite=true minratio=0.40 fastareadlen=6000 ambiguous=best minscaf=100 startpad=10000 stoppad=10000 midpad=6000 in=/data/*****/PacBio/****/consensus_isoforms.fasta out=/data/***/pac88.sam ref=/data/shared/research_pipelines_reference_data/human/RNA/genome.fa

Note both of the things in bold: "/opt/compsci/bbmap/33.73/current/" is probably correct, but maybe not, depending on how it was installed. I'm guessing that somewhere under "/opt/compsci/bbmap/33.73/" is a subdirectory named "current" - that is what you want to point it to.

For "out=/data/***/pac88.sam", I just added the ".sam" suffix. All of the programs in BBTools are extension-sensitive, so if you don't include an extension, it will just output the default format (which in this case is sam, but it's always best to specify it).
Brian Bushnell is offline   Reply With Quote
Old 11-07-2014, 03:21 PM   #98
sdriscoll
I like code
 
Location: San Diego, CA, USA

Join Date: Sep 2009
Posts: 438
Default

I apologize if this has been covered already but I've just downloaded and run bbmap a few times (seems great so far and I like the control provided via the available command line options) and I have a quick question. When building the index is it possible for bbmap to ignore parts of the reference names past the first blank space? For example I'm using the Ensembl human release and the FASTA reference names look something like this:

>1 dna_sm:chromosome chromosome:GRCh37:1:1:249250621:1 REF

Currently when I run alignments on this FASTA file that entire name shows up in the SAM file while instead I should only see "1". Other aligners (STAR, bowtie, bwa, etc) throw out everything after the first space when indexing leaving just the main reference name. This is particularly important for using aligners as part of a gene expression pipeline using something like eXpress since some transcriptome extraction tools (programs that parse a GTF and extract transcripts from an accompanying genome FASTA file) usually name the transcripts as >TRANSCRIPT_ID GENE_SYMBOL. eXpress will throw out everything after the first space however if the alignments contain the full reference names then no alignments are matched up to transcripts in the output.
__________________
/* Shawn Driscoll, Gene Expression Laboratory, Pfaff
Salk Institute for Biological Studies, La Jolla, CA, USA */
sdriscoll is offline   Reply With Quote
Old 11-07-2014, 03:26 PM   #99
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Yes, I put in an option for that because some programs require it. The flag is trd=t (or trimreaddescriptions=true), which will truncate all names (both read names and reference names) at the first whitespace symbol.
Brian Bushnell is offline   Reply With Quote
Old 11-07-2014, 03:43 PM   #100
sdriscoll
I like code
 
Location: San Diego, CA, USA

Join Date: Sep 2009
Posts: 438
Default

rad. thanks for that.
__________________
/* Shawn Driscoll, Gene Expression Laboratory, Pfaff
Salk Institute for Biological Studies, La Jolla, CA, USA */
sdriscoll 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 06:12 AM.


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