[Brian Bushnell]

1) I have not looked at Nanopore reads in several years now so I don't have any further experience in mapping them. But from what I read and hear the quality is gradually creeping upward (and my initial exposure was to very low quality data) so some of the settings can probably be relaxed; e.g., increasing to k=10 and minratio=0.25.

2) "Read Depth" (DP field) is the depth of all reads covering (or spanning) the event, while "Allele Depth" (AD field) is the same but only for those reads exhibiting the variation.

Coverage is added to an array tracking per-reference-base coverage at each position when each read is processed (and deleted bases are considered covered by a read containing that deletion). Similarly, coverage is added to each variation when a read containing that variation is processed.
The reference coverage for a variant is calculated by summing the coverage per-base in the array over the length of the event. For substitutions this is straightforward. For deletions, this is the deleted reference bases. For insertions, this is the average coverage of the bases on either side of the insertion.
It is possible for the allele coverage to be slightly higher than the total coverage using this calculation method (for example, as a result of a sudden drop in coverage on one side of an insertion event) so DP reports max(ref coverage, allele coverage) while COV reports the raw calculated coverage (they will virtually always be the same, though). AF is calculated as DP/AD to prevent AF over 1.0.

3) "Allele Fraction" is "Allele Depth"/"Read Depth". As you say, this is straightforward for SNPs but less so for indels, particularly with short reads and long deletions. It's not entirely clear what AF should really mean for a deletion, because it describes a ratio of two classes of reads - those that cover a variant location and do not indicate the variant, and those that cover a variant location and do indicate the variant. But when deletions are longer than reads there is a third case; those that land in the middle of the deletion but do not touch either end. You could, for example, have 100x coverage on either end of a deletion event, with 50x indicating the deletion and 50x indicating no deletion, and calculate an AF of 0.5. But if you had 1000x coverage somewhere in the middle of the deletion event from reads not spanning the junction on either end, that's a 3rd category and really they should not contribute to the AF (since the deleted allele would never be observed in that case), but currently they do. I'll look into changing the way that long deletion AF and DP are calculated; thanks for making me think about it.
For insertions, there is "Alternate Allele Frequency" which is allele frequency adjusted to compensate for the fact that insertions will not be called on reads that do not span the full event. So if you have a 50bp insertion and 100bp reads, only 1/3rd of the reads with bases in the inserted region would actually span the full insertion and thus possibly yield an insertion call, so a 100% true allele frequency might be reported as 33%. "Alternate Allele Frequency" would adjust the 33% observed back up to 100% based on the relative average read length and insertion length.

[jweger1988]

Hi Brian,

I have a couple of questions about using callvariants.sh and mapPacbio.sh. I am using mapPacbio.sh to align reads from a Nanopore minion to a viral sequence and looking for deletions with callvariants.sh. My questions are

1. Do you have a suggested method for mapping minion reads. I used a suggestion you provided in 2014 but I don't know if you have since updated this. Here is what I'm running.

mapPacBio.sh k=8 in=reads.fq ref=ref.fa maxreadlen=1000 minlen=200 sam=1.4 idtag ow int=f qin=33 indelhist=indelhist1.txt out=mapped.sam minratio=0.15 ignorequality slow ordered maxindel=12000 outputunmapped=f bs=bs1.sh

2. I am then calling variants with the desire to look for deletions with callvariants.sh.

callvariants.sh in=mapped.sam rarity=0.0001 out=deletions.tsv vcf=deletions.vcf ref=ref.fa extended=t border=20 minscore=10

In the output can you tell me what the meaning of "coverage" or "read Depth" or "allele depth" are? What is the difference between these? Also, what does the allele fraction mean in this case? For a SNP it's self-explanatory but for a indel it seems like it's not as straight forward.

Any insight you can provide would be extremely useful. Thanks again for all you do.

[Brian Bushnell]

Hi Dario and Genomax,

I though I had gotten rid of that a while ago, but I guess not - I'll investigate and fix it. It's harmless and due to a race condition for a thread finishing after it was prematurely shut down because enough reads were generated, but it looks scary.

[GenoMax]

I would like to add that reformat.sh also exhibits this "issue".

$ reformat.sh in=P1.fastq.gz out=stdout.fa reads=5

Input:                          5 reads                 1094 bases
Output:                         5 reads (100.00%)       1094 bases (100.00%)

Time:                           0.178 seconds.
Reads Processed:           5    0.03k reads/sec
Bases Processed:        1094    0.01m bases/sec
Exception in thread "main" java.lang.RuntimeException: ReformatReads terminated in an error state; the output may be corrupt.
        at jgi.ReformatReads.process(ReformatReads.java:1148)
        at jgi.ReformatReads.main(ReformatReads.java:45)

[dariober]

bbduk.sh: "reads" parameter does not exit clean

Hi- Not sure this thread is the best place to post this...

It seems to me that the `reads` parameter makes bbduk exits with an error when the input is not fully processed.

These runs are both ok (input file has less then 1M reads):

bbduk.sh in=test.fq.gz out=test.out.fq reads=-1 overwrite=true
bbduk.sh in=test.fq.gz out=test.out.fq reads=1000000 overwrite=true

However reads=10 gives:

bbduk.sh in=test.fq.gz out=test.out.fq reads=10 overwrite=true

java -Djava.library.path=/home/db291g/applications/BBmap/bbmap/jni/ -ea -Xmx15023m -Xms15023m -cp /home/db291g/applications/BBmap/bbmap/current/ jgi.BBDukF in=test.fq.gz out=test.out.fq reads=10 overwrite=true
Executing jgi.BBDukF [in=test.fq.gz, out=test.out.fq, reads=10, overwrite=true]

BBDuk version 37.54
NOTE: No reference files specified, no trimming mode, no min avg quality, no histograms - read sequences will not be changed.
Initial:
Memory: max=15097m, free=14703m, used=394m

Input is being processed as paired
Started output streams:	0.082 seconds.
Processing time:   		0.017 seconds.

Input:                  	20 reads 		3011 bases.
Total Removed:          	0 reads (0.00%) 	0 bases (0.00%)
Result:                 	20 reads (100.00%) 	3011 bases (100.00%)

Time:   			0.114 seconds.
Reads Processed:          20 	0.18k reads/sec
Bases Processed:        3011 	0.03m bases/sec
Exception in thread "main" java.lang.RuntimeException: jgi.BBDukF terminated in an error state; the output may be corrupt.
	at jgi.BBDukF.process(BBDukF.java:913)
	at jgi.BBDukF.main(BBDukF.java:72)

echo $?
1

I guess the output is still ok but the non-zero exit code causes problems inside a pipeline.

Am I missing something here or should this be fixed?

Thanks

Dario

[boulund]

Yes, that's correct. This is especially important for BBMerge, where insert size can only be computed from reads that overlap; so if the average calculated insert size is say 270bp for 2x150bp reads, but "PercentOfPairs" is only 20%, then presumably the actual insert size is substantially higher but those pairs did not overlap so it could not be computed. It's less important for BBMap but still useful in some cases - if PercentOfPairs is very low, something probably went wrong and the insert size is not trustworthy. For example, if you have an extremely fragmented reference with contig length shorter than insert size, such that most paired reads map to different contigs, then the pairing rate will be low and the calculated insert size will be untrustworthy (shorter than the actual insert size).

Hi, thanks for that excellent response!
I'm one of the authors of a BBMap module for the MultiQC bioinformatics log aggregation/summarization tool, and we'd love your input/help with small snippets of explanatory text for some of the outputs produced by the different BBMap tools. Would you be willing to help out? The first reasonably finished version of the BBMap module for MultiQC was finished and merged yesterday, but it unfortunately still lacks explanatory help text for most plot types. Information at the level of your reply regarding PercentOfPairs is perfect for the help text in the MultiQC reports, to aid people unfamiliar with BBMap output to interpret the plotted results included in a MultiQC report. Perhaps you even have suggestions on how to better visualize aggregations of some of the different outputs. Let me know what you think, and we'll set up an Issue in the MultiQC Github issue tracker to organize such work if you're interested.

[Brian Bushnell]

Hi,

I have a question about the PercentOfPairs output found in the insert size histogram output (ihist). What does it mean, percent of pairs? Is it the proportion of read pairs where an insert size could be computed?

Yes, that's correct. This is especially important for BBMerge, where insert size can only be computed from reads that overlap; so if the average calculated insert size is say 270bp for 2x150bp reads, but "PercentOfPairs" is only 20%, then presumably the actual insert size is substantially higher but those pairs did not overlap so it could not be computed. It's less important for BBMap but still useful in some cases - if PercentOfPairs is very low, something probably went wrong and the insert size is not trustworthy. For example, if you have an extremely fragmented reference with contig length shorter than insert size, such that most paired reads map to different contigs, then the pairing rate will be low and the calculated insert size will be untrustworthy (shorter than the actual insert size).

[boulund]

Hi,

I have a question about the PercentOfPairs output found in the insert size histogram output (ihist). What does it mean, percent of pairs? Is it the proportion of read pairs where an insert size could be computed?

[luc]

Hi Brian,

thanks for the details -- there is always a another twist to it.

Hi Luc,

To make this very clear, it's the estimate of the size of sequenced genomic molecule. So, for example, if you have 2 reads like this, where R1 maps to the plus strand and R2 maps to the minus strand:

   111111111111111
-------------------------------------------------------------------
                                          22222222222222222

...then it would generally be the distance from the leftmost 1 to the rightmost 2. But if either read contains insertions, those are added to the insert size, and deletions are subtracted from the insert size. Otherwise, the insert sizes would be dramatically overestimated in spliced RNA-seq data.

[Brian Bushnell]

That's correct for single-ended reads, but I do not recommend setting "secondarysitescoreratio = 1.0" for paired reads because one pairing is currently slightly favored over others due to the difficulty of tracking all possible pairings. So, just set "secondary=t ambig=all maxsites=100". In other words, it's not currently possible to ensure you get exactly and only all reads that have the exact same (best) alignment score unless they are processed unpaired.

[darthsequencer]

"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.

Hi - So if I set secondarysitescoreratio to 1.0 do I need to set secondary = T if ambiguous=all and increase the maxsites setting?

I guess what I'm really asking is what do I need to set to ensure that I'm getting all mapping locations of read with all the same top score.

[Brian Bushnell]

Hi Luc,

To make this very clear, it's the estimate of the size of sequenced genomic molecule. So, for example, if you have 2 reads like this, where R1 maps to the plus strand and R2 maps to the minus strand:

   111111111111111
-------------------------------------------------------------------
                                          22222222222222222

...then it would generally be the distance from the leftmost 1 to the rightmost 2. But if either read contains insertions, those are added to the insert size, and deletions are subtracted from the insert size. Otherwise, the insert sizes would be dramatically overestimated in spliced RNA-seq data.

[luc]

I am sorry I can't find the info - I am sure it is somewhere.

How does BBmap define insert size ( e.g of the insert size histograms)?
Distance between start point of the reads (forward and reverse) or between the ends of the reads?

#############
I just checked with a test file; as expected the insert size seems to be the distance between the start point of the reads.

[boulund]

Oh, that was easy. Thanks a lot for your quick response, both of you. Sorry for not spotting that option!

@boulund: Can you include the exact starswrapper.sh command you are running?

The command I used was simple:

statswrapper.sh in=assembly_1.fasta in=assembly_2.fasta

With the "format=3" option it produces exactly the output I was expecting.

[Brian Bushnell]

Hi Boulund,

Please add the flag "format=3" to get one line per file.