[minja]

Hello Brain,

I'm trying to make use of bbmapskimmer.
Here is a simple fasta example (read.fasta):

>1
GAACATGATCCCCTTGTACTA[B]C[/B]TACATTATCACTAGCTTTATGTTTTCTA
>2
TAGAAAACATAAAGCTAGTGATAATGTAGTAGTACAAGGGGATCATGTTC
>3
GAACATGATCCCCTTGTACTA[B]T[/B]TACATTATCACTAGCTTTATGTTTTCTA

Read 1 and read 3 are only different in single "T/C" letter, read 2 is reverse-complement of 1.

I'm aligning to mm10 genome:

~/Distr/bbmap/bbmapskimmer.sh in=~/tmp/HiC_repeats/input/temp/read.fasta out=test_bbmap1.sam ref=~/tmp/HiC_repeats/fasta/mm10.fa ambiguous=all vslow k=8 minid=0.6 maxsites2=8000 maxsites=8000 local

Here is test_bbmap1.sam:

@HD    VN:1.4    SO:unsorted
@SQ    SN:chr10    LN:130694993
....
@SQ    SN:chrM    LN:16299
@SQ    SN:chrX    LN:171031299
@SQ    SN:chrY    LN:91744698
@PG    ID:BBMap    PN:BBMap    VN:37.02    CL:java -Djava.library.path=/mnt/storage/home/vsfishman/Distr/bbmap/jni/ -ea -Xmx98546m align2.BBMapPacBioSkimmer build=1 overwrite=true minratio=0.40 fastareadlen=6000 ambig=all minscaf=100 startpad=10000 stoppad=10000 midpad=6000 in=/mnt/storage/home/vsfishman/tmp/HiC_repeats/input/temp/read.fasta out=test_bbmap1.sam ref=/mnt/storage/home/vsfishman/tmp/HiC_repeats/fasta/mm10.fa ambiguous=all vslow k=8 minid=0.6 maxsites2=8000 maxsites=8000 local
1    0    chr16    39205299    43    50=    *    0    0    GAACATGATCCCCTTGTACTACTACATTATCACTAGCTTTATGTTTTCTA    *    NM:i:0    AM:i:43    NH:i:2
1    272    chr10    38707112    39    28=1X21=    *    0    0    *    *    NM:i:1    AM:i:39    NH:i:2
2    16    chr16    39205299    43    50=    *    0    0    GAACATGATCCCCTTGTACTACTACATTATCACTAGCTTTATGTTTTCTA    *    NM:i:0    AM:i:43    NH:i:2
2    256    chr10    38707112    39    28=1X21=    *    0    0    *    *    NM:i:1    AM:i:39    NH:i:2
3    16    chr10    38707112    3    50=    *    0    0    TAGAAAACATAAAGCTAGTGATAATGTAATAGTACAAGGGGATCATGTTC    *    XT:A:R    NM:i:0    AM:i:3    NH:i:80
3    272    chr10    72154315    3    50=    *    0    0    *    *    NM:i:0    AM:i:3    NH:i:80
3    256    chr10    77806115    3    50=    *    0    0    *    *    NM:i:0    AM:i:3    NH:i:80
3    272    chr10    127797927    3    50=    *    0    0    *    *    NM:i:0    AM:i:3    NH:i:80
3    256    chr13    4007635    3    50=    *    0    0    *    *    NM:i:0    AM:i:3    NH:i:80
3    256    chr14    17647057    3    50=    *    0    0    *    *    NM:i:0    AM:i:3    NH:i:80
3    256    chr14    40019581    3    50=    *    0    0    *    *    NM:i:0    AM:i:3    NH:i:80
3    256    chr15    50069753    3    50=    *    0    0    *    *    NM:i:0    AM:i:3    NH:i:80
3    272    chr15    72113945    3    50=    *    0    0    *    *    NM:i:0    AM:i:3    NH:i:80
3    256    chr15    75259455    3    50=    *    0    0    *    *    NM:i:0    AM:i:3    NH:i:80
3    256    chr17    69579150    3    50=    *    0    0    *    *    NM:i:0    AM:i:3    NH:i:80
3    256    chr19    60830365    3    50=    *    0    0    *    *    NM:i:0    AM:i:3    NH:i:80
3    256    chr1    95272134    3    50=    *    0    0    *    *    NM:i:0    AM:i:3    NH:i:80
3    256    chr1    102465771    3    50=    *    0    0    *    *    NM:i:0    AM:i:3    NH:i:80
....

(I reduced some header lines and some alinments of read 3)

And, just to note, without "local" option I got same results.

As you can see, there are many perfect alignmens for read3, but only 2 alignments of reads 1 and 2. As I understand, with minid=0.6 all alignments of read3 should be also reported for reads 1 and 2?

Could you, please, help to clarify this?
And one more question, why reads 1 and 2 do not have "XT:A:" tag?

Best,
Minja

[GenoMax]

While you wait for @Brian to respond (he is on California Daylight Time) you could turn off "pigz=f unpigz=f" and see if that allows the job to move forward.

[sebastian_1]

Hi Brian,

Thank you for your quick reply. I've tried to run BBNorm with prefilter=t and bits=16 on our server, but I still run out of memory. I have a server with 128GB of RAM, and I've tried also on a 256 GB but it didn't work. But I saw that the high CPU usage (more than defined in threads=x) is just at the beginning, for a short time and afterwards java limits itself to the required number of cores.
Luckily, I've found that the PBS system divides the CPU hours to the walltime, and thus gets the CPU usage. The workaround was to ask the job as a STDIN, and be idle for around 30 minutes, which will create a "buffer", and thus even if I start with much higher initial CPU usage, the time buffer is enough to compensate, and thus the job doesn't get killed by the scheduler.

However now I get an error with pigz:

Table creation time: 1667.979 seconds.
Started output threads.
pigz: write error code 122
pigz: abort: write error on <stdout>
Exception in thread "Thread-53" java.lang.RuntimeException: java.io.IOException: Stream closed
at stream.ReadStreamByteWriter.run(ReadStreamByteWriter.java:31)
Caused by: java.io.IOException: Stream closed
at java.lang.ProcessBuilder$NullOutputStream.write(ProcessBuilder.java:433)
at java.io.OutputStream.write(OutputStream.java:116)
at java.io.BufferedOutputStream.write(BufferedOutputStream.java:122)
at stream.ReadStreamByteWriter.writeFastq(ReadStreamByteWriter.java:451)
at stream.ReadStreamByteWriter.processJobs(ReadStreamByteWriter.java:96)
at stream.ReadStreamByteWriter.run2(ReadStreamByteWriter.java:41)
at stream.ReadStreamByteWriter.run(ReadStreamByteWriter.java:27)
Exception in thread "Thread-64" java.lang.RuntimeException: Writing to a terminated thread.
at stream.ConcurrentGenericReadOutputStream.write(ConcurrentGenericReadOutputStream.java:207)
at stream.ConcurrentGenericReadOutputStream.addOrdered(ConcurrentGenericReadOutputStream.java:193)
at stream.ConcurrentGenericReadOutputStream.add(ConcurrentGenericReadOutputStream.java:98)
at jgi.KmerNormalize$ProcessThread.normalizeInThread(KmerNormalize.java:3129)
at jgi.KmerNormalize$ProcessThread.run(KmerNormalize.java:2801)
Output buffer became full; key 95560 waiting on 95304.

LE: I set the $TMPDIR on the cluster to my current working directory and i got rid of the pigz error:-)

[minja]

bbmapskimmer does not find all results

Hello Brain,

I'm trying to make use of bbmapskimmer.
Here is a simple fasta example (read.fasta):

>1
GAACATGATCCCCTTGTACTA[B]C[/B]TACATTATCACTAGCTTTATGTTTTCTA
>2
TAGAAAACATAAAGCTAGTGATAATGTAGTAGTACAAGGGGATCATGTTC
>3
GAACATGATCCCCTTGTACTA[B]T[/B]TACATTATCACTAGCTTTATGTTTTCTA

Read 1 and read 3 are only different in single "T/C" letter, read 2 is reverse-complement of 1.

I'm aligning to mm10 genome:

~/Distr/bbmap/bbmapskimmer.sh in=~/tmp/HiC_repeats/input/temp/read.fasta out=test_bbmap1.sam ref=~/tmp/HiC_repeats/fasta/mm10.fa ambiguous=all vslow k=8 minid=0.6 maxsites2=8000 maxsites=8000 local

Here is test_bbmap1.sam:

@HD	VN:1.4	SO:unsorted
@SQ	SN:chr10	LN:130694993
....
@SQ	SN:chrM	LN:16299
@SQ	SN:chrX	LN:171031299
@SQ	SN:chrY	LN:91744698
@PG	ID:BBMap	PN:BBMap	VN:37.02	CL:java -Djava.library.path=/mnt/storage/home/vsfishman/Distr/bbmap/jni/ -ea -Xmx98546m align2.BBMapPacBioSkimmer build=1 overwrite=true minratio=0.40 fastareadlen=6000 ambig=all minscaf=100 startpad=10000 stoppad=10000 midpad=6000 in=/mnt/storage/home/vsfishman/tmp/HiC_repeats/input/temp/read.fasta out=test_bbmap1.sam ref=/mnt/storage/home/vsfishman/tmp/HiC_repeats/fasta/mm10.fa ambiguous=all vslow k=8 minid=0.6 maxsites2=8000 maxsites=8000 local
1	0	chr16	39205299	43	50=	*	0	0	GAACATGATCCCCTTGTACTACTACATTATCACTAGCTTTATGTTTTCTA	*	NM:i:0	AM:i:43	NH:i:2
1	272	chr10	38707112	39	28=1X21=	*	0	0	*	*	NM:i:1	AM:i:39	NH:i:2
2	16	chr16	39205299	43	50=	*	0	0	GAACATGATCCCCTTGTACTACTACATTATCACTAGCTTTATGTTTTCTA	*	NM:i:0	AM:i:43	NH:i:2
2	256	chr10	38707112	39	28=1X21=	*	0	0	*	*	NM:i:1	AM:i:39	NH:i:2
3	16	chr10	38707112	3	50=	*	0	0	TAGAAAACATAAAGCTAGTGATAATGTAATAGTACAAGGGGATCATGTTC	*	XT:A:R	NM:i:0	AM:i:3	NH:i:80
3	272	chr10	72154315	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr10	77806115	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	272	chr10	127797927	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr13	4007635	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr14	17647057	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr14	40019581	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr15	50069753	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	272	chr15	72113945	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr15	75259455	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr17	69579150	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr19	60830365	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr1	95272134	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr1	102465771	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
....

(I reduced some header lines and some alinments of read 3)

As you can see, there are many perfect alignmens for read3, but only 2 alignments of reads 1 and 2. As I understand, with minid=0.6 all alignments of read3 should be also reported for reads 1 and 2?

Could you, please, help to clarify this?

Best,
Minja

[minja]

bbmapskimmer does not find all results

Hello Brain,

I'm trying to make use of bbmapskimmer.
Here is a simple fasta example (read.fasta):

>1
GAACATGATCCCCTTGTACTA[B]C[/B]TACATTATCACTAGCTTTATGTTTTCTA
>2
TAGAAAACATAAAGCTAGTGATAATGTAGTAGTACAAGGGGATCATGTTC
>3
GAACATGATCCCCTTGTACTA[B]T[/B]TACATTATCACTAGCTTTATGTTTTCTA

Read 1 and read 3 are only different in single "T/C" letter, read 2 is reverse-complement of 1.

I'm aligning to mm10 genome:

~/Distr/bbmap/bbmapskimmer.sh in=~/tmp/HiC_repeats/input/temp/read.fasta out=test_bbmap1.sam ref=~/tmp/HiC_repeats/fasta/mm10.fa ambiguous=all vslow k=8 minid=0.6 maxsites2=8000 maxsites=8000 local

Here is test_bbmap1.sam:

@HD	VN:1.4	SO:unsorted
@SQ	SN:chr10	LN:130694993
....
@SQ	SN:chrM	LN:16299
@SQ	SN:chrX	LN:171031299
@SQ	SN:chrY	LN:91744698
@PG	ID:BBMap	PN:BBMap	VN:37.02	CL:java -Djava.library.path=/mnt/storage/home/vsfishman/Distr/bbmap/jni/ -ea -Xmx98546m align2.BBMapPacBioSkimmer build=1 overwrite=true minratio=0.40 fastareadlen=6000 ambig=all minscaf=100 startpad=10000 stoppad=10000 midpad=6000 in=/mnt/storage/home/vsfishman/tmp/HiC_repeats/input/temp/read.fasta out=test_bbmap1.sam ref=/mnt/storage/home/vsfishman/tmp/HiC_repeats/fasta/mm10.fa ambiguous=all vslow k=8 minid=0.6 maxsites2=8000 maxsites=8000 local
1	0	chr16	39205299	43	50=	*	0	0	GAACATGATCCCCTTGTACTACTACATTATCACTAGCTTTATGTTTTCTA	*	NM:i:0	AM:i:43	NH:i:2
1	272	chr10	38707112	39	28=1X21=	*	0	0	*	*	NM:i:1	AM:i:39	NH:i:2
2	16	chr16	39205299	43	50=	*	0	0	GAACATGATCCCCTTGTACTACTACATTATCACTAGCTTTATGTTTTCTA	*	NM:i:0	AM:i:43	NH:i:2
2	256	chr10	38707112	39	28=1X21=	*	0	0	*	*	NM:i:1	AM:i:39	NH:i:2
3	16	chr10	38707112	3	50=	*	0	0	TAGAAAACATAAAGCTAGTGATAATGTAATAGTACAAGGGGATCATGTTC	*	XT:A:R	NM:i:0	AM:i:3	NH:i:80
3	272	chr10	72154315	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr10	77806115	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	272	chr10	127797927	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr13	4007635	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr14	17647057	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr14	40019581	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr15	50069753	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	272	chr15	72113945	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr15	75259455	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr17	69579150	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr19	60830365	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr1	95272134	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
3	256	chr1	102465771	3	50=	*	0	0	*	*	NM:i:0	AM:i:3	NH:i:80
....

(I reduced some header lines and some alinments of read 3)

As you can see, there are many perfect alignmens for read3, but only 2 alignments of reads 1 and 2. As I understand, with minid=0.6 all alignments of read3 should be also reported for reads 1 and 2?

Could you, please, help to clarify this?

Best,
Minja

[Brian Bushnell]

Hi Sebastian,

You can reduce the memory needs of BBNorm by adding the flags "prefilter" and "bits=16", which might allow it to run on a lower-memory computer, depending on the complexity of the input data.

But as for reducing the CPU utilization... that's tricky. You can add "pigz=f unpigz=f" to disable pigz, so you will only be left with java threads. Java thread scheduling is nondeterministic due to garbage collection and various other background stuff, and I'm not sure if it's possible to ensure the process will be capped at a certain number of total threads without binding the process to specific cores. But you can try manually specifying the number of GC threads:

-XX:ParallelGCThreads=2

This command may or may not work, depending on your java version. You'd have to skip (or modify) the shellscript so that the actual java command would look something like this:

java -ea -Xmx200g -Xms200g -XX:ParallelGCThreads=2 -cp /path/to/bbmap/current/ jgi.KmerNormalize in1=st1c_R1.fastq.gz in2=st1c_R2.fastq.gz out1=st1c_R1_norm.fastq.gz out2=st1c_R2_norm.fastq.gz threads=12 target=60 mindepth=2

[sebastian_1]

Hello,

I am new to this forum so I apologize if I posted in the wrong place.

Currently I am trying to assemble a single cell genome from some Illumina Hiseq 2x250 bp reads.
One of my first steps is to do normalization of the data. Previously I was using bbnorm.sh because it's fast and very good. However, with this data I cannot use bbnorm.sh on our server due to memory issues. (I have around 100 million PE reads).
So I tried to use bbnorm on a cluster where I allocated 16 cores and 700GB of RAM.

I run the software using this command line:
bbnorm.sh in1=st1c_R1.fastq.gz in2=st1c_R2.fastq.gz out1=st1c_R1_norm.fastq.gz out2=st1c_R2_norm.fastq.gz threads=12 target=60 mindepth=2

The software runs fine, but the process will be killed by the scheduler because of high CPU usage.

=>> PBS: job killed: ncpus 59.52 exceeded limit 16 (sum)

I tried to increase the number of threads requested to 30, lower the number of threads in the command line to 10, but still I have the same problem.

The cluster has installed BBMap version 36.92. It looks like for some reason the number of threads which I define is not taken into account, even if during the run it tells me that the number of threads was set to 12.

Settings:
threads: 12
k: 31
deterministic: true
toss error reads: false
passes: 1
bits per cell: 16
cells: 145.37B
hashes: 3
base min quality: 5
kmer min prob: 0.5

target depth: 240
min depth: 2
max depth: 300
min good kmers: 15
depth percentile: 64.8
ignore dupe kmers: true
fix spikes: false

Any suggestions? I really cannot ask for the entire node just for me.
I understand that due to other components like pigz, it can use more than 12 cores, but here it looks like uses 60 cores, or all available cores.

Thank you,
Sebastian

[shimingt]

Dear Brian,

Thanks for your reply.

This is from one of the mapping files (.csv) for mmgenome.

tanshiming@S620100019205:~/Downloads$ head C13.11.14.csv
"Name","Consensus length","Total read count","Single reads","Reads in pairs","Average coverage","Reference sequence","Reference length","Reference common name","Reference Latin name"
"1 mapping","5621","103","19","84","1.443","1","8264","",""
"2 mapping","370","6","2","4","0.625","2","1027","",""
"3 mapping","1665","198","52","146","13.536","3","1665","",""
"4 mapping","94","1","1","0","0.01","4","9056","",""
"5 mapping","3042","95","19","76","3.205","5","3343","",""
"6 mapping","901","8","2","6","9.663E-3","6","98207","",""
"7 mapping","5788","147","11","136","2.612","7","6480","",""
"8 mapping","14602","381","47","334","2.782","8","15790","",""
"9 mapping","1403","1056","282","774","85.135","9","1403","",""
tanshiming@S620100019205:~/Downloads$

[Brian Bushnell]

BBMap can directly output coverage using the "covstats" or "basecov" flags, or it can be generated from sam files with pileup.sh. Can you post a few lines from one of these mmgenome csv files so I can see what the format looks like?

[shimingt]

Converting a coverage file into a csv format for mmgenome

Hello there,

Is there a way I could convert the mapping files that bbmap produced into csv coverage file that can be used for the mmgenome software (link: http://madsalbertsen.github.io/mmgen...ad_data.html)?

Thanks!

[Marmottontaine]

Thank you very much GenoMax and Brian for your helpful answers.

[Brian Bushnell]

And for the parameters : maxindel, qtrim and untrim, they are both set as the default then as I asked in the same command line. Will it be problematic?

Later parameters override earlier parameters, so that's not a problem.

As for 0.65 versus 0.56, that looks like either a typo, or else I changed the default at some point but forgot to update the help info; I'll fix that. Thanks! In the meantime you can always manually add "minratio=0.65" to override it.

Edit: To clarify, 0.56 is the intended behavior, and the description is in error.

[GenoMax]

I see. @Brian hopefully will be along later today with an official answer but my guess is that it probably has to do with ambigous2=toss option that you have selected.

[Marmottontaine]

Thank you GenoMax for your answer.
I let minratio at the default value. So it should be 0.65 as said in bbsplit help. here it's 0.56.
And for the parameters : maxindel, qtrim and untrim, they are both set as the default then as I asked in the same command line. Will it be problematic?
Many thanks.

[GenoMax]

Which exact parameters are you referring to?