[Thias]

Some questions regarding blacklisting

Hello Brian,

Sorry for bothering you again, but I have a question regarding blacklisting and its (correct) usage.

As a favor for a collegue, I quickly mapped some genomic PE data (ChIP-seq) to the human genome. The mapping statistics looked good and I handed the data over to him didn't give it a moment's thought anymore.

However he assessed the coverage profiles and now suspects a contamination of his library with DNA fragments originating from certain BACs, which are used in his laboratory. He showed me some example chromosomal regions and indeed, there are are several ~150kb big regions with remarkable increase in coverage overlapping the annotated BAC.

I offered to exclude reads from those BAC areas of the genome and try to identify real ChIP-signal somewhere else in the data.

Thus I restarted the alignment from FASTQ and specified a Multi-FASTA file containing the genomic sequences of the BACs as blacklist. bbmap indeed now generated the specified outb sam file, but not a single read pair went there - for every sample the file was always empty except the header.

Genome:                 3
Key Length:             13
Max Indel:              16000
Minimum Score Ratio:    0.56
Mapping Mode:           normal
Reads Used:             148059914       (14335498222 bases)

Mapping:                51822.670 seconds.
Reads/sec:              1428.52
kBases/sec:             276.63


Pairing data:           pct reads       num reads       pct bases          num bases

mated pairs:             98.5211%        72935102        98.8062%        14164366056
bad pairs:                1.3022%          964043         1.2593%          180529886
insert size avg:          425.12


Read 1 data:            pct reads       num reads       pct bases          num bases

mapped:                  99.8531%        73921244        99.8812%         7174307831
unambiguous:             92.3162%        68341677        92.4104%         6637694291
ambiguous:                7.5369%         5579567         7.4708%          536613540
low-Q discards:           0.0000%               0         0.0000%                  0

perfect best site:       78.4885%        58104980        78.7959%         5659785119
semiperfect site:        78.5457%        58147328        78.8410%         5663027631
rescued:                  0.2712%          200743

Match Rate:                   NA               NA        97.9450%         7152758588
Error Rate:              21.3317%        15768653         2.0538%          149983878
Sub Rate:                19.9830%        14771656         0.2769%           20218648
Del Rate:                 1.1258%          832171         1.7599%          128524996
Ins Rate:                 0.8954%          661856         0.0170%            1240234
N Rate:                   0.0807%           59648         0.0012%              90361


Read 2 data:            pct reads       num reads       pct bases          num bases

mapped:                  99.8967%        73953500        99.9086%         7146116046
unambiguous:             92.1534%        68221148        92.2105%         6595499668
ambiguous:                7.7433%         5732352         7.6981%          550616378
low-Q discards:           0.0000%              13         0.0000%                470

perfect best site:       68.2412%        50518933        68.6293%         4908820958
semiperfect site:        68.3924%        50630855        68.7725%         4919060978
rescued:                  0.3138%          232341

Match Rate:                   NA               NA        96.7175%         7114752246
Error Rate:              20.8947%        15452343         3.1485%          231610349
Sub Rate:                19.4650%        14395026         0.2754%           20259503
Del Rate:                 1.4361%         1062021         2.8562%          210107455
Ins Rate:                 0.8087%          598072         0.0169%            1243391
N Rate:                  13.1625%         9734162         0.1340%            9860906

Total time:             53846.780 seconds

To cross-validate, I indexed the same Multi-FASTA file and used bbmap to align the reads with exactly the same parameter settings against it. This are the results from mapping the same sample as above against the BAC collection:

Genome:                 1
Key Length:             11
Max Indel:              16000
Minimum Score Ratio:    0.56
Mapping Mode:           normal
Reads Used:             148059914       (14335498222 bases)

Mapping:                17929.866 seconds.
Reads/sec:              4128.86
kBases/sec:             799.53


Pairing data:           pct reads       num reads       pct bases          num bases

mated pairs:              0.2606%          192935         0.2584%           37041616
bad pairs:                0.2608%          193052         0.2578%           36960674
insert size avg:          1540.38


Read 1 data:            pct reads       num reads       pct bases          num bases

mapped:                   2.9434%         2179014         2.8388%          203905833
unambiguous:              1.6845%         1247069         1.6340%          117369776
ambiguous:                1.2589%          931945         1.2048%           86536057
low-Q discards:           0.0000%               0         0.0000%                  0

perfect best site:        0.0508%           37607         0.0407%            2926546
semiperfect site:         0.0514%           38023         0.0409%            2939958
rescued:                  0.0097%            7213

Match Rate:                   NA               NA         6.7979%          181724888
Error Rate:              98.2727%         2141375        93.2020%         2491518976
Sub Rate:                98.2038%         2139875         0.7663%           20485726
Del Rate:                28.9776%          631426        92.3724%         2469340316
Ins Rate:                20.6261%          449445         0.0633%            1692934
N Rate:                   0.0796%            1735         0.0001%               2285


Read 2 data:            pct reads       num reads       pct bases          num bases

mapped:                   4.8730%         3607461         4.6617%          333437379
unambiguous:              2.7381%         2026992         2.6422%          188984411
ambiguous:                2.1349%         1580469         2.0196%          144452968
:

My questions now are:

• Shouldn't I have just by random some pairs going to the outb sam file ? (since this are BACs which contain fragments of the human genome, against which I align)
• Does the latter mapping results point into the direction of a BAC contamination at all? I am a bit puzzling, because of the low number of mated pairs.

[Brian Bushnell]

I just verified that "outb" and "blacklist" still work, but their usage may be a bit unclear. I don't use them anymore as they have been superseded by BBSplit which allows any number of reference-specific output streams. I think in your case you may not have built the BAC fasta into the index.

How to use the blacklist to remove e.coli contaminants from human reads:

Concatenate the references:
cat ecoli.fa human.fa > both.fa

index the concatenated file:
bbmap.sh ref=both.fa

map to the combined index and specify a blacklist:
bbmap.sh in=reads.fq outm=good.fq outb=bad.fq blacklist=ecoli.fa

The "blacklist" flag does not index anything; it only parses the headers in the fasta file, whose contents must have already been indexed. If a read can map to two places equally well, by default, BBMap will assign it to the first genomic location (unless you alter the "ambig" flag). That is why the blacklisted files must be first when concatenating files together.

The way you would do this with BBSplit is simpler:

bbsplit.sh ref=ecoli.fa,human.fa
bbsplit.sh in=reads.fq pattern=out_%.fq

That would create two output files, out_ecoli.fq and out_human.fq (or .sam if you want). You can determine the behavior of reads that map to both references equally well in bbsplit with the "ambiguous2" flag; default behavior is to go to the output stream for the first reference. The other options are "ambig2=toss" (discard anything that maps to multiple references equally well), "ambig2=all" (output it to the stream for each reference it maps to equally well), and "ambig2=split" (write the multi-mapping reads to separate files).

[Thias]

I think in your case you may not have built the BAC fasta into the index.

*Sigh*. You know your cardboard homies ("Pappenheimer" ) too well ...

I was indeed not rebuilding the index of the reference genome for that purpose, but just added the blacklist parameter.

Thanks a lot for the fast response, the detailed description of the procedure and the inherent logic.

The way you would do this with BBSplit is simpler:

bbsplit.sh ref=ecoli.fa,human.fa
bbsplit.sh in=reads.fq pattern=out_%.fq

Great! This possibility will also aid to solve another question in a different context. Thanks for pointing this out ...

[Brian Bushnell]

*Sigh*. You know your cardboard homies ("Pappenheimer" ) too well ...

That was an interesting reference; I enjoyed reading it! In fact, it spurred me to go see if Mount&Blade II was near completion, but alas, it still seems to be in permanent development.

Since I develop BBMap at a facility where people use it in unexpected ways, I'm used to getting error reports caused by both bugs and misunderstandings. Both are valuable! I don't like difficult-to-use tools, but I can only make BBMap more user-friendly with feedback, so please don't feel shy to give feedback whether or not there's a bug.

In other news - BBMap and Pileup now have an "rpkm" flag. This will directly report RPKM and FPKM counts (as well as coverage depth and raw count) for each reference sequence, in the same format as Seal. Essentially, it's for RNA-seq quantification, and I primarily made it to make it easier to compare BBMap with Seal (since the RPKM/FPKM output only makes sense with a transcriptome reference, not a genome reference) - but it's there if anyone wants it.

[rkanwar]

Strange behavior for some of the aligned PE reads

Hello Brian,

I am using BBMap to align DNASeq PE data [100x2 bp reads]. The alignment works fine. But I am seeing strange alignments for a few pairs. For examples consider the SAM records for the following PE read:

R0266248:292:C5C0TACXX:4:1101:6482:13924 1:N:0:CTGAAGCTTAATCTTA	153	1	16877330	44	101M	=	16877330	0	GAAAGCCCCAGTGCCATCAGACAGCCACACCTCATCCTCATCAGTGACACTATGAGGTGAAGACCCCTCCAGGGTGTCAAGAGCTCTCAGCTTCCAGGGTC	4DDDBDDDDDDDDDDEDDDDBDDDDCAACFFFFHHHHHHIGIJIFJIHGIIIJJJJJJJJIHIIIJJIJJJJJJJIIJJJIJJJJJJJHHHHHFFFFFCCC	RG:Z:AS_N+AS_N_1+FCC5C0TACXX+L4+ICTGAAGCT-TAATCTTA	AM:i:0	NM:i:0

R0266248:292:C5C0TACXX:4:1101:6482:13924 1:N:0:CTGAAGCTTAATCTTA	101	1	16877330	0	*	=	16877330	0	TCCTCAGCCAGCTTCCCTTCTAGTCAATTTTTTGGGAGCTGGGCCTCCAGCTGGGATACTCTCTGGATGAGACTCTTCAGGTCCTTTTTGGCCTGAAGTCC	CCCFFFFFHHHHHJJJJJJJJJIJJJGJJJJJJIJJJJJJJJJJJJJJJJJJJJJGHJJHHHHHHFFFFFEEEEECDDDDDCDDDDDDC?@CCDDDDDDDD	RG:Z:AS_N+AS_N_1+FCC5C0TACXX+L4+ICTGAAGCT-TAATCTTA

The first end aligns properly. The other end does not align and has an alignment score of 0, which I think means that there were multiple alignments for this end. When I BLAT these two sequences then this is what I get for the reads:

END1
browser details YourSeq          101     1   101   101 100.0%     1   +   16877330  16877430    101
browser details YourSeq           99     1   101   101  99.1%     1   -   17183501  17183601    101
browser details YourSeq           97     1   101   101  98.1%     1   +   17008345  17008445    101
browser details YourSeq           95     1   101   101  97.1%     1   +   21738088  21738188    101
browser details YourSeq           91     1   101   101  95.1%     1   +  144879169 144879269    101

END2
browser details YourSeq          101     1   101   101 100.0%     1   -   17183611  17183711    101
browser details YourSeq          101     1   101   101 100.0%     1   +   16877220  16877320    101
browser details YourSeq           93     1   101   101  96.1%     1   +  144879059 144879159    101
browser details YourSeq           92     1   101   101  97.0%     1   +   21737977  21738078    102
browser details YourSeq           74    26   101   101  98.7%     1   +   17008260  17008335     76
browser details YourSeq           25     9    34   101 100.0%    10   +   69768146  69768175     30

BBMap aligns the first end properly with a start position of 16877330. The other end has 2 perfect hits, one of those hits [16877220] is in the neighborhood of of the first end. This neighboring hit should have been reported as the alignment for this end as it best explained the whole fragment [my insert sizes are around 250bp]. Here is the command I am using. I want to get only those alignment where the PE reads aligns uniquely on the reference genome. Can you please help me to set the right flags to get these reads align properly ? Thanks!

bbmap_cmd="${bbmap} -Xms10G -Xmx31G \
           path=${bbmap_ref} \
           in=${bbduk_fq1} \
           in2=${bbduk_fq2} \
           out=stdout \
           ambiguous=toss \
           rgid=${ID} \
           rgsm=${SM} \
           rglb=${LIB} \
           rgpl=${PL} \
           pairlen=1000 \
           threads=12 \
           slow=t \
           local=t \
           reads=10000 \
           samversion=1.3"

regards,
Rahul

[Brian Bushnell]

Hi Rahul,

In this case, the problem is "ambiguous=toss". If you omit that flag, the second read will be mapped. You can always filter ambiguous mappings later by excluding anything with a mapq of under 4.

However, this situation is a bit strange. your BLAT alignments show both reads hitting the same strand. Did you reverse-complement it prior to BLAT, or something? By default, BBMap requires read 1 and read 2 to be on opposite strands, so if they really are both mapping to the same strand, they won't be considered a proper pair. You can override that with "rcs=f" (requirecorrectstrand=false), but in this case, it's not clear to me that they should be paired; rather, it seems like that might be a chimeric molecule or some kind of structural rearrangement in the genome. Or a misassembly.

[rkanwar]

Hello Brian,

Thank you for your reply! I was using the read sequence from the SAM record for BLAT, which might explain the same strand.
Is there some easy way to filter out properly mapped PE reads [which have the right orientation and are within the fragment length distribution] which are uniquely aligned ?

Rahul

[Brian Bushnell]

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.

[rkanwar]

Hello Brian,

I removed the flag "ambiguous=toss" as you had suggested. I was able to get the correct alignments for the previous PE read and many other. But there are still a few PE reads which are not aligning properly, even though I can find the correct alignment for them using BLAT. I used the following command for aligning the data:

bbmap_cmd="${bbmap} -Xms10G -Xmx31G \
           path=${bbmap_ref} \
           in=${bbduk_fq1} \
           in2=${bbduk_fq2} \
           out=stdout \
           rgid=${ID} \
           rgsm=${SM} \
           rglb=${LIB} \
           rgpl=${PL} \
           pairlen=1000 \
           threads=12 \
           slow=t \
           local=t \
           maxsites2=2000 \
           reads=10000 \
           samversion=1.3"

One example of PE read that still behaves odd is the following. I can BLAT both the sequences and can see that there is a unique alignment which explains this read fragment. The % of such reads is very small [0.06%]. I am using grhc37 fasta file that comes with the GATK bundle as the reference. I am using v34.48 of BBMap.

R0266248:292:C5C0TACXX:4:1101:11723:6183 1:N:0:CTGAAGCTTAATCTTA	101	1	40799223	0	*	=	40799223	0	GAATTATAATCCCCATAATCTCCACCTGTCAAGGGTGGAACCAGATGTAGGTTATTTCCCTTGATGTCTACAATATCACCTTTCTTATAGATTCGCATATA	<<==<=<=>>==>=<===>>>===>=>=<=;>=>=<=>==>======<======>>=====>>=<==>>==<>====<====<>===<====<========	MC:Z:101M	BD:Z:MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMLMMMMMMMMMMMMMMMMMMMMMMMMMLMMMMMMMMMMMMMMMMMM	PG:Z:MarkDuplicates	RG:Z:AS_N+AS_N_1+FCC5C0TACXX+L4+ICTGAAGCT-TAATCTTA	BI:Z:OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
R0266248:292:C5C0TACXX:4:1101:11723:6183 1:N:0:CTGAAGCTTAATCTTA	153	1	40799223	44	101M	=	40799223	0	TTTGTCATTTTGGCGAATTACTGGAAGGTGGCGTTTCTGGCCGTAAGGCTGATATTCTTCATTTTTAAACCCAATATAACCATCCTTTCCTGTAGACTCTT	========<<3===3===<333===<=====<====><==<========<======>>=>=>==><<=<=>>>==<=>=>=========><====<>====	BD:Z:LMMMMMMLLMMMMMMMMMMMMMMMMMMMMMMMMLMMMMMMMMMMMMMMMMMMMMMMMMMMMLLLMMLMMMMMMMMMMMMMMMMMMLMMMMMMMMMMMMMMM	PG:Z:MarkDuplicates	RG:Z:AS_N+AS_N_1+FCC5C0TACXX+L4+ICTGAAGCT-TAATCTTA	BI:Z:OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO	AM:i:0	NM:i:0

Rahul

[Brian Bushnell]

BBMap is a global aligner. If you use the local flag, it will trim global alignments to become local alignments, but it will not find terrible global alignments in the first place. When run at high sensitivity, this is the output:

@PG	ID:BBMap	PN:BBMap	VN:34.59	CL:java -Djava.library.path=/usr/common/jgi/utilities/bbtools/prod-v34.59/lib -ea -Xmx31g align2.BBMapPacBio build=1 overwrite=true minratio=0.40 fastareadlen=6000 ambiguous=best minscaf=100 startpad=10000 stoppad=10000 midpad=6000 -Xmx31g t=1 k=11 in=pair3.fq out=pb.sam path=/global/projectb/sandbox/gaag/bbtools/hg191 reads=4 int=t ignorequality minratio=0.1 maxindel=40 vslow ambig=all customtag idtag ow reads=1
1 	99	chr1	40799086	4	3=3X2=3I4=40I46=	=	40799223	238	GAATTATAATCCCCATAATCTCCACCTGTCAAGGGTGGAACCAGATGTAGGTTATTTCCCTTGATGTCTACAATATCACCTTTCTTATAGATTCGCATATA	<<==<=<=>>==>=<===>>>===>=>=<=;>=>=<=>==>======<======>>=====>>=<==>>==<>====<====<>===<====<========	NM:i:46	AM:i:4	NH:i:2	YI:f:54.46	X1:Z:$1,0,40807085,40807142,0,00,32,3075,4058,6517,6517,,mmmSSSmmIIImmmmIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIImmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm$6,1,154470997,154471095,0,00,33,3754,4697,0,4697,,mSSSmmIIIImmImmSmmmIIIIIIIIIIIIIIImmmSmmSmmmmDDDDDDDDDDDDDDDDDDmmmmmmmmSmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm	X2:Z:m3S3mmI3m4I40m46	X3:i:4058	X4:i:101	X5:Z:0	X6:i:2054
1 	369	chr5	154462998	5	1=3X2=4I2=1I2=1X3=15I3=1X2=1X4=18D8=1X47=	chr1	40799223	0	*	*	NM:i:45	AM:i:5	NH:i:2	YI:f:62.18	X1:Z:$1,0,40807085,40807142,0,00,32,3075,4058,6517,6517,,mmmSSSmmIIImmmmIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIImmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm$6,1,154470997,154471095,0,00,33,3754,4697,0,4697,,mSSSmmIIIImmImmSmmmIIIIIIIIIIIIIIImmmSmmSmmmmDDDDDDDDDDDDDDDDDDmmmmmmmmSmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm	X2:Z:mS3mmI4mmImmSm3I15m3SmmSm4D18m8Sm47	X3:i:4697	X4:i:101	X5:Z:0	X6:i:34823
1 	147	chr1	40799223	44	101=	=	40799086	-238	TTTGTCATTTTGGCGAATTACTGGAAGGTGGCGTTTCTGGCCGTAAGGCTGATATTCTTCATTTTTAAACCCAATATAACCATCCTTTCCTGTAGACTCTT	========<<3===3===<333===<=====<====><==<========<======>>=>=>==><<=<=>>>==<=>=>=========><====<>====	NM:i:0	AM:i:40	NH:i:1	YI:f:100	X1:Z:$1,1,40807222,40807322,0,11,81,10090,10090,11079,11079,,mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm	X2:Z:m101	X3:i:10090	X4:i:101	X5:Z:0	X6:i:6159

The two alignments for the first read are at 62% and 54% identity, both of which are below what BBMap looks for at default settings.

[rkanwar]

Thanks, that helps a lot !

[holmrenser]

Hi Brian,

I am trying to extract reads from a bbmap sam-file using picard SamToFastq and get an error about unpaired reads. Whether I tell BBmap to output paired only or not, picard SamToFastq claims there are reads flagged as paired that are missing a pair.

BBmap command:
bbmap.sh ref=chloroplast.reference.fasta in=filter.temp.filtered.R1.fastq in2=filter.temp.filtered.R2.fastq threads=10 outm=folder/prefix.sam covstats=folder/prefix.covstats.txt basecov=folder/prefix.basecov.txt -Xmx10G fast=t keepnames=t machineout=t pairedonly=t

picard SamToFastq command:
picard.jar SamToFastq INPUT=folder/prefix.sam F=folder/prefix.R1.fastq F2=folder/prefix.R2.fastq FU=folder/prefix.unpaired.fastq

Error:
Exception in thread "main" htsjdk.samtools.SAMFormatException: SAM validation error: ERROR: Found 29230 unpaired mates

When I run it without the pairedonly=t option I get the following error:
Exception in thread "main" htsjdk.samtools.SAMFormatException: SAM validation error: ERROR: Found 61616 unpaired mates

Any suggestion?

Thanks in advance

[Brian Bushnell]

That's odd. Your command looks fine. Unfortunately I will be gone until Monday and won't be able to look into this until then, but would you be able to email me the fastqs and reference? Or, ideally, a small subset of the reads that triggers the same problem.

Thanks,
Brian

[susanklein]

Hi Brian,

is there any way to use bbmap to map very long reads from e.g. nanopore (minION).

Also I read on another thread that dedup.sh can be used to merge assemblies. Can it produce a consensus sequence?

Thanks very much.

S

[GenoMax]

BBMap can map reads of up to 6kb: http://seqanswers.com/forums/showpos...43&postcount=2