SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
BBMap (aligner for DNA/RNAseq) is now open-source and available for download. Brian Bushnell Bioinformatics 556 09-07-2017 03:42 PM
BBMap for BitSeq dietmar13 Bioinformatics 1 04-30-2015 08:40 AM
Please help my BBMap cannot remove Illumina adapter TofuKaj Bioinformatics 4 04-28-2015 08:53 AM
BBMap Error Phage Hunter Bioinformatics 5 01-14-2015 04:34 AM
Introducing BBMap, a new short-read aligner for DNA and RNA Brian Bushnell Bioinformatics 24 07-07-2014 09:37 AM

Reply
 
Thread Tools
Old 04-12-2016, 03:20 AM   #101
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,439
Default

Quote:
Originally Posted by susanklein View Post
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).
GenoMax is offline   Reply With Quote
Old 05-28-2016, 02:08 AM   #102
JeanLove
Junior Member
 
Location: croatia

Join Date: Jan 2016
Posts: 5
Default

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.
JeanLove is offline   Reply With Quote
Old 05-28-2016, 04:45 AM   #103
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,439
Default

Quote:
Originally Posted by JeanLove View Post
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

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

Code:
$ 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
GenoMax is offline   Reply With Quote
Old 05-28-2016, 07:26 AM   #104
JeanLove
Junior Member
 
Location: croatia

Join Date: Jan 2016
Posts: 5
Default

Quote:
Originally Posted by GenoMax View Post
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..


Quote:
Originally Posted by GenoMax View Post
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!
JeanLove is offline   Reply With Quote
Old 05-28-2016, 07:38 AM   #105
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,668
Default

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

Quote:
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
Brian Bushnell is offline   Reply With Quote
Old 05-29-2016, 12:59 AM   #106
JeanLove
Junior Member
 
Location: croatia

Join Date: Jan 2016
Posts: 5
Default

Quote:
Originally Posted by Brian Bushnell View Post
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.
JeanLove is offline   Reply With Quote
Old 05-29-2016, 04:01 AM   #107
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,668
Default

Can you please post the exact command you used, and the entire screen output?
Brian Bushnell is offline   Reply With Quote
Old 05-29-2016, 05:06 AM   #108
JeanLove
Junior Member
 
Location: croatia

Join Date: Jan 2016
Posts: 5
Default

Quote:
Originally Posted by Brian Bushnell View Post
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.

Last edited by JeanLove; 05-29-2016 at 05:08 AM. Reason: forgot to write something
JeanLove is offline   Reply With Quote
Old 06-03-2016, 12:26 PM   #109
Canadian_philosophy
Junior Member
 
Location: Canada

Join Date: Apr 2010
Posts: 9
Default 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?
Canadian_philosophy is offline   Reply With Quote
Old 06-03-2016, 01:05 PM   #110
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,668
Default

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).
Brian Bushnell is offline   Reply With Quote
Old 06-03-2016, 01:20 PM   #111
Canadian_philosophy
Junior Member
 
Location: Canada

Join Date: Apr 2010
Posts: 9
Default

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 :
Code:
 
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:
Code:
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
Canadian_philosophy is offline   Reply With Quote
Old 06-03-2016, 01:24 PM   #112
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,668
Default

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.
Brian Bushnell is offline   Reply With Quote
Old 06-03-2016, 01:33 PM   #113
Canadian_philosophy
Junior Member
 
Location: Canada

Join Date: Apr 2010
Posts: 9
Default

Hmmm.. so I pasted that from the stderr/out file of the program, maybe that is how you formatted it there? This is how my command looks
Code:
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 interleaved=True basename=o5%_#.fq ref=/data/davelab/bix/resources/genomes/hg19/ucsc.hg19.fasta,/data/davelab/bix/resources/genomes/mm10/mm10_gatkorder.fa refstats=refstats.txt ambiguous2=toss
The only comma is in between ref files. Does this look okay to you?

Also, thanks for fixing the default input.
Canadian_philosophy is offline   Reply With Quote
Old 06-03-2016, 01:40 PM   #114
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,668
Default

Hmmm, that command looks fine. Just to clarify, is that the specific command that produced sam output?
Brian Bushnell is offline   Reply With Quote
Old 06-06-2016, 06:23 AM   #115
Canadian_philosophy
Junior Member
 
Location: Canada

Join Date: Apr 2010
Posts: 9
Default

Hi Brian,
This command below worked for me. It looks like if I name the output files, that doesn't work (e.g out1_hg19=,out2_hg19=,out1_mm10= etc)
.
If I use the basename command, the fastqs look fine!

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 interleaved=True basename=o5%_#.fq ref=/data/davelab/bix/resources/genomes/hg19/ucsc.hg19.fasta,/data/davelab/bix/resources/genomes/mm10/mm10_gatkorder.fa refstats=refstats.txt ambiguous2=toss
Canadian_philosophy is offline   Reply With Quote
Old 07-22-2016, 05:36 AM   #116
moistplus
Member
 
Location: Germany

Join Date: Feb 2016
Posts: 40
Default

Is it possible to compute the insert size (and make a histogram distribution) from a .bam file ?
moistplus is offline   Reply With Quote
Old 07-22-2016, 06:02 AM   #117
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,439
Default

Quote:
Originally Posted by moistplus View Post
Is it possible to compute the insert size (and make a histogram distribution) from a .bam file ?
I think so (use unsorted bam). If reads are spliced you will get some odd results. @Brian may have a suggestion later on how to handle that.

Make sure samtools is available in your $PATH.

Code:
reformat.sh in=your_file.bam ihist=hist.txt
GenoMax is offline   Reply With Quote
Old 07-22-2016, 08:59 AM   #118
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,668
Default

As long as the file is unsorted or name-sorted, that command will work. Spliced reads shouldn't cause a problem, unless the splice junction is in the unsequenced area between the reads, of course.
Brian Bushnell is offline   Reply With Quote
Old 09-07-2016, 10:15 AM   #119
JVGen
Member
 
Location: East Coast

Join Date: Jul 2016
Posts: 37
Default

Quote:
Originally Posted by Brian Bushnell View Post
As long as the file is unsorted or name-sorted, that command will work. Spliced reads shouldn't cause a problem, unless the splice junction is in the unsequenced area between the reads, of course.
Hi Brian and Geno, I didn't want to start a new thread, so I figured I'd quote and post to get your attention.

I started using BBDeDupe and Normalization in Geneious on Merged Paired Reads. These reads are not yet assembled, but will be with De Novo after DeDupe/Norm. I'm doing it because I have about 50,000 reads that cover a 10kb PCR amplicon, and I think the depth is causing assembly issues.

I wasn't able to find much information on DeDupe. Does it remove all 100% matching reads that are encompassed in another read, thereby reducing my coverage to 1? If so, is there a way to increase my coverage, to say 50-100 instead?

I'm getting a few areas with low coverage, less than 10, and I generally require that I have a minimum coverage of 15 to call a base in the consensus sequence, otherwise I call a gap. I generally have >500-1000 coverage in all places, but occasionally I get some places that have <10 coverage. I don't think this is amplicon bias, rather I think it is contaminating sequence. I use a core service for the illumina sequencing, but nearly everyone that submits samples is working on HIV, hence my need for a coverage threshold.

Thoughts on this?

Thanks,
Jake
JVGen is offline   Reply With Quote
Old 09-07-2016, 10:56 AM   #120
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,668
Default

Hi Jake,

By default, Dedupe will remove all but one copy of reads that are fully contained within other reads. This may not be a good idea prior to assembly of merged reads (since many will be contained by others due to the variable insert size); just using normalization alone with a target of 100x may work better. You can set Dedupe to only remove reads that are identical to others, though, by adding the flag "ac=f".

If you suspect your data is contaminated, I suggest you try to decontaminate it prior to normalization/assembly. If you are working on non-HIV and it is contaminated with HIV, that's fairly straightforward; you can use BBSplit or Seal (if you have both references) or just BBMap (if you only have the HIV reference but nothing for your amplicon) to separate the reads by mapping. Alternately, if you are confident that all of the low-coverage stuff is contamination, you can assemble the low-coverage portion of your data using Tadpole, which allows coverage thresholds (e.g. "mincr=1 maxcr=15" will make it assemble only those contigs with coverage between 1x and 15x). Then use the resulting assembly to filter all of the reads via mapping, keeping only the unmapped reads.

Incidentally, Tadpole tends to do a good job assembling very-high-coverage short sequences.
Brian Bushnell is offline   Reply With Quote
Reply

Tags
bbmap

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 02:15 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO