[Brian Bushnell]

I think you found a bug ...
When you look at the given line 625 in the file ./current/shared/TrimRead.java, it evaluates the CIGAR string and asserts that it ends with a match (M or =). This should be the case because the trimReadWithMatchFast() function is only invoked for 100% match or unmapped reads:

char c=sl.cigar.charAt(sl.cigar.length()-1);
assert(c=='M' || c=='=') : c+"; "+sl.cigar+"\n"+sl;

But in your case, you have the hard-trimmed right part, so c is a H and the assertion fails, even though the function is rightfully invoked, since your remaining 22 bases are all matching.

Yes, that is exactly correct; good analysis. I need to add support for H symbols, which I did not expect to encounter in that function.

[Thias]

I think you found a bug ...
When you look at the given line 625 in the file ./current/shared/TrimRead.java, it evaluates the CIGAR string and asserts that it ends with a match (M or =). This should be the case because the trimReadWithMatchFast() function is only invoked for 100% match or unmapped reads:

char c=sl.cigar.charAt(sl.cigar.length()-1);
assert(c=='M' || c=='=') : c+"; "+sl.cigar+"\n"+sl;

But in your case, you have the hard-trimmed right part, so c is a H and the assertion fails, even though the function is rightfully invoked, since your remaining 22 bases are all matching.

[CharlesARoy]

Hi Brian Bushnell

I recently installed BBMap_38.98.tar.gz from SourceForge on a machine running Ubuntu 20.04.4 LTS.

Per the usage guide, I checked my current version of Java with "java -Xmx90m -version" which returned the following:

openjdk version "11.0.15" 2022-04-19
OpenJDK Runtime Environment (build 11.0.15+10-Ubuntu-0ubuntu0.20.04.1)
OpenJDK 64-Bit Server VM (build 11.0.15+10-Ubuntu-0ubuntu0.20.04.1, mixed mode, sharing)

I tried running reformat on a sam file with a command like:

reformat.sh in=in.sam out=out.sam forcetrimleft=10

And got the following error:

java -ea -Xms300m -cp /data/usr/croy/tools/bbmap/current/ jgi.ReformatReads in=in.sam out=out.sam forcetrimleft=10
Executing jgi.ReformatReads [in=in.sam, out=out.sam, forcetrimleft=10]

Input is being processed as unpaired
Waiting on header to be read from a sam file.
Exception in thread "main" java.lang.AssertionError: H; 22M123H
FS10001085:123:BRL95608-3321:1:1101:7680:1000 2195 chr3 50608025 18 22M123H = 50607939 -98 AGTGGTTTTCACTGACAGCGTG F,:FFFFFFFF::FFF:F,F,F NM:i:0 MD:Z:22 MC:Z:93M21D52M AS:i:22 XS:i:0 SA:Z:chr3,50608053,-,17S128M,60,1;
at shared.TrimRead.trimReadWithMatchFast(TrimRead.java:625)
at shared.TrimRead.trimReadWithMatch(TrimRead.java:684)
at shared.TrimRead.trimByAmount(TrimRead.java:279)
at shared.TrimRead.trimToPosition(TrimRead.java:240)
at jgi.ReformatReads.process(ReformatReads.java:919)
at jgi.ReformatReads.main(ReformatReads.java:53)

The sam file was aligned to hg38 using a recent version of BWA MEM (v0.7.17-r1188). The CIGAR string seems to follow the SAM specs as far as I can tell, so I'm not sure where things are going wrong.

Please let me know if you need any more information.

Best,

Charles

[Omri.huji]

Is bbmap missing mappings of bbmapskimmer ?

Hello Brian,

I will be glad to get your help with the next issue.
I mapped a large DB of reads(over 100K) against the human genome, and found different results between bbmap's and bbmapskimmer's mappings.

I used bbmap to map reads against hg19, then I was interested finding all the mappings over the genome- so I used bbmapskimmer.
Except this concept I do not know any main differences between bbmap and bbmapskimmer codes (am I right?).

And still, bbmap mapped 98,425 reads and bbmapskimmer mapped 107,708. While 96,987 reads are shared.
Means bbmap mapped 1,438 reads that wasn't mapped by bbmapskimmer, which did mapped another 10,721 that wasn't mapped by bbmap- some of them with perfect minid score of 1.0.
Is there any possible reason for that?

I attached SAM files, with the relevant mappings for reads with maximum minid score in each algorithem (best minid for 1,438 reads of bbmap is 0.96, and for bbmapskimmer is 1.0).


- I used the next reference genome: hg19_main_mask_ribo_animal_allplant_allfungus.fa.gz

- All mappings were done with next parameters: minid=0.9 ambig=all

[sghignone]

filter reads based on GC% content

Hello Brian and everyone,
I'm wondering whether it could be possible (and with which tool) to perform a read filter based on its GC% content.
For example, working with a metagenome with fungal and bacterial DNA fragments, how can I retain only paired-end reads with a GC > 50% and isolate bacterial reads?
Thanks for the attention
Stefano

[ellalalalala]

Error running BBMap

Hi!

I already used BBMap successfully in the past to map RNA-seq data to a multi-FASTA transcriptome. Now I was trying to map to the human genome (downloaded from human genome browser of NCBI). I used the following command:

bbmap.sh ref=$mapping_ref build="$BBmap_ref_ID" in=$BBduked_reads1 in2=$BBduked_reads2 \
outu=unmapped_"$name1"_"$date".sam outm=mapped_"$name1"_"$date".sam \
maxindel=200k mdtag=true sam=1.4 subfilter=3 pairedonly=true

Unfortunately I´m running into an error:

[...]
Retaining first best site only for ambiguous mappings.
Writing reference.
Executing dna.FastaToChromArrays2 [/scratch/hpc-prf-ptma2/ptma2001/bbmap/GRCh38_genome_16032022.fna, 6, writeinthread=false, genscaffoldinfo=true, retain, waitforwriting=false, gz=true, maxlen=536670912, writechroms=true, minscaf=1, midpad=300, startpad=8000, stoppad=8000, nodisk=false]

Set genScaffoldInfo=true
Writing chunk 1
Writing chunk 2
Exception in thread "main" java.lang.NumberFormatException: For input string: "ӌ�jHN�H��¬���l��,ΚK��/S��ͩm���T
�!�l��X��o�ã’*s$�s^��m�h
J@�
�m��$֧����@�H,���xI%_ Uz�?� S�B�2[*i�_����)i�NEր�ٜ�?Z�i£�� "
at java.lang.NumberFormatException.forInputString(NumberFormatException.java:65)
at java.lang.Integer.parseInt(Integer.java:580)
at java.lang.Integer.parseInt(Integer.java:615)
at dna.Data.setGenome2(Data.java:1016)
at dna.Data.setGenome(Data.java:769)
at align2.BBMap.loadIndex(BBMap.java:316)
at align2.BBMap.main(BBMap.java:32)

Any idea what causes the problem and how I could fix it?

Thanks a lot in advance and have a great day

Ella

[Ahmer]

BBmap Error

I am getting this error how to solve it.

java -ea -Xmx4393m -Xms4393m -cp /home/szweda/bbmap/current/ align2.BBMap build=1 overwrite=true fastareadlen=500 in=PY53_contigs.fasta out=PY53_bbmapped.sam bamscript=bs.sh
Executing align2.BBMap [build=1, overwrite=true, fastareadlen=500, in=PY53_contigs.fasta, out=PY53_bbmapped.sam, bamscript=bs.sh]
Version 38.95

Retaining first best site only for ambiguous mappings.
Found samtools 1.10
Exception in thread "main" java.lang.RuntimeException: Can't find file ref/genome/1/summary.txt
at fileIO.ReadWrite.getRawInputStream(ReadWrite.java:933)
at fileIO.ReadWrite.getInputStream(ReadWrite.java:898)
at fileIO.TextFile.open(TextFile.java:280)
at fileIO.TextFile.<init>(TextFile.java:123)
at dna.Data.setGenome2(Data.java:823)
at dna.Data.setGenome(Data.java:769)
at align2.BBMap.loadIndex(BBMap.java:316)
at align2.BBMap.main(BBMap.java:32)

[gtwa-bio]

Entropy filtering

Hello,
I am working with shotgun metagenomic sequencing data from gut microbiome samples. We're planning to do taxonomic abundance estimates. For data preprocessing we are going to trim and filter sequences for quality and adapter content with bbduk, as well as remove host sequences with bbsplit. We are also considering an additional entropy filtering step with bbduk after our other preprocessing steps. I was hoping I could ask you for some more information about how this entropy filtering process works, and if you might have a recommendation in our case.

• In some of our samples, we have an overrepresentation of G homopolymers. We’re confident that these are technical artifacts from the NextSeq sequencing protocol. I know we can filter these out with an entropy threshold of 0.1. However, I’ve seen in some metagenomic studies they filter out repetitive sequences that are not sequencing artifacts. Would you recommend raising this entropy threshold in our case, and if so to what new value?
• On the more technical side of how bbduk implements entropy filtering…
  • How does the sliding window traverse a read? Is it one base at a time, or does it move by the whole window size?
  • If a read has a region of low complexity sequences at the beginning/end, are only these sections filtered or is the entire read removed?
  • How might a read with an internal region of low complexity be treated?

[Franz]

How is coverage calculated in bbmap?

Hi Brian and everyone else,

I am using bbmap for my metagenomes. Currently, I am comparing the coverage of different genes across metagenomes. BBmap has the flag -rpkm, which calculates fold coverage, RPKM and FPKM. But I couldn´t find the respective formulas, how exactly it is calculated. When I tried to re-calculate these parameters, I came to different results.
Thank you very much in advance! I hope this thread is at the right position in SEQanswers.

Best,
Franz

[turnersd]

Those genotype calls are based on the read data, which are summarized in the INFO fields of the VCF. If you change the genotypes, those field data will be inconsistent.

It would be useful if you explained more clearly what you're trying to accomplish.

Yes, completely correct about GT not matching up with INFO fields. I'll probably end up removing all the info fields altogether because GT is really the only thing I'm interested in.

I've used a published simulation framework to simulate pedigrees having genotype data, and what I'm trying to do is specify a rate at which I want to mutate genotypes, i.e., to mimic allele drop-out or drop-in, each at specified rates. Mimicking a het to hom-ref would be easy if variant alleles were all I cared about (ie just deleting random rows from the VCF), but I _do_ care about genotypes at invariant sites.

My plan was to write some custom bash/python/bcftoolsy things to get out genotypes and mutate at a specified rate then stick them back into the VCF. I was just wondering if there was some existing tool to do something like this in bbmap (or elsewhere). I can't seem to find anything.

[HESmith]

Not exactly. I literally want to replace the 0/1s or 0/0s or 1/1s in the VCF, not generate reads.

Those genotype calls are based on the read data, which are summarized in the INFO fields of the VCF. If you change the genotypes, those field data will be inconsistent.

It would be useful if you explained more clearly what you're trying to accomplish.

[turnersd]

You can do it in two steps: 1) use mutate.sh on the reference genome; 2) use randomreads.sh of the mutated reference to generate the data (I assume that's what you mean by sample).

Not exactly. I literally want to replace the 0/1s or 0/0s or 1/1s in the VCF, not generate reads.

[HESmith]

You can do it in two steps: 1) use mutate.sh on the reference genome; 2) use randomreads.sh of the mutated reference to generate the data (I assume that's what you mean by sample).

[turnersd]

Insert mutations in a VCF

I see mutate.sh takes a vcf and inserts mutations in a reference genome. I need something a little different. I want to specify a mutation rate and potentially a sample, and insert mutations at that specified rate in the specified sample. Anything in bbmap that might help with this?

[reliscu]

All &quot;N&quot; reads in paired reads not being filtered if other read is &quot;good&quot;

Version 38.94

This thread details the "bug":
https://www.biostars.org/p/9493307/#9493339