[Brian Bushnell]

Oh, looks like the problem is the commas. There should not be any commas in the command line. BBTools is interpreting "SL118119_S5_L001_R2_001.fastq.gz," as a filename, and it does not know what format "fastq.gz," is so it defaults to sam. I'll change it to default to fastq for that program.

[Canadian_philosophy]

Hi Brian,

Thanks for responding.

1) That is what I assumed, wanted to hear it from someone else here

2) These are things I am trying with bbsplit:

Trial 1 :

 
java -Djava.library.path=/data/davelab/dr153/softwares/bbmap/jni/ -ea -Xmx56G -cp /data/davelab/dr153/softwares/bbmap/current/ align2.BBSplitter ow=t fastareadlen=500 minhits=1 minratio=0.56 maxindel=20 qtrim=rl untrim=t trimq=6
in1=/data/davelab/projects/Xenomousie/fastq_link/SL118119_S5_L001_R1_001.fastq.gz, in2=/data/davelab/projects/Xenomousie/fastq_link/SL118119_S5_L001_R2_001.fastq.gz, out_hg19=SL118119_t4_hg19_gatkorder_R1.fastq, out_hg19=SL118119_t4_hg19_gatkorder_R2.fastq, out_mm10=SL118119_t4_mm10_gatkorder_R1.fastq, out_mm10=SL118119_t4_mm10_gatkorder_R2.fastq, ref=/data/davelab/bix/resources/genomes/hg19/ucsc.hg19.fasta,/data/davelab/bix/resources/genomes/mm10/mm10_gatkorder.fa, ambiguous2=toss

Trial 2:

java -Djava.library.path=/data/davelab/dr153/softwares/bbmap/jni/ -ea -Xmx56G -cp /data/davelab/dr153/softwares/bbmap/current/ align2.BBSplitter ow=t fastareadlen=500 minhits=1 minratio=0.56 maxindel=20 qtrim=rl untrim=t trimq=6 in1=/data/davelab/projects/Xenomousie/fastq_link/SL118119_S5_L001_R1_001.fastq.gz, in2=/data/davelab/projects/Xenomousie/fastq_link/SL118119_S5_L001_R2_001.fastq.gz, out_hg19_1=SL118119_t3_hg19_gatkorder_R1.fastq, out_hg19_2=SL118119_t3_hg19_gatkorder_R2.fastq, out_mm10_1=SL118119_t3_mm10_gatkorder_R1.fastq, out_mm10_2=SL118119_t3_mm10_gatkorder_R2.fastq, ref=/data/davelab/bix/resources/genomes/hg19/ucsc.hg19.fasta,/data/davelab/bix/resources/genomes/mm10/mm10_gatkorder.fa, ambiguous2=toss

3) Thanks for clarifying the bam issue. I saw the fastq as recommended output

Thanks,
Deepthi

[Brian Bushnell]

1) The difference is only in highly conserved regions. If you first map to mouse, then map the unmapped reads to human, human reads from regions with high identity to mouse will already have been removed, because they mapped to the mouse. When mapping to both genomes at once, the reads will go to the organism they actually came from because they (should) match that one better.

2) Can you give your entire command line? I just tried it and verified that I got fastq output when I used this command:

bbsplit.sh in=reads.fq basename=o_%.fq

3) Sorry, that's kind of inevitable. Actually I do not recommend using sam output from BBSplit, only fastq output. Sam output will have the ambiguously-mapped reads listed as mapped to only one organism if you are in ambig=all mode. For example, if a read mappes to human and mouse equally well, it will be output in both sam files, but in each case it will be listed as mapped to the same location (which is either human or mouse).

[Canadian_philosophy]

BBMap

1) I am interested in understanding how BBmap is simultaneously aligning reads to multiple genomes? Is this different from using bwa to align to multiple genomes sequencially?

2) I used Bbsplit to split my fastq to hg19 reads and mouse reads. The output was only in samformat ( although I tried to name the o/p as _1.fastq.gz and 2.fastq.gz). Is it possible to directly get fastq outputs after aligning to multiple genomes?
with something like this :
out_hg19_1=SL118119_t3_hg19_gatkorder_R1.fastq, out_hg19_2=SL118119_t3_hg19_gatkorder_R2.fastq, out_mm10_1=SL118119_t3_mm10_gatkorder_R1.fastq, out_mm10_2=SL118119_t3_mm10_gatkorder_R2.fastq

3) The sam output header for both the .mm10.sam and .hg19.sam had both mouse and human chromosomes. I had to manually edit this for samtools to do post processing. Is there a way out of this?

[JeanLove]

Can you please post the exact command you used, and the entire screen output?

getreads.sh in=af_svi_0.05Perc_22633.fasta out=bla.fasta id=6272,19594

java -ea -Xmx200m -cp /common/WORK/bbmap/current/ jgi.GetReads in=af_svi_0.05Perc_22633.fasta out=bla.fasta id=6272,19594
Executing jgi.GetReads [in=af_svi_0.05Perc_22633.fasta, out=bla.fasta, id=6272,19594]

Input is unpaired
java.lang.AssertionError: Attempting to read from a closed file. Current header: MISEQ:70:000000000-A4FTJ:1:2111:9262:20597 1:N:0:
at stream.FastaReadInputStream.nextBases(FastaReadInputStream.java:357)
at stream.FastaReadInputStream.generateRead(FastaReadInputStream.java:237)
at stream.FastaReadInputStream.fillList(FastaReadInputStream.java:210)
at stream.FastaReadInputStream.nextList(FastaReadInputStream.java:121)
at stream.ConcurrentGenericReadInputStream$ReadThread.readLists(ConcurrentGenericReadInputStream.java:656)
at stream.ConcurrentGenericReadInputStream$ReadThread.run(ConcurrentGenericReadInputStream.java:635)
Time: 0.167 seconds.
Reads Processed: 19595 117.09k reads/sec
Bases Processed: 4621k 27.62m bases/sec
Exception in thread "main" java.lang.RuntimeException: GetReads terminated in an error state; the output may be corrupt.
at jgi.GetReads.<init>(GetReads.java:321)
at jgi.GetReads.main(GetReads.java:37)

I get the output file and it seems normal, but I don't know why do I get this error.

[Brian Bushnell]

Can you please post the exact command you used, and the entire screen output?

[JeanLove]

If you know the sequence numbers, you can use getreads.sh:

I have tried it, but i get the java error:

Exception in thread "main" java.lang.RuntimeException: GetReads terminated in an error state; the output may be corrupt.
at jgi.GetReads.<init>(GetReads.java:321)
at jgi.GetReads.main(GetReads.java:37)

what could be the cause? the input file is normal .fasta file with 22000 reads from Illumina MiSeq.

[Brian Bushnell]

If you know the sequence numbers, you can use getreads.sh:

Usage: getreads.sh in=<file> id=<number,number,number...> out=<file>

The first read (or pair) has ID 0, the second read (or pair) has ID 1, etc.

Parameters:
in=<file> Specify the input file, or stdin.
out=<file> Specify the output file, or stdout.
id= Comma delimited list of numbers or ranges, in any order.
For example: id=5,93,17-31,8,0,12-13

[JeanLove]

It is not clear if you know their identifiers (or just their positions). If you have the identifiers then you can do the following[/url]

I do not know identifiers..

In case BBMap does not work you can try faSomeRecords utility from Jim Kent @UCSC which is described here: http://seqanswers.com/forums/showthread.php?t=64004

thanks for your help, i'll look faSomeRecords utility up!

[GenoMax]

hi! is there a way to extract specific sequences from fasta file, regarding their positions in the file using BBMap tools? I know the numbers of their name lines and want to extract whole sequences to a file.

It is not clear if you know their identifiers (or just their positions). If you have the identifiers then you can do the following

By default, "filterbyname" discards reads with names in your name list, and keeps the rest. To include them and discard the others, do this:

$ filterbyname.sh in=003.fastq out=filter003.fq names=names003.txt include=t

Even though this example is for fastq files I am reasonsably certain that it would work for fasta. Make sure your file names end in .fa to facilitate that.

In case BBMap does not work you can try faSomeRecords utility from Jim Kent @UCSC which is described here: http://seqanswers.com/forums/showthread.php?t=64004

[JeanLove]

hi! is there a way to extract specific sequences from fasta file, regarding their positions in the file using BBMap tools? I know the numbers of their name lines and want to extract whole sequences to a file.

[GenoMax]

Hi,

can you suggest how to get raw reads mapped per gene from bbmap?

Thanks,

S.

Based on your other posts you appear to have discovered htseq-count and featureCounts. BBMap suite does not have a read counting program (though @Brian has one on the wishlist of things to add to BBMap).

[susanklein]

reads mapped per gene

Hi,

can you suggest how to get raw reads mapped per gene from bbmap?

Thanks,

S.

[GenoMax]

Unless your scheduler is buffering the entire output in some temp location you should have seen trimmed files right away (have you ever looked before if a job does that)? Did you try top -H to see all ~40 java threads?

Based on your last comment that must be the case.

[jazz710]

Based on Brian's comments, I bet I have run out of memory. My reference is very large, and a 3*K increase in memory would put me over my server capacity. I will try to re-run with a smaller reference.