[Carol_B]

You can filter proper pairs only with the "po" flag; unpaired reads will then be classified as unmapped. You can use the 'outm' stream to just get mapped reads. You can also filter with samtools after the fact.

From these, you could filter only the pairs where at least one read is uniquely aligned (mapq>3 or some higher threshold for greater confidence), but this can't be done in BBMap currently; you could, however, do this with a custom script or with samtools.

Hi Brian,

I'm also struggling with the selection of unique alignments. I had previously used the option ambiguous=toss, but for some pair-end reads I only get the forward to be considered aligned. I saw your suggestions of removing the "ambiguos=toss" flag and later filter the alignment based on their mapq value. Nevertheless, I am curious on how to select a correct threshold for my filtering. I am using local=t alignment, is the maximum mapq score different when you use local as compared to global alignment? why do you suggest 4 as a good threshold?

Thanks in advance !!!

Carol

[GenoMax]

Hi,

I would like to calculate the average coverage of some genomes. I used:

bbmap.sh in1=reads_R1.fq in2=reads_R2.fq ref=scaffolds.fasta nodisk covstats=stats.txt covhist=histogram.txt

I get nice output files with coverage etc. for each scaffold, but the average for the whole genome is shown only in the terminal. How could I output it to a file, and preferably in a way that I can get tens of genomes coverages summarized in the same file?

Thanks!

Capture the standard out/error streams to a file (would depend on your shell but examples for bash are here) and you will have that information in a file.

[sebl]

coverage stats

Hi,

I would like to calculate the average coverage of some genomes. I used:

bbmap.sh in1=reads_R1.fq in2=reads_R2.fq ref=scaffolds.fasta nodisk covstats=stats.txt covhist=histogram.txt

I get nice output files with coverage etc. for each scaffold, but the average for the whole genome is shown only in the terminal. How could I output it to a file, and preferably in a way that I can get tens of genomes coverages summarized in the same file?

Thanks!

[yaanlpc]

hey, Brian, I am not able to enable the multiple hits feature with bbmap. I actually want to study the indels. I constructed a ref with 3 templates, one original and two with ins and del for 2 and 3 bp, respectively. however, the output indicates that read can only have one hit (the original), even I set maxindel=10 maxsites=10.

my command:

bbmap.sh ref=../input/demo.genome.fa in=../input/query.fastq outm=test.sam threads=1 maxindel=10 maxsites=10

what's the way out?

ref seq:
>chr1_ins
atcggatCCAAAATCTAGCCGGTCCTCTCACAATAGCAAATTGGGCTTGCAAAATggCCCAAATGCTCTCTCCACATCTTTTCTAGCCGCCGCATGAGCATTGTGGAAGTttgac(insertion gg from >chr1)
>chr1
atcggatCCAAAATCTAGCCGGTCCTCTCACAATAGCAAATTGGGCTTGCAAAATCCCAAATGCTCTCTCCACATCTTTTCTAGCCGCCGCATGAGCATTGTGGAAGTttgac
>chr2_del
atcggatCCAAAATCTAGCCGGTCCTCTCACAATAGCAAATTGGGCTTGCAAAATAAATGCTCTCTCCACATCTTTTCTAGCCGCCGCATGAGCATTGTGGAAGTttgac(deletion ccc from >chr1)


input:
@HWI-ST822:1690U49ACXX:8:1101:2952:2236 1:N:0:CAGATC
CCAAAATCTAGCCGGTCCTCTCACAATAGCAAATTGGGCTTGCAAAATCCCAAATGCTCTCTCCACATCTTTTCTAGCCGCCGCATGAGCATTGTGGAAGT
+
@C@FFDFFHHFFDIJDDEHHGDEHIIIBC>FEIADGAHHGHADHEDGIIHGDEGGEGGGGGIGIIEHEECD@CD>?CCED@@;B;9<>:@@CC:@C84:>>

outputonly one hit)
HWI-ST822:1690U49ACXX:8:1101:2952:2236 1:N:0:CAGATC 0 chr1 8 38 101= * 0 0 CCAAAATCTAGCCGGTCCTCTCACAATAGCAAATTGGGCTTGCAAAATCCCAAATGCTCTCTCCACATCTTTTCTAGC
CGCCGCATGAGCATTGTGGAAGT @C@FFDFFHHFFDIJDDEHHGDEHIIIBC>FEIADGAHHGHADHEDGIIHGDEGGEGGGGGIGIIEHEECD@CD>?CCED@@;B;9<>:@@CC:@C84:>> NM:i:0 AM:i:38

[Shini Sunagawa]

Hi Brian,

This sounds really great!

Just FYI, I ran another large clustering job in less then 5 h, which required >1 month (!) using alternative tools, with similar output. As I said, it is only that dedupe does not provide the containment information... .

Best wishes,
Shini

[Brian Bushnell]

Hi Shini,

I will certainly do this if I get a chance. Right now I'm working on 3 new pieces of software, which have priority... but I have more time for it than I did in January, so hopefully I'll be able to get to it sometime soon.

-Brian

[Shini Sunagawa]

Hi Brian,

Since I saw increased activity lately again, I was wondering if you might have thought about the issue we discussed back in January (~post #300). It was about dedupe not writing out exact matched and contained sequence identifiers.

As mentioned before, solving this would make this tool very competitive to existing ones, due to the immense speed-up.

Thanks for your consideration!

Best wishes,
Shini

[GenoMax]

@moistplus: If you were to use bbmap.sh to do the alignments then you would get that information in the alignment report along with the bam file (as long as you have samtools available in $PATH).

[GenoMax]

@lankage: You don't have to align to the short amplicon regions. You could align to the genome (and find out if you have any non-specific amplification along the way).

[Brian Bushnell]

I agree with Susan. BBMap is a global aligner, and not really designed to map reads to substantially shorter reference sequences. But you could try with the flags "minid=0 local", which might work. Note that "semiperfectmode" will not allow a single mismatch or indel, so it's really only useful in special situations; "local" is more appropriate in this situation.

[susanklein]

Hi,

couldn't you just do it the other way around, Have the pacbio as ref and may your short refs to it?

Although I don't understand why you refs are so short.


S.

[lankage]

Hi, Brian

We recently increased our PacBio amplicon size from ~1100 to 3kb. With the smaller amplicon size we were able to map reads to our allele reference sequence library of non-full length allele sequences using "semiperfectmode" to allow for soft-clipping. Im now looking to map ~3kb read sequences obtained from gDNA sequencing to exon reference sequences of ~270 bases a piece and not able to tune the settings to get any mapping results. Is there a way to tune mapPacBio.sh to get hits for regions within long reads to short exon sequences that perfectly match?

[Brian Bushnell]

"ambiguous=best" is a bit misleading, but it means the genomically first location with a maxmimum score will be used. "ambiguous=all" will report all locations within the ambiguity threshold of the first. This does not mean they need exactly the same score; it means that they are very close, so much so that none can be confidently determined to be the correct mapping location. Normally they're identical, but if for example one mapping had a single 1bp deletion and another mapping had two 1bp substitutions, the scores would be different, but would be close enough to be both reported. But if there was a third potential mapping with, say, 5 substitutions, that would be excluded. This can be controlled with the "secondarysitescoreratio" flag; if you set it to 1.0, only mappings with identical scores to the best score will be reported.

[darthsequencer]

Hi I have a couple questions on the terminology used for retaining ambiguous sites using bbmap.

If "ambiguous=best" this means that if there are a bunch of reads all the with the sam score only the first match will be retained? Or does it mean that of all the reads mapping above a score cutoff the first one will be picked?

Along the same lines - for "ambiguous=all" does this mean that if say 5 locations all share the same highest score that they will be reported or does it mean that all locations above the score cutoff will be retained?

[Thias]

Thanks a lot for your answers! Your exemplified replies were really helpful for some more insight.