[Brian Bushnell]

No, BBMap does not output chimerically-mapped reads to a separate file, though it can be used to separate properly-paired reads from improper (possibly chimeric) pairs.

[Dario1984]

But can it output chimerically mapped reads to a separate output file, like STAR can? There's no mention of chimeric reads in the user guide, so I'm not sure if it's even suitable for that case.

[Brian Bushnell]

Yes, it's random.

[GenoMax]

With the new version, I don't have anymore negative coverage!

Another question:
Using this command :



in1 and in2 are Forward reads and Reverse reads right ? So it keeps the mated pair at the end ?

It should do that.

[Brian Bushnell]

Indeed, I see 3 entries with negative coverage, which should not happen. However, v32.15 was released about 2.5 years ago and there have been thousands of changes since then, including hundreds of bug fixes. Could you try the latest version and see if that fixes the problem?

[Brian Bushnell]

I've ran pileup for contig coverage estimation after assembly.

The output is :

ID Avg_fold Length Ref_GC Base_Coverage Read_GC

1) The coverage is the Avg_fold right ?

Yes, that's correct.

2) If yes, I have some negative or 0 coverage ... How can it be possible ?

Zero coverage is certainly possible, since mapping and assembly results can differ. Negative coverage should be impossible. Can you send me the output file, and post your exact command line?

Thanks,
Brian

[Brian Bushnell]

Thanks Brian,

I took the de novo approach to try to eliminate contaminating sequences. We're all working on HIV, and it's highly mutated: it's hard to remove contaminants by DNA sequence alone. However, since I know the approximate size of my PCR amplicons (it's different in each case because of DNA deletions), I know what size my contig should be after de novo assembly. By using stringent assembly parameters, that require ~30 bp overlap with maybe 1 mismatch, I can force contaminants to be assembled into a separate contig. I can then extract consensus sequences from each of these contigs, and filter them so that only contigs greater than ~150 bp are used in a subsequent map to reference, which should remove any minor contaminants that were present. I can then extract the consensus from this alignment, which should represent the original PCR amplicon - any contaminants that might have made it into my contig should be lost as they are outnumber by the 100s.

Does this workflow make sense for my application? I'm working on a desktop Mac, so I have limited options. I've been told that Spades might be a better assembler for me, but I think I'd need to purchase a server. With the little coding experience I have, I'm a bit nervous to invest the money, lest we never get the thing to work.

Are you available for hire as a consultant?

Thanks,
Jake

Hi Jake,

It might be helpful in this case if you could generate a kmer frequency histogram to see whether the contaminant and non-contaminant sequence is easily separable by depth alone. If so, there are a couple of easy ways to remove it. You can generate a kmer-frequency histogram with kmercountexact or BBNorm; just attach the text file to this thread. Normally I look at it in a log-log plot.

What assembler are you currently using, by the way? I've had poor results with Spades on viruses, and better results with Tadpole. But this was raw viral sequence and amplicon sequence may give different results.

As for consultancy, I've sent you a pm.

-Brian

[JVGen]

Hi Jake,

By default, Dedupe will remove all but one copy of reads that are fully contained within other reads. This may not be a good idea prior to assembly of merged reads (since many will be contained by others due to the variable insert size); just using normalization alone with a target of 100x may work better. You can set Dedupe to only remove reads that are identical to others, though, by adding the flag "ac=f".

If you suspect your data is contaminated, I suggest you try to decontaminate it prior to normalization/assembly. If you are working on non-HIV and it is contaminated with HIV, that's fairly straightforward; you can use BBSplit or Seal (if you have both references) or just BBMap (if you only have the HIV reference but nothing for your amplicon) to separate the reads by mapping. Alternately, if you are confident that all of the low-coverage stuff is contamination, you can assemble the low-coverage portion of your data using Tadpole, which allows coverage thresholds (e.g. "mincr=1 maxcr=15" will make it assemble only those contigs with coverage between 1x and 15x). Then use the resulting assembly to filter all of the reads via mapping, keeping only the unmapped reads.

Incidentally, Tadpole tends to do a good job assembling very-high-coverage short sequences.

Thanks Brian,

I took the de novo approach to try to eliminate contaminating sequences. We're all working on HIV, and it's highly mutated: it's hard to remove contaminants by DNA sequence alone. However, since I know the approximate size of my PCR amplicons (it's different in each case because of DNA deletions), I know what size my contig should be after de novo assembly. By using stringent assembly parameters, that require ~30 bp overlap with maybe 1 mismatch, I can force contaminants to be assembled into a separate contig. I can then extract consensus sequences from each of these contigs, and filter them so that only contigs greater than ~150 bp are used in a subsequent map to reference, which should remove any minor contaminants that were present. I can then extract the consensus from this alignment, which should represent the original PCR amplicon - any contaminants that might have made it into my contig should be lost as they are outnumber by the 100s.

Does this workflow make sense for my application? I'm working on a desktop Mac, so I have limited options. I've been told that Spades might be a better assembler for me, but I think I'd need to purchase a server. With the little coding experience I have, I'm a bit nervous to invest the money, lest we never get the thing to work.

Are you available for hire as a consultant?

Thanks,
Jake

[Brian Bushnell]

Hi Jake,

By default, Dedupe will remove all but one copy of reads that are fully contained within other reads. This may not be a good idea prior to assembly of merged reads (since many will be contained by others due to the variable insert size); just using normalization alone with a target of 100x may work better. You can set Dedupe to only remove reads that are identical to others, though, by adding the flag "ac=f".

If you suspect your data is contaminated, I suggest you try to decontaminate it prior to normalization/assembly. If you are working on non-HIV and it is contaminated with HIV, that's fairly straightforward; you can use BBSplit or Seal (if you have both references) or just BBMap (if you only have the HIV reference but nothing for your amplicon) to separate the reads by mapping. Alternately, if you are confident that all of the low-coverage stuff is contamination, you can assemble the low-coverage portion of your data using Tadpole, which allows coverage thresholds (e.g. "mincr=1 maxcr=15" will make it assemble only those contigs with coverage between 1x and 15x). Then use the resulting assembly to filter all of the reads via mapping, keeping only the unmapped reads.

Incidentally, Tadpole tends to do a good job assembling very-high-coverage short sequences.

[JVGen]

As long as the file is unsorted or name-sorted, that command will work. Spliced reads shouldn't cause a problem, unless the splice junction is in the unsequenced area between the reads, of course.

Hi Brian and Geno, I didn't want to start a new thread, so I figured I'd quote and post to get your attention.

I started using BBDeDupe and Normalization in Geneious on Merged Paired Reads. These reads are not yet assembled, but will be with De Novo after DeDupe/Norm. I'm doing it because I have about 50,000 reads that cover a 10kb PCR amplicon, and I think the depth is causing assembly issues.

I wasn't able to find much information on DeDupe. Does it remove all 100% matching reads that are encompassed in another read, thereby reducing my coverage to 1? If so, is there a way to increase my coverage, to say 50-100 instead?

I'm getting a few areas with low coverage, less than 10, and I generally require that I have a minimum coverage of 15 to call a base in the consensus sequence, otherwise I call a gap. I generally have >500-1000 coverage in all places, but occasionally I get some places that have <10 coverage. I don't think this is amplicon bias, rather I think it is contaminating sequence. I use a core service for the illumina sequencing, but nearly everyone that submits samples is working on HIV, hence my need for a coverage threshold.

Thoughts on this?

Thanks,
Jake

[Brian Bushnell]

As long as the file is unsorted or name-sorted, that command will work. Spliced reads shouldn't cause a problem, unless the splice junction is in the unsequenced area between the reads, of course.

[GenoMax]

Is it possible to compute the insert size (and make a histogram distribution) from a .bam file ?

I think so (use unsorted bam). If reads are spliced you will get some odd results. @Brian may have a suggestion later on how to handle that.

Make sure samtools is available in your $PATH.

reformat.sh in=your_file.bam ihist=hist.txt

[Canadian_philosophy]

Hi Brian,
This command below worked for me. It looks like if I name the output files, that doesn't work (e.g out1_hg19=,out2_hg19=,out1_mm10= etc)
.
If I use the basename command, the fastqs look fine!

java -Djava.library.path=/data/davelab/dr153/softwares/bbmap/jni/ -ea -Xmx56G -cp /data/davelab/dr153/softwares/bbmap/current/ align2.BBSplitter ow=t fastareadlen=500 minhits=1 minratio=0.56 maxindel=20 qtrim=rl untrim=t trimq=6 in1=/data/davelab/projects/Xenomousie/fastq_link/SL118119_S5_L001_R1_001.fastq.gz in2=/data/davelab/projects/Xenomousie/fastq_link/SL118119_S5_L001_R2_001.fastq.gz interleaved=True basename=o5%_#.fq ref=/data/davelab/bix/resources/genomes/hg19/ucsc.hg19.fasta,/data/davelab/bix/resources/genomes/mm10/mm10_gatkorder.fa refstats=refstats.txt ambiguous2=toss

[Brian Bushnell]

Hmmm, that command looks fine. Just to clarify, is that the specific command that produced sam output?

[Canadian_philosophy]

Hmmm.. so I pasted that from the stderr/out file of the program, maybe that is how you formatted it there? This is how my command looks

java -Djava.library.path=/data/davelab/dr153/softwares/bbmap/jni/ -ea -Xmx56G -cp /data/davelab/dr153/softwares/bbmap/current/ align2.BBSplitter ow=t fastareadlen=500 minhits=1 minratio=0.56 maxindel=20 qtrim=rl untrim=t trimq=6 in1=/data/davelab/projects/Xenomousie/fastq_link/SL118119_S5_L001_R1_001.fastq.gz in2=/data/davelab/projects/Xenomousie/fastq_link/SL118119_S5_L001_R2_001.fastq.gz interleaved=True basename=o5%_#.fq ref=/data/davelab/bix/resources/genomes/hg19/ucsc.hg19.fasta,/data/davelab/bix/resources/genomes/mm10/mm10_gatkorder.fa refstats=refstats.txt ambiguous2=toss

The only comma is in between ref files. Does this look okay to you?

Also, thanks for fixing the default input.