[GenoMax]

Since you asked for 30K reads you got 30K paired-end reads.

I would sugegst trying without the metagenome=t option for a single small input fasta. If you do have the entire genome then use that in your randomread generation step along with the plasmid.

[bartn]

randomreads.sh

Hi!

I want to created a (small) RNAseq test set from a small bacterial plasmid. So I used :

randomreads.sh reads=30000 paired=t metagenome=t  ref=transcripts.fasta out1=test_RNA_1.fq out2=test_RNA_2.fq

(Then I get 30k reads for each pair. No problem, I guess that's because the default is single end)

Then I mapped back the reads with BBmap

bbmap.sh ref=Rhizobium_Leguminosarum_plasmid_NZ_CP025014.1_transcripts.fasta  in1=Unlock_test_RNA_1.fq in2=Unlock_test_RNA_2.fq covstats=RNA.cov

Results:

Reads:                               	60000
Mapped reads:                        	54225
Mapped bases:                        	8063372
Ref scaffolds:                       	562
Ref bases:                           	521537

Percent mapped:                      	90.375
Percent proper pairs:                	64.822
Average coverage:                    	15.461
Standard deviation:                    	51.768
Percent scaffolds with any coverage: 	59.96
Percent of reference bases covered:  	45.87

The mapping percentages seem very low to me. Mapping to the genome gives similar results.
I am using randomreads.sh in the wrong way for this?

[bartn]

Hi!

I want to create a (small) bacterial RNAseq testset from transcripts using randomreads.sh:

randomreads.sh reads=30000 ref=transcripts.fasta metagenome=t out1=test_RNA_1.fq out2=test_RNA_2.fq paired=t

Mapping it back to the transcripts:

bbmap.sh ref=transcripts.fasta  in1=test_RNA_1.fq in2=test_RNA_2.fq covstats=RNA.cov

Results

Reads:                               	60000
Mapped reads:                        	54225
Mapped bases:                        	8063372
Ref scaffolds:                       	562
Ref bases:                           	521537

Percent mapped:                      	90.375
Percent proper pairs:                	64.822
Average coverage:                    	15.461
Standard deviation:                    	51.768
Percent scaffolds with any coverage: 	59.96
Percent of reference bases covered:  	45.87

Mapping it to the genome:

Reads:                               	60000
Mapped reads:                        	54104
Mapped bases:                        	8061481
Ref scaffolds:                       	1
Ref bases:                           	592529

Percent mapped:                      	90.173
Percent proper pairs:                	64.573
Average coverage:                    	13.605
Standard deviation:                    	48.628
Percent scaffolds with any coverage: 	100.00
Percent of reference bases covered:  	41.02

Why are the mapping percentages so low, what I missing with either generating the reads or mapping them back?

(bowtie2 gives an even worse result..
15482 (51.61%) aligned concordantly 0 times
14351 (47.84%) aligned concordantly exactly 1 time
167 (0.56%) aligned concordantly >1 times
)

[GenoMax]

@Noobie: I don't think so. You will have to do them independently. What is the use case here, if I may ask?

[noobie]

I want to subsample 10 million reads from both an interleaved reads file and a singletons reads file at the same time. Can reformat.sh do this?

[bwubb]

Have you tried adding "VALIDATION_STRINGENCY=LENIENT" to your picard command? Default value is STRICT and sometime it leads to picard throwing errors like this.
...

Yes indeed, but it will still be an ERROR.

However, after more tinkering I was able to come up with a combination that worked.

reformat.sh -Xmx5g in=qSortA.bam out=qSortA_reads.fq
reformat.sh -Xmx5g int=t addcolon=t uniquenames=t in=qSortA_reads.fq out=out_reads.fq
repair.sh -Xmx5g in=out_reads.fq out1=reads1.fq out2=reads2.fq

Taking those reads into alignment and so forth will pass ValidateSamFile. Im guessing these actions are allowing me to completely rewrite the pairing status, which I am perfectly ok with.

Thank you for the comments and help!

[GenoMax]

Have you tried adding "VALIDATION_STRINGENCY=LENIENT" to your picard command? Default value is STRICT and sometime it leads to picard throwing errors like this.

It does looks like reformat.sh took out some of the reads so the problem reads have been hopefully addressed.

[bwubb]

@bwbubb: Option for reformat.sh that you are looking for is "pairedonly=f Toss reads that are not mapped as proper pairs". Set pairedonly=t and process your BAM files.

EDIT:

Thank you @GenoMax. This actually does not work for for this case.

[bwubb@node061 disambres]$ reformat.sh pairedonly=t in=input.bam out=fix.bam
java -ea -Xmx200m -cp /home/bwubb/software/bbmap/current/ jgi.ReformatReads pairedonly=t in=input.bam out=fix.bam
Executing jgi.ReformatReads [pairedonly=t, in=input.bam, out=fix.bam]

Found samtools 1.9
Input is being processed as unpaired
Input:                          13726088 reads                  2041390865 bases
Output:                         12583399 reads (91.68%)         1880228501 bases (92.11%)

Time:                           349.142 seconds.
Reads Processed:      13726k    39.31k reads/sec
Bases Processed:       2041m    5.85m bases/sec

[bwubb@node061 disambres]$ samtools index fix.bam
[bwubb@node061 disambres]$ java -jar ~/software/picard/2.20.2/picard.jar ValidateSamFile I=fix.bam MODE=SUMMARY
INFO    2019-06-13 08:00:16     ValidateSamFile

...

## HISTOGRAM    java.lang.String
Error Type      Count
ERROR:MATE_NOT_FOUND    2

[Thu Jun 13 08:02:01 EDT 2019] picard.sam.ValidateSamFile done. Elapsed time: 1.74 minutes.
Runtime.totalMemory()=2075918336
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp

Not sure how things work internally, but from my understanding the "problem" reads think they are paired.


I toiled with this more yesterday and I believe I can use reformat.sh to make an interleaved fastq file, and then use repair.sh to make read1.fq and read2.fq and take it back to alignment. All attempts at correcting the bam directly seem to fail. I have yet to have total success, which in my case would be passing picardtools ValidateSamFile.

I am aligning with bwa mem -M, having some other issue in the validation that maybe is unrelated or will be fixed with the right reformat/repair arguments.

[GenoMax]

@bwbubb: Option for reformat.sh that you are looking for is "pairedonly=f Toss reads that are not mapped as proper pairs". Set pairedonly=t and process your BAM files.

[bwubb]

Greetings,

I have a seemingly simple problem, that Im sure bbmap can help with, but Im unclear what steps must be taken for success.

I have several bams from paired end sequencing where some number of mates have been removed. The remaining read still says it is paired so every tool seems to throw an error. Samtools or picard fixmate does not correct it, at best they yell the read names in question.

I am trying to use reformat.sh to remove those mates that remain, but I also seem to get a java exception, even when just trying to convert to fastq.

[bwubb@node107 disambres]$ reformat.sh -Xmx5g in=qsortA.bam out1=reads1.fq out2=reads2.fq ow
java -ea -Xmx5g -cp /home/bwubb/software/bbmap/current/ jgi.ReformatReads -Xmx5g in=qsortA.bam out1=reads1.fq out2=reads2.fq ow
Executing jgi.ReformatReads [-Xmx5g, in=qsortA.bam, out1=reads1.fq, out2=reads2.fq, ow]

Set INTERLEAVED to true
Found samtools 1.9
Input is being processed as paired
Exception in thread "Thread-2" java.lang.AssertionError: K00315:202:H577WBBXY:7:1101:1438:38767 13      -1      -       57328513        57328662        1000000100000000011     1       0  60       TCCAGAAAGAGAGTTCTTGAAAGAAAGAGGCGCTATCATTTTGACACAGATGGNAAGGGCTCGATTCACGATCAAAAAGGCTCCAAAAANAAAAAAAATTCTGTTTGTTNTNCTTTTTCCTCCNGATTCTAGTTTTTTTNTTATNTTTTG  AAFFFFFAJJJJJ7FFFJJ7JJ7FJJJJJFFJJJJJFJJJJJF-A-AAFFJFA!AF---AAJJ-AJF<FF-AAFAF<AAAJJJJFA-7A!--77-<A<JJ<A--7A<FJ!J!<FJJJJA<FF7!)7AJF---7AFJJFJ!JF--!7-7<F      .       30      C68m28Nm53      .

K00315:202:H577WBBXY:7:1101:1438:38767  113     19      57328582        60      68S82M  1       176081441       0       CAAAANATAANAAAAAAACTAGAATCNGGAGGAAAAAGNANAACAAACAGAATTTTTTTTNTTTTTGGAGCCTTTTTGATCGTGAATCGAGCCCTTNCCATCTGTGTCAAAATGATAGCGCCTCTTTCTTTCAAGAACTCTCTTTCTGGA      F<7-7!--FJ!JFJJFA7---FJA7)!7FF<AJJJJF<!J!JF<A7--A<JJ<A<-77--!A7-AFJJJJAAA<FAFAA-FF<FJA-JJAA---FA!AFJFFAA-A-FJJJJJFJJJJJFFJJJJJF7JJ7JJFFF7JJJJJAFFFFFAA      NM:i:1  MD:Z:28C53      MC:Z:150M       AS:i:80 XS:i:21 SA:Z:1,176081357,+,81S62M7S,54,7;       RG:Z:FGC2033.7
flag=1110001

        at stream.SamReadInputStream.toReadList(SamReadInputStream.java:136)
        at stream.SamReadInputStream.fillBuffer(SamReadInputStream.java:89)
        at stream.SamReadInputStream.hasMore(SamReadInputStream.java:53)
        at stream.ConcurrentGenericReadInputStream$ReadThread.readLists(ConcurrentGenericReadInputStream.java:643)
        at stream.ConcurrentGenericReadInputStream$ReadThread.run(ConcurrentGenericReadInputStream.java:635)

I thought I could make read1.fq reads2.fq and use repair.sh Clearly it finds the two problem reads. I have tried other arguments with the reformat.sh, and a few other tools, but none want to correct the MATE_NOT_FOUND problem. It would seem impossible that I am the first person to run into this issue. Thank you very much for any advice/help.

-bwubb

EDIT: Sorry if this post shows up multiple times. I literally tried making it three times and nothing showed up for a full day. Tried a quick reply and it showed up immediately.

[bwubb]

Greetings,

I have a seemingly simple problem, that Im sure bbmap can help with, but Im unclear what steps must be taken for success.

I have several bams from paired end sequencing where some number of mates have been removed. The remaining read still says it is paired so every tool seems to throw an error. Samtools or picard fixmate does not correct it, at best they yell the read names in question.

I am trying to use reformat.sh to remove those mates that remain, but I also seem to get a java exception, even when just trying to convert to fastq.

[bwubb@node107 disambres]$ reformat.sh -Xmx5g in=qsortA.bam out1=reads1.fq out2=reads2.fq ow
java -ea -Xmx5g -cp /home/bwubb/software/bbmap/current/ jgi.ReformatReads -Xmx5g in=qsortA.bam out1=reads1.fq out2=reads2.fq ow
Executing jgi.ReformatReads [-Xmx5g, in=qsortA.bam, out1=reads1.fq, out2=reads2.fq, ow]

Set INTERLEAVED to true
Found samtools 1.9
Input is being processed as paired
Exception in thread "Thread-2" java.lang.AssertionError: K00315:202:H577WBBXY:7:1101:1438:38767 13      -1      -       57328513        57328662        1000000100000000011     1       0  60       TCCAGAAAGAGAGTTCTTGAAAGAAAGAGGCGCTATCATTTTGACACAGATGGNAAGGGCTCGATTCACGATCAAAAAGGCTCCAAAAANAAAAAAAATTCTGTTTGTTNTNCTTTTTCCTCCNGATTCTAGTTTTTTTNTTATNTTTTG  AAFFFFFAJJJJJ7FFFJJ7JJ7FJJJJJFFJJJJJFJJJJJF-A-AAFFJFA!AF---AAJJ-AJF<FF-AAFAF<AAAJJJJFA-7A!--77-<A<JJ<A--7A<FJ!J!<FJJJJA<FF7!)7AJF---7AFJJFJ!JF--!7-7<F      .       30      C68m28Nm53      .

K00315:202:H577WBBXY:7:1101:1438:38767  113     19      57328582        60      68S82M  1       176081441       0       CAAAANATAANAAAAAAACTAGAATCNGGAGGAAAAAGNANAACAAACAGAATTTTTTTTNTTTTTGGAGCCTTTTTGATCGTGAATCGAGCCCTTNCCATCTGTGTCAAAATGATAGCGCCTCTTTCTTTCAAGAACTCTCTTTCTGGA      F<7-7!--FJ!JFJJFA7---FJA7)!7FF<AJJJJF<!J!JF<A7--A<JJ<A<-77--!A7-AFJJJJAAA<FAFAA-FF<FJA-JJAA---FA!AFJFFAA-A-FJJJJJFJJJJJFFJJJJJF7JJ7JJFFF7JJJJJAFFFFFAA      NM:i:1  MD:Z:28C53      MC:Z:150M       AS:i:80 XS:i:21 SA:Z:1,176081357,+,81S62M7S,54,7;       RG:Z:FGC2033.7
flag=1110001

        at stream.SamReadInputStream.toReadList(SamReadInputStream.java:136)
        at stream.SamReadInputStream.fillBuffer(SamReadInputStream.java:89)
        at stream.SamReadInputStream.hasMore(SamReadInputStream.java:53)
        at stream.ConcurrentGenericReadInputStream$ReadThread.readLists(ConcurrentGenericReadInputStream.java:643)
        at stream.ConcurrentGenericReadInputStream$ReadThread.run(ConcurrentGenericReadInputStream.java:635)

I thought I could make read1.fq reads2.fq and use repair.sh Clearly it finds the two problem reads. I have tried other arguments with the reformat.sh, and a few other tools, but none want to correct the MATE_NOT_FOUND problem. It would seem impossible that I am the first person to run into this issue. Thank you very much for any advice/help.

-bwubb

[bwubb]

Greetings,

I have a seemingly simple problem, that Im sure bbmap can help with, but Im unclear what steps must be taken for success.

I have several bams from paired end sequencing where some number of mates have been removed. The remaining read still says it is paired so every tool seems to throw an error. Samtools or picard fixmate does not correct it, at best they yell the read names in question.

I am trying to use reformat.sh to remove those mates that remain, but I also seem to get a java exception, even when just trying to convert to fastq.

[bwubb@node107 disambres]$ reformat.sh -Xmx5g in=qsortA.bam out1=reads1.fq out2=reads2.fq ow
java -ea -Xmx5g -cp /home/bwubb/software/bbmap/current/ jgi.ReformatReads -Xmx5g in=qsortA.bam out1=reads1.fq out2=reads2.fq ow
Executing jgi.ReformatReads [-Xmx5g, in=qsortA.bam, out1=reads1.fq, out2=reads2.fq, ow]

Set INTERLEAVED to true
Found samtools 1.9
Input is being processed as paired
Exception in thread "Thread-2" java.lang.AssertionError: K00315:202:H577WBBXY:7:1101:1438:38767 13      -1      -       57328513        57328662        1000000100000000011     1       0  60       TCCAGAAAGAGAGTTCTTGAAAGAAAGAGGCGCTATCATTTTGACACAGATGGNAAGGGCTCGATTCACGATCAAAAAGGCTCCAAAAANAAAAAAAATTCTGTTTGTTNTNCTTTTTCCTCCNGATTCTAGTTTTTTTNTTATNTTTTG  AAFFFFFAJJJJJ7FFFJJ7JJ7FJJJJJFFJJJJJFJJJJJF-A-AAFFJFA!AF---AAJJ-AJF<FF-AAFAF<AAAJJJJFA-7A!--77-<A<JJ<A--7A<FJ!J!<FJJJJA<FF7!)7AJF---7AFJJFJ!JF--!7-7<F      .       30      C68m28Nm53      .

K00315:202:H577WBBXY:7:1101:1438:38767  113     19      57328582        60      68S82M  1       176081441       0       CAAAANATAANAAAAAAACTAGAATCNGGAGGAAAAAGNANAACAAACAGAATTTTTTTTNTTTTTGGAGCCTTTTTGATCGTGAATCGAGCCCTTNCCATCTGTGTCAAAATGATAGCGCCTCTTTCTTTCAAGAACTCTCTTTCTGGA      F<7-7!--FJ!JFJJFA7---FJA7)!7FF<AJJJJF<!J!JF<A7--A<JJ<A<-77--!A7-AFJJJJAAA<FAFAA-FF<FJA-JJAA---FA!AFJFFAA-A-FJJJJJFJJJJJFFJJJJJF7JJ7JJFFF7JJJJJAFFFFFAA      NM:i:1  MD:Z:28C53      MC:Z:150M       AS:i:80 XS:i:21 SA:Z:1,176081357,+,81S62M7S,54,7;       RG:Z:FGC2033.7
flag=1110001

        at stream.SamReadInputStream.toReadList(SamReadInputStream.java:136)
        at stream.SamReadInputStream.fillBuffer(SamReadInputStream.java:89)
        at stream.SamReadInputStream.hasMore(SamReadInputStream.java:53)
        at stream.ConcurrentGenericReadInputStream$ReadThread.readLists(ConcurrentGenericReadInputStream.java:643)
        at stream.ConcurrentGenericReadInputStream$ReadThread.run(ConcurrentGenericReadInputStream.java:635)

I thought I could make read1.fq reads2.fq and use repair.sh Clearly it finds the two problem reads. I have tried other arguments with the reformat.sh, and a few other tools, but none want to correct the MATE_NOT_FOUND problem. It would seem impossible that I am the first person to run into this issue. Thank you very much for any advice/help.

-bwubb

[bwubb]

Fixing Mate not found error

Greetings,

I have a seemingly simple problem, that Im sure bbmap can help with, but Im unclear what steps must be taken for success.

I have several bams from paired end sequencing where some number of mates have been removed. The remaining read still says it is paired so every tool seems to throw an error. Samtools or picard fixmate does not correct it, at best they yell the read names in question.

I am trying to use reformat.sh to remove those mates that remain, but I also seem to get a java exception, even when just trying to convert to fastq.

[bwubb@node107 disambres]$ reformat.sh -Xmx5g in=qsortA.bam out1=reads1.fq out2=reads2.fq ow
java -ea -Xmx5g -cp /home/bwubb/software/bbmap/current/ jgi.ReformatReads -Xmx5g in=qsortA.bam out1=reads1.fq out2=reads2.fq ow
Executing jgi.ReformatReads [-Xmx5g, in=qsortA.bam, out1=reads1.fq, out2=reads2.fq, ow]

Set INTERLEAVED to true
Found samtools 1.9
Input is being processed as paired
Exception in thread "Thread-2" java.lang.AssertionError: K00315:202:H577WBBXY:7:1101:1438:38767 13      -1      -       57328513        57328662        1000000100000000011     1       0  60       TCCAGAAAGAGAGTTCTTGAAAGAAAGAGGCGCTATCATTTTGACACAGATGGNAAGGGCTCGATTCACGATCAAAAAGGCTCCAAAAANAAAAAAAATTCTGTTTGTTNTNCTTTTTCCTCCNGATTCTAGTTTTTTTNTTATNTTTTG  AAFFFFFAJJJJJ7FFFJJ7JJ7FJJJJJFFJJJJJFJJJJJF-A-AAFFJFA!AF---AAJJ-AJF<FF-AAFAF<AAAJJJJFA-7A!--77-<A<JJ<A--7A<FJ!J!<FJJJJA<FF7!)7AJF---7AFJJFJ!JF--!7-7<F      .       30      C68m28Nm53      .

K00315:202:H577WBBXY:7:1101:1438:38767  113     19      57328582        60      68S82M  1       176081441       0       CAAAANATAANAAAAAAACTAGAATCNGGAGGAAAAAGNANAACAAACAGAATTTTTTTTNTTTTTGGAGCCTTTTTGATCGTGAATCGAGCCCTTNCCATCTGTGTCAAAATGATAGCGCCTCTTTCTTTCAAGAACTCTCTTTCTGGA      F<7-7!--FJ!JFJJFA7---FJA7)!7FF<AJJJJF<!J!JF<A7--A<JJ<A<-77--!A7-AFJJJJAAA<FAFAA-FF<FJA-JJAA---FA!AFJFFAA-A-FJJJJJFJJJJJFFJJJJJF7JJ7JJFFF7JJJJJAFFFFFAA      NM:i:1  MD:Z:28C53      MC:Z:150M       AS:i:80 XS:i:21 SA:Z:1,176081357,+,81S62M7S,54,7;       RG:Z:FGC2033.7
flag=1110001

        at stream.SamReadInputStream.toReadList(SamReadInputStream.java:136)
        at stream.SamReadInputStream.fillBuffer(SamReadInputStream.java:89)
        at stream.SamReadInputStream.hasMore(SamReadInputStream.java:53)
        at stream.ConcurrentGenericReadInputStream$ReadThread.readLists(ConcurrentGenericReadInputStream.java:643)
        at stream.ConcurrentGenericReadInputStream$ReadThread.run(ConcurrentGenericReadInputStream.java:635)

I thought I could make read1.fq reads2.fq and use repair.sh Clearly it finds the two problem reads. I have tried other arguments with the reformat.sh, and a few other tools, but none want to correct the MATE_NOT_FOUND problem. It would seem impossible that I am the first person to run into this issue. Thank you very much for any advice/help.

-bwubb

[GenoMax]

I think you are looking for the following options that are accessible from "reformat.sh".

Sam and bam processing options:

mappedonly=f            Toss unmapped reads.
unmappedonly=f          Toss mapped reads.
pairedonly=f            Toss reads that are not mapped as proper pairs.
unpairedonly=f          Toss reads that are mapped as proper pairs.
primaryonly=f           Toss secondary alignments.  Set this to true for sam to fastq conversion.
minmapq=-1              If non-negative, toss reads with mapq under this.
maxmapq=-1              If non-negative, toss reads with mapq over this.
requiredbits=0          (rbits) Toss sam lines with any of these flag bits unset.  Similar to samtools -f.
filterbits=0            (fbits) Toss sam lines with any of these flag bits set.  Similar to samtools -F.

[raquelnb]

Hello Brian, it is just a question that pop up when using bbmap. For the outm, what type of flags do you use to both take pairs where both reads mappedand those where just one read of the pair mapped? Or how do you get read of both 77 and 141 flags?

Thank you and best