[rc630]

Sounds like the files are probably misformatted. Can you run "head 1_S1_L001_R1_001.fastq" and "head 1_S1_L001_R2_001.fastq" and post the output here?

head 1_S1_L001_R1_001.fastq
@M01936:151:000000000-AL2M3:1:1101:16824:1851 1:N:0:1
TCGATCGCTTGGTCATTTAGAGGAAGTAAAAGTCGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCACTAGAGATATCGCCCTTCGGGGCTCTCATCCCTTCTATCACACCCCTGTGCACTTTGGCCACCGCTTCATTGCGGTGTGTCTT
+
BCCACFCCFCCCGGGGGGGGGGGHGHHHGHGHHHGGGAE53CGHHHHFGHGGGHHHGHHHHHGGGGGGHHHHHFB3335FHHHGGGAEHFA11E/EGHHHHHHGHGFHHHHFHHFDEDFAFGHHHHHHFEGCFFCFGGGGHHHGFFGDDDDHHHH
@M01936:151:000000000-AL2M3:1:1101:18467:1911 1:N:0:1
TCGATCGCTTGGTCATTTAGAGGAAGTAAAAGTCGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCACTAGAGATATCGCCCTTCGGGGCTCTCATCCCTTCTATCACACCCCTGTGCACTTTGGCCACCGCTTCATTGCGGTGTGTCTT
+
BCCDCFCCFDCDGGGGGGGGGGHHHHHHHHHHHHGGHHGHGHHHHHHHGHHGGHHHHHHHHHGFGGGGHHHHFHHHE35GHHHFGGGGHHEE1E/EHHHHHHFHHGHHHHHHHHHGGCGEGHHHHHHHHGHHHGHGGGGGHHHGHGGGGGGHHHH
@M01936:151:000000000-AL2M3:1:1101:14252:2010 1:N:0:1
TCGATCGCTTGGTCATTTAGAGGAAGTAAAAGTCGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCACTAGAGATATCGCCCTTCGGGGCTCTCATCCCTTCTATCACACCCCTGTGCACTTTGGCCACCGCTTCATTGCGGTGTGTCTT

1_S1_L001_R2_001.fastq 
@M01936:151:000000000-AL2M3:1:1101:16824:1851 2:N:0:1
CGGACTTGATGTACGAGCTGCGTTCTTCATCGATGCGAGAACCTAGAGATCCGTTGTTGAAAGTTATTTTAGTTTATAAGTTTTTACATTCAATGACTTGTGTTTATAGGTATGGTATAGTAAAAAGACACACCGCAATGAAGCGGTGGCCAAAG
+
33>>ABFFFFFFGGGGGGGGGFGGGGHGGHGFBGHGGCGFGB2B3333DFFEGFH1AFF555BFHHHHHG53@BF5FGGHHHHHEHDGHHHHBHHFHGGHHFHHHBFHFFHHHHFGFHEGF?GFFGEFHGFFFECCEGDFFGHGG?FCDCFG/FG
@M01936:151:000000000-AL2M3:1:1101:18467:1911 2:N:0:1
CGGACTCGATGTACGAGCTGCGTTCTTCATCGATGCGAGAACCTAGAGATCCGTTGTTGAAAGTTATTTTAGTTTATAAGTTTTTACATTCAATGACTTGTGTTTATAGGTATGGTATAGTAAAAAGACACACCGCAATGAAGCGGTGGCCAAAG
+
3>>AABFCBFFFGGGGGGGGGGGGGGHHHHGHHHHGGGGGGGFHGH335GHHHFHCGHGHG5B5DHHHHHHFHFHGGHGHFHHHGHHHFFHHFHHHHHHHHFHHHHHHHGFGHHHHHHHHBGFHHGEGHHHHGFG?DGCHHHEHGCGGGFHHAGB
@M01936:151:000000000-AL2M3:1:1101:14252:2010 2:N:0:1
CGGACTTGATGTACGAGCTGCGTTCTTCATCGATGCGAGAACCTAGAGATCCGTTGTTGAAAGTTATTTTAGTTTATAAGTTTTTACATTCAATGACTTGTGTTTATAGGTATGGTATAGTAAAAAGACACACCGCAATGAAGCGGTGGCCAAAG

[Brian Bushnell]

Sounds like the files are probably misformatted. Can you run "head 1_S1_L001_R1_001.fastq" and "head 1_S1_L001_R2_001.fastq" and post the output here?

[rc630]

"Created read stream of 0 reads"

I am trying to use BBMerge with the two files I have for one sample of MiSeq data. I am getting this error however:

"Executing jgi.BBMerge [in1=1_S1_L001_R1_001.fastq, in2=1_S1_L001_R2_001.fastq, out=S1merged, reads, outu1=S1unmerged1, outu2=S1unmerged2]

BBMerge version 8.9
Writing mergable reads merged.
Unspecified format for output S1merged; defaulting to fastq.
Unspecified format for output S1unmerged1; defaulting to fastq.
Unspecified format for output S1unmerged2; defaulting to fastq.
Started output threads.
crisG: Warning - created a read stream for 0 reads.
Exception in thread "main" java.lang.AssertionError
at stream.ConcurrentGenericReadInputStream.<init>(ConcurrentGenericReadInputStream.java:129)
at stream.ConcurrentReadInputStream.getReadInputStream(ConcurrentReadInputStream.java:129)
at stream.ConcurrentReadInputStream.getReadInputStream(ConcurrentReadInputStream.java:51)
at jgi.BBMerge.runPhase(BBMerge.java:957)
at jgi.BBMerge.process(BBMerge.java:721)
at jgi.BBMerge.main(BBMerge.java:48)"

Why might it be saying there are 0 reads? I have opened it with FastQC and there are reads in the file.

Thanks

[Brian Bushnell]

For interleaved reads:

bbmerge.sh in=reads.fq out=merged.fq outu=unmerged.fq

For reads in two files:

bbmerge.sh in1=read1.fq in2=read2.fq out=merged.fq outu1=unmerged1.fq outu2=unmerged2.fq

[tarias]

Hi,
I have been using this program to merge reads but have not been able to generate a file with only the unmerged reads, any ideas on how to do this? thanks

[Brian Bushnell]

That's an excellent question. I've never used Trinity. A post in their forum claimed that it will ignore the extra reads if you concatenate them to the end of one of the paired files. However, their faq claims the opposite:

If you have additional singletons, add them to the .fq file that they correspond to based on the sequencing method used (if they’re equivalent to the left.fq entries, add them there, etc).

So it sounds like probably you should, indeed, concatenate them at the end of the left reads. You can easily test this by using a single pair of reads (one read in the left and one in the right file), then concatenating all the unmerged reads onto the left file, then assembling. If if finishes almost immediately and nothing assembles, then it's ignoring the merged reads

[jazz710]

How to use merged reads

Brian,

Cool tool (as always)! I'm loving the Geneious integration for these tools. A+ move on your part.

I'm performing my first merges on some raw and trimmed datasets and I had a very basic question.

After merging my PE data, I will be left with some merged reads, and some unmerged (original PE) reads, correct? For de novo assembly, would I simply concatenate the merged reads onto the bottom of the 'Left' unmerged reads, then run Trinity with L (PE + merged) and R (PE R only)?

Best,
Bob

[Brian Bushnell]

BBMerge now has a new flag - "outa" or "outadapter". This allows you to automatically detect the adapter sequence of reads with short insert sizes, in case you don't know what adapters were used. It works like this:

bbmerge.sh in=reads.fq outa=adapters.fa reads=1m

Of course, it will only work for paired reads! The output fasta file will look like this:

>Read1_adapter
GATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG
>Read2_adapter
GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGATGTATCTCGTATGCCGTCTTCTGCTTG

If you have multiplexed things with different barcodes in the adapters, the part with the barcode will show up as Ns, like this:

GATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG

[pmiguel]

Hi Phillip,

[...]

There is probably a better way to derive the quality than the simple equation I am using, but I like it because it is simple and gives useful results that can generally be used by tools that are designed for raw Illumina quality scores. Deriving an equation that correctly models all factors, or does a substantially better job overall, while still being simple enough that a person can understand the relationship between input and output, would be extremely difficult, I think.

Well, you've convinced me! Thanks for the response, Brian.

--
Phillip

[Brian Bushnell]

Hi Brian,
That is interesting. Any reason you chose to add only 1/4th the quality value of one of the pair? The old phrap method was to just add the quality values together if bases from different strands agreed.

That seemed reasonable enough -- if chance of one basecall being wrong was 1 in 100 (QV 20) but the basecall from the other strand agreed and also was QV 20, then the chance they are both wrong seems like it should be 1 in 10,000 (QV 40). Right?

--
Phillip

Hi Phillip,

Mathematically, if you have two completely independent observations, with perfectly accurate quality scores, the correct formula would be to add the scores together if the bases agree, and subtract them if they differ. However, I don't really like that approach since the quality scores are not perfectly accurate and the observations are not completely independent. For one thing, if a cluster is close enough to another cluster that there is interference during read 1, there will also be interference during read 2, with a similar nonrandom bias. Or if a cluster is near the edge of a lane for read 1 and thus slightly out of focus, it will be for read 2 as well. So even if the quality scores are accurate, both are affected by a similar bias rather than unrelated random biases, and thus strictly adding them is not appropriate.

Furthermore, a lot of downstream programs or analyses may be calibrated to assume that the dominant source of error is sequencing reading error, and ignore other error modes, as they tend to be smaller. Thus Q40 reads may yield Q40 variants. But when you have a Q80 read, the error mode will be dominated by other things like PCR errors or unwanted chemical reactions - there's no way merging two Q40 reads will yield bases with a 1/100,000,000 error rate, even if you can absolutely ensure that the overlap frame was correct, which you can't.

There is probably a better way to derive the quality than the simple equation I am using, but I like it because it is simple and gives useful results that can generally be used by tools that are designed for raw Illumina quality scores. Deriving an equation that correctly models all factors, or does a substantially better job overall, while still being simple enough that a person can understand the relationship between input and output, would be extremely difficult, I think.

[HESmith]

Hi Brian,

Thanks for the prompt response. Given the lower q scores of disconcordant basecalls, I can easily filter those from consideration.

Take care,
Harold

[pmiguel]

Hi Harold,
When there is a mismatch, the base chosen is the one with the higher quality value, or N if they are equal. The resulting quality value of overlapping bases b1 and b2 with qualities q1 and q2:

if b1 equals b2:
q=max(q1, q2)+min(q1, q2)/4
(capped at 41 by default, though this can be overridden)

Hi Brian,
That is interesting. Any reason you chose to add only 1/4th the quality value of one of the pair? The old phrap method was to just add the quality values together if bases from different strands agreed.

That seemed reasonable enough -- if chance of one basecall being wrong was 1 in 100 (QV 20) but the basecall from the other strand agreed and also was QV 20, then the chance they are both wrong seems like it should be 1 in 10,000 (QV 40). Right?

--
Phillip

[Brian Bushnell]

Hi Harold,

BBMerge will reject overlaps that have too many mismatches; by default that number is 3, though you can override it. When there is a mismatch, the base chosen is the one with the higher quality value, or N if they are equal. The resulting quality value of overlapping bases b1 and b2 with qualities q1 and q2:

if b1 equals b2:
q=max(q1, q2)+min(q1, q2)/4
(capped at 41 by default, though this can be overridden)

if b1 does not equal b2:
q=max(q1, q2)-min(q1, q2)

As a result, the quality value will be very low if the bases did not match, so any resulting erroneous variant calls would also be very low quality, and thus should be easy to filter.

There is currently no provision for using a reference base where the reads don't match, as BBMerge does not take the reference into account.

[HESmith]

A colleague is attempting to detect low-frequency (<0.1%) mutations in plasmid amplicons by PE sequencing and error correction in the overlap. I thought BBMerge might be suitable for this application, and have a couple of questions:

1) how does BBMerge treat base-call discrepancies in the overlap?
2) for discrepancies, is there a method to impose the reference base on the merge, assuming one of the two reads matches the reference? I suppose this could be handled at the alignment or consensus-building step, if the discrepancies are flagged in some way.

Thanks Brian,
Harold

[Marisa_Miller]

Hi Marisa,

BBMerge will only calculate the insert size distribution for paired reads, but the insert size of the merged reads will be equal to their actual length. So you can calculate the distribution after merging like this:

readlength.sh in=merged.fastq out=hist.txt

-Brian

Thanks Brian!