[Brian Bushnell]

Woody,

Yep, the problem was error-correcting the two files independently, which will reorder the reads on a multithreaded machine. So, this command will solve it:

ecc.sh threads=8 -Xmx14g in=~/reads/n_R1.fastq in2=~/reads/n_R2.fastq out=~/reads/n_eccR1.fastq out2=~/reads/n_eccR2.fastq

[woodydon]

My command was:

~/bbmap/ecc.sh threads=8 -Xmx14g in=~/reads/n_R1
.fastq out=~/reads/n_eccR1.fastq

Ecc.sh's log:

Settings:
threads:                8
k:                      31
deterministic:          false
toss error reads:       false
passes:                 1
bits per cell:          16
cells:                  3147.19M
hashes:                 3
prefilter bits:         2
prefilter cells:        13.56B
prefilter hashes:       2
base min quality:       5
kmer min prob:          0.5

target depth:           40
min depth:              6
max depth:              40
min good kmers:         15
depth percentile:       54.0
ignore dupe kmers:      true
fix spikes:             false

Changed from ASCII-33 to ASCII-64 on input quality 94 while prescanning.
Made prefilter:         hashes = 2       mem = 3.15 GB          cells = 13.54B 
       used = 3.739%
Made hash table:        hashes = 3       mem = 5.86 GB          cells = 3144.33M   
     used = 6.197%

Estimated kmers of depth 1-3:   203133322
Estimated kmers of depth 4+ :   54945575
Estimated unique kmers:         258078898

Table creation time:            1887.618 seconds.
Started output threads.
Table read time:                1510.451 seconds.       6757.71 kb/sec
Total reads in:                 136095860       100.000% Kept
Total bases in:                 10207189500     100.000% Kept
Error reads in:                 12059846        8.861%
Error type 1:                   3939140         2.894%
Error type 2:                   6350473         4.666%
Error type 3:                   9031666         6.636%

During Error Correction:
Errors Suspected:               31378616
Errors Corrected:               15385594
Errors Marked:                  0

Total kmers counted:            6085312047
Total unique kmer count:        408064601

Includes forward kmers only.
The unique kmer estimate can be more accurate than the unique count, if the tables are
 very full.
The most accurate value is the greater of the two.

Percent unique:                  6.71%
Depth average:                  14.91   (unique kmers)
Depth median:                   1       (unique kmers)
Depth standard deviation:       535.10  (unique kmers)

Depth average:                  6532.20 (all kmers)
Depth median:                   2141    (all kmers)
Depth standard deviation:       26317.70        (all kmers)

Approx. read depth median:      3568.33

Total time:                     3399.627 seconds.       3002.44 kb/sec

Hope it helps.

Woody

[Brian Bushnell]

Woody,

That's very strange. The only explanation I can come up with is that the order of the reads is getting changed so read pairs are no longer always together.

If you ran ecc.sh twice on different files, that would make sense; unless you add the "ordered" flag, the reads will not necessarily come out in the same order. But it's best to correct all the reads together so that there is more coverage. So, assuming your reads are in two files (rather than interleaved), you should run it like this:

ecc.sh in1=read1.fq in2=read2.fq out1=corrected1.fq out2=corrected2.fq

And if your reads are interleaved in 1 file, you can add the "int=t" flag to ensure the reads are treated as interleaved. This is autodetected but the autodetection will fail if the read names are different than expected.

ecc.sh in=reads.fq out=corrected.fq int=t

If you do not think that was the problem, then please post the command line and standard out / standard error output of ecc. Thanks!

[woodydon]

Woody,

bbmerge does not do error-correction; for that, you need to use the script "ecc.sh", e.g.

ecc.sh in=reads.fq out=corrected.fq

And as Genomax implied, all of my programs are sensitive to file name extensions. If you say "out=reads.fasta" you'll get fasta; if you say "out=reads.sam.gz", you'll get gzipped sam, and so forth.

Hi Brain,

I used ecc.sh to correct my raw .fq files and then use bbmerge.sh again. I found the assemble rate is much lower:

Reads:          136095860
Joined:         13084199        9.614%
Ambiguous:      31319           0.023%
No Solution:    122968112       90.354%
Too Short:      12230           0.009%
Avg Insert:                     120.9
Standard Deviation:             12.9

Insert range:           26 - 138
90th percentile:        134
50th percentile:        123
10th percentile:        107

Without ecc.sh, the joined reads were about 30~40%. Now, it is only 9%...

Woody

[woodydon]

Woody,

bbmerge does not do error-correction; for that, you need to use the script "ecc.sh", e.g.

ecc.sh in=reads.fq out=corrected.fq

And as Genomax implied, all of my programs are sensitive to file name extensions. If you say "out=reads.fasta" you'll get fasta; if you say "out=reads.sam.gz", you'll get gzipped sam, and so forth.

I see~ thanks for GenoMax's and your explanations!

Woody

[Brian Bushnell]

I can't find error-correcting option in bbmerge.sh. Also, in bbmap.sh, "out=" seems outputs reads instead of sam files. Could you please give me a hint? Thanks!

Woody,

bbmerge does not do error-correction; for that, you need to use the script "ecc.sh", e.g.

ecc.sh in=reads.fq out=corrected.fq

And as Genomax implied, all of my programs are sensitive to file name extensions. If you say "out=reads.fasta" you'll get fasta; if you say "out=reads.sam.gz", you'll get gzipped sam, and so forth.

[GenoMax]

I can't find error-correcting option in bbmerge.sh. Also, in bbmap.sh, "out=" seems outputs reads instead of sam files. Could you please give me a hint? Thanks!

With bbmap.sh out=filename.sam should produce a sam file (http://seqanswers.com/forums/showpos...69&postcount=1).

[woodydon]

Woody,

Did you try error-correcting and using the 'vloose' setting? Both of those can substantially improve the merge rate.

If you want to see what the actual, unbiased insert size distribution of your data, then you can generate a histogram by mapping:

bbmap.sh in1=read1.fq in2=read2.fq ref=reference.fa ihist=ihist.txt nodisk

Then plot a graph of 'ihist.txt', so you can see how many reads should be expected to assemble.

I can't find error-correcting option in bbmerge.sh. Also, in bbmap.sh, "out=" seems outputs reads instead of sam files. Could you please give me a hint? Thanks!

[woodydon]

Hi Brian,

Thank you for the quick reply. It indeed was the problem of not finding java. I have consulted with the administrator about this since I used java before without a problem.

Sincerely,
Woody

[Brian Bushnell]

Woody,

It sounds like Java is not installed on your cluster. Try running this command:

java -version

It should print some kind of message, like this:

java version "1.6.0_07"
Java(TM) SE Runtime Environment (build 1.6.0_07-b06)
Java HotSpot(TM) Client VM (build 10.0-b23, mixed mode)

If it just says "command not found" then java is not installed. You should contact the system administrator, or download it from here:

http://www.oracle.com/technetwork/java/javase/downloads/jre8-downloads-2133155.html

...and put 'java' in your path. Or explicitly change the path of 'java' in my shellscript to wherever you install it.

[woodydon]

Hi Brain,

I am using a computing cluster to run BBMerge and found the following error:

java -ea -Xmx48g -cp /home/xxx/tools/bbmap/current/ align2.BBMap build=1 overwrite=true fastareadlen=500 -Xmx48g in1=/home/xxx/reads/n_R1.fastq in2=/home/xxx/reads/n_R2.fastq ref=/home/xxx/ref/hg19c/genome_hg19.fa ihist=/home/xxx/bbmap/n_ihist.txt out=/home/xxx/bbmap/n_merged.fq
/home/xxx/tools/bbmap/bbmap.sh: line 176: java: command not found

my command was:

/home/xxx/tools/bbmap/bbmap.sh -Xmx48g threads=10 in1=/home/xxx/reads/n_R1.fastq in2=/home/xxx/reads/n_R2.fastq ref=/home/xxx/ref/hg19c/genome_hg19.fa ihist=/home/xxx/bbmap/n_ihist.txt out=/home/xxx/bbmap/n_merged.fq

Thanks,
Woody

[woodydon]

Woody,

Did you try error-correcting and using the 'vloose' setting? Both of those can substantially improve the merge rate.

If you want to see what the actual, unbiased insert size distribution of your data, then you can generate a histogram by mapping:

bbmap.sh in1=read1.fq in2=read2.fq ref=reference.fa ihist=ihist.txt nodisk

Then plot a graph of 'ihist.txt', so you can see how many reads should be expected to assemble.

Great. That's an neat option.

[Brian Bushnell]

Thanks for the reply. It turns out that I confused your "insert size" with "inner distance" of pair-end reads. Since the ~25 million reads were exome-seq data, 44.594% assemble rate is a bit lower than what I expected. If the assemble is perfect, it means 50% of the fragments were actually longer than 188. I hope you could improve those fragments whose "insert size" is longer than 188 in the future.

Woody,

Did you try error-correcting and using the 'vloose' setting? Both of those can substantially improve the merge rate.

If you want to see what the actual, unbiased insert size distribution of your data, then you can generate a histogram by mapping:

bbmap.sh in1=read1.fq in2=read2.fq ref=reference.fa ihist=ihist.txt nodisk

Then plot a graph of 'ihist.txt', so you can see how many reads should be expected to assemble.

[woodydon]

Woody,

2x100bp reads can have a 74bp insert or a 188bp insert size. 74bp insert means that the molecule being sequenced was shorter than read length, and as a result the data collection continued off the end of the genomic sequence and into the adapter. So, before merging, the reads each contained 26bp of junk. And a 188bp insert size means that the reads overlapped by 12 base pairs. BBMerge does not look for overlaps shorter than 12bp in the default mode; the shorter the overlap, the more likely that it's a false positive.

-Brian

Thanks for the reply. It turns out that I confused your "insert size" with "inner distance" of pair-end reads. Since the ~25 million reads were exome-seq data, 44.594% assemble rate is a bit lower than what I expected. If the assemble is perfect, it means 50% of the fragments were actually longer than 188. I hope you could improve those fragments whose "insert size" is longer than 188 in the future.

[Brian Bushnell]

Ive been using it over the last couple of weeks and it has been really quick. Im actually using it to run adapter trimming (with small inserts, all of your reads should overlap), so adding in some options for adapter matching would be handy (although I do know that you already have a trimming script in the bbtools suite)

jimmybee,

I added some relevant flags and modes, in version 32.32. Thanks for the suggestions!

First, BBMerge now has a new flag, "tbo" (trimbyoverlap). Rather than joining reads, it will just trim the adapter sequence off the ends of reads that have an apparent insert size shorter than read length.

Second, I added the ability to overlap to BBDuk also, with the same flag, "tbo". So if you are trimming adapters fragment libraries with the normal read orientation, you can trim both by specifying an adapter sequence file and by overlap at the same time, which drastically reduces the false negative rate (reads that did not have adapters completely removed) by around a factor of 12. I highly recommend using both for fragment libraries, as the combination outperforms either one alone:

bbduk.sh in=reads.fq out=trimmed.fq tbo ktrim=r ref=truseq.fa k=25 mink=1 hdist=1

The 'tbo' flag is NOT recommended for long mate pair libraries which may have the adapters in totally different places and are not expected to overlap, which is why it is disabled by default.

-Brian