[darthsequencer]

Hi Brian - Is there a bbtool that will keep only mapped reads and if the mate doesn't map output the unmapped also from sam/bam files?

Thanks

[susanklein]

Ok I thought this was solved, but see below, still getting sam file alignments with MAPQ=3, even though they actually look like good alignments. Is '3' given because you can't calculate a perfect hit MAPQ? i.e. MAPQ=-10*log10(zero)

@PG ID:BBMap PN:BBMap VN:37.93 CL:java -Djava.library.path=/bin/bbmap/jni/ -ea -Xmx43986m align2.BBMap build=1 overwrite=true fastareadlen=500 ref=genome/genome.fasta ambig=best vslow perfectmode maxsites=10000 nzo=t mappedonly=t outu=test/unmapped/testKT45.unmapped.fasta scafstats=test/mapped/testKT45.map in=testKT45.fastq out=test/mapped/testKT45.sam -idfilter=0.9
700666F:325:CCC50ANXX:6:2211:7729:2233 1:N:0:GTCGAG 0 ctg3849_RHIMB__RM170330.3849 Backbone_3849|arrow|pilon 111580 3 31= * 0 0 GCATTGGTGGTTCAGTGGTAGAATTCTCGCC CCCCCGGGGGGGGGGGGGGGGGGGGGGGGGG XT:A:R NM:i:0 AM:i:3
700666F:325:CCC50ANXX:6:2211:7914:2047 1:N:0:GTAGAG 16 ctg18961_RHIMB__RM170330.18961 Backbone_18961|arrow|pilon 20097 3 31= * 0 0 TCACAATGATAGGAAGAGCCGACATCGAAGG GGGGGGGGGGGGGGGGGGGGGGGGGGCCCCC XT:A:R NM:i:0 AM:i:3
700666F:325:CCC50ANXX:6:2211:7905:2066 1:N:0:GTCGAG 16 ctg18961_RHIMB__RM170330.18961 Backbone_18961|arrow|pilon 34735 3 31= * 0 0 GGTAGGCGCAGAAAGTACCATCGAAAGTTGA GGGGGFGGGGGGGGGGGGGGGGGGGGCCCCC XT:A:R NM:i:0 AM:i:3
700666F:325:CCC50ANXX:6:2211:7970:2099 1:N:0:GTCGAG 0 ctg16585_RHIMB__RM170330.16585 Backbone_16585|arrow|pilon 33974 3 23= * 0 0 GGGAATACCAGGTGCTGTAGGCT CCCCCGEGGGGGGGGGGGGGGGG XT:A:R NM:i:0 AM:i:3
700666F:325:CCC50ANXX:6:2211:7937:2145 1:N:0:GTCGAG 16 ctg13544_RHIMB__RM170330.13544 Backbone_13544|arrow|pilon 7627 3 39= * 0 0 GGTGAGAGCGCCGAATCCTAACCACTAGACCACCAGGGA GGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGCCCCC XT:A:R NM:i:0 AM:i:3

[susanklein]

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).

ambig=best worked, Thanks.

S.

[GenoMax]

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).

[susanklein]

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.

[GenoMax]

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.

java -ea -Xmx200m -cp /path/to/bbmap/current/ jgi.BBDukF in=reads.fq out=processed.fq

[GenoMax]

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?

[JenBarb]

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

[Coyk]

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?

[kevin199011]

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

Thank you!

Kevin

[GenoMax]

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

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.

[GenoMax]

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.

[bwubb@node063 GRCh37]$ samtools view 4721.ready.bam 1
[main_samview] region "1" specifies an unknown reference name. Continue anyway.


[bwubb@node063 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

>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.

[luc]

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.

[bwubb]

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.

[bwubb@node063 GRCh37]$ samtools view 4721.ready.bam 1
[main_samview] region "1" specifies an unknown reference name. Continue anyway.


[bwubb@node063 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

>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

[tabwalsh]

Thanks, I've submitted a ticket here: https://sourceforge.net/p/bbmap/tickets/7/