SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Bowtie, an ultrafast, memory-efficient, open source short read aligner Ben Langmead Bioinformatics 513 05-14-2015 02:29 PM
Introducing BBMap, a new short-read aligner for DNA and RNA Brian Bushnell Bioinformatics 24 07-07-2014 09:37 AM
Miso's open source joyce kang Bioinformatics 1 01-25-2012 06:25 AM
Targeted resequencing - open source stanford_genome_tech Genomic Resequencing 3 09-27-2011 03:27 PM
EKOPath 4 going open source dnusol Bioinformatics 0 06-15-2011 01:10 AM

Reply
 
Thread Tools
Old 05-11-2017, 12:42 PM   #501
jweger1988
Member
 
Location: Paris, France

Join Date: Apr 2017
Posts: 34
Default

Brian,

Thanks for the response. My issue was exactly that. It seems to be working fine now by only using one at a time. The other issue was a hamming distance thing that I also figured out.

Will you be incorporating the soft-clipped base trimming in this next update?

Thanks again.
jweger1988 is offline   Reply With Quote
Old 05-11-2017, 03:55 PM   #502
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,669
Default

Quote:
Originally Posted by Mchicken View Post
Thanks a lot for your helpful replies. It is now working perfectly . One more question to the "xstag". I see that if the gene is on the + strand, R2 maps on+ and R1 maps on -. If I use "xstag=ss" I get for spliced R1 and R2 reads XS=+. So I think this is the correct library type. Am I right?
I think so. I've always found Tophat's description of those terms to be very confusing, so I have trouble remember for more than a few days which is the correct one to use and when.

Quote:
Originally Posted by jweger1988
Will you be incorporating the soft-clipped base trimming in this next update?
Yes, I have added that as of v37.23 which I just released. The flag is "trimclip", for example:

bbduk.sh in=mapped.sam out=trimmed.sam trimclip ordered

You can also output them as fastq if you want.
Brian Bushnell is offline   Reply With Quote
Old 05-16-2017, 04:01 AM   #503
jweger1988
Member
 
Location: Paris, France

Join Date: Apr 2017
Posts: 34
Default

Hi Brian,

I have a couple of questions.

I've been comparing different alignment programs and then some different variant calling software. For some reason when I use a bam file produced with bbmap and try to call variants with lofreq I get an error message, but not with bam files produced with other programs (bowtie, mosaic, etc...). Even after sorting the bam I get the same message. Any thoughts on why this might be?

Also, I'm wondering if it's possible to use dedupe with the input as either paired fastq files (i.e. in1 and in2) or a bam file?

Thanks,
James
jweger1988 is offline   Reply With Quote
Old 05-16-2017, 04:29 AM   #504
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,440
Default

Not sure if this is applicable but BBMap produces SAM v.1.4 CIGAR strings by default. If your variant calling program expects v.1.3 format then you can either add sam=1.3 option when you align (or reformat.sh in=your.bam out=new.bam sam=1.3).

Dedupe can be used with PE files but only if reads are interleaved. You may want to use clumpify instead of dedupe.
GenoMax is offline   Reply With Quote
Old 05-16-2017, 08:11 AM   #505
jweger1988
Member
 
Location: Paris, France

Join Date: Apr 2017
Posts: 34
Default

Genomax,

Awesome, changing CIGAR fixed the issue. Thanks.

About Clumpify vs Dedupe, that's a good idea. I remember Brian saying somewhere that clumpify is not quite as good as Dedupe? What would be your recommended settings?

Thanks
jweger1988 is offline   Reply With Quote
Old 05-16-2017, 08:29 AM   #506
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,440
Default

If you wanted to use dedupe with PE files you could do it this way
Code:
reformat.sh in1=R1.fq in2=R2.fq out=stdout.fq | dedupe.sh in=stdin.fq out=stdout.fq | reformat.sh in=stdin.fq out1=new1.fq out2=new2.fq
Clumpify is more flexible compared to dedupe. You can select what happens to the duplicates. You may be interested in this option:
Code:
dedupe=t optical=f
All duplicates are detected, whether optical or not.  All copies except one are removed for each duplicate.
GenoMax is offline   Reply With Quote
Old 05-17-2017, 01:07 PM   #507
sk8bro
Member
 
Location: boston

Join Date: Feb 2012
Posts: 26
Default minid and idfilter

Hi Brian,

Thank you for making the deterministic flag. Results in that regard look great.

I am confused by how minid and idfilter are working, however. I will email you a small dataset that recreates the 'unexpected' behavior I am encountering if that is helpful.

My expectation is that this command should only allow alignments where BOTH the forward and reverse read, independently align with at least 75% pairwise identity to the reference.... ie overlap in the reads is NOT taken into account and a contig pairwise identity to the reference is NOT what the setting is referring to. Note: I usually output 4 fastqs for mapped/unmapped pairs but have changed to sam output and specified the idtag flag to give more info.
Code:
bbsplit.sh\
 -Xmx23g\
 averagepairdist=200\
 deterministic=t\
 k=8\
 minid=0.75\
 idfilter=0.75\
 maxindel=20\
 nzo=f\
 po=t\
 ambiguous2=split\
 ref=test_ref\
 idtag=t\
 in=Input_1P.fastq\
 in2=Input_2P.fastq\
 outm=Mapped.sam\
 outu=Unmapped.sam\
Here is the head of Mapped.sam... how are these mapping? YI:f:73.10, YI:f:64.14 etc...
Code:
@HD     VN:1.4  SO:unsorted
@SQ     SN:006_Koxy_tonB_PC     LN:317
@SQ     SN:022_NDM_PC   LN:267
@PG     ID:BBMap        PN:BBMap        VN:37.22        CL:java -Djava.library.path=/install/bbmap/jni/ -ea -Xmx23g align2.BBSplitter ow=t fastareadlen=500 minhits=1 minratio=0.56 maxindel=20 qtrim=rl untrim=t trimq=6 -Xmx23g averagepairdist=200 deterministic=t k=8 minid=0.75 idfilter=0.75 maxindel=20 nzo=f po=t ambiguous2=split ref=test_ref idtag=t in=Input_1P.fastq in2=Input_2P.fastq outm=Mapped.sam outu=Unmapped.sam
gi|828959694|gb|CP011636.1|_Klebsiella_oxytoca_strain_CAV1374,_complete_genome_(reversed)-1-146 83      006_Koxy_tonB_PC        209     18      4=2X2=37I100=   =       6       -311    ATCAAGCTGTTTGCCGGGAATAACAACCAGCGCGGGGCGGCGTTTGACGTTACCGGGGCGCTGGACGATAACGATCGCGTGGCGGCGCGCTTAAGCGGCATGACCCGCTATGCAGACTCGCAGTTTGATACCTTAAAAGAGCAGC  ?????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????  NM:i:39 AM:i:18 YI:f:73.10
gi|828959694|gb|CP011636.1|_Klebsiella_oxytoca_strain_CAV1374,_complete_genome_(reversed)-1-146 163     006_Koxy_tonB_PC        6       9       87=1X2=47I1X1=1X3=2X    =       209     311     AATGATGGGCGATACCAACTCGCACAGCTCGCTGGTGGTCGATCCGTGGTTCCTGGAAAATATCGAAGTGGTGCGCGGCCCGGCCTCAGTGCTGTACGGCCGCTCTTCGCCCGGCGGCATCGTCGCCCTCACCTCGCGTAAACCC  ?????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????  NM:i:52 AM:i:9  YI:f:64.14
gi|828959694|gb|CP011636.1|_Klebsiella_oxytoca_strain_CAV1374,_complete_genome_(reversed)-1-146 83      006_Koxy_tonB_PC        209     18      4=2X2=37I100=   =       6       -311    ATCAAGCTGTTTGCCGGGAATAACAACCAGCGCGGGGCGGCGTTTGACGTTACCGGGGCGCTGGACGATAACGATCGCGTGGCGGCGCGCTTAAGCGGCATGACCCGCTATGCAGACTCGCAGTTTGATACCTTAAAAGAGCAGC  ?????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????  NM:i:39 AM:i:18 YI:f:73.10
gi|828959694|gb|CP011636.1|_Klebsiella_oxytoca_strain_CAV1374,_complete_genome_(reversed)-1-146 163     006_Koxy_tonB_PC        6       9       87=1X2=47I1X1=1X3=2X    =       209     311     AATGATGGGCGATACCAACTCGCACAGCTCGCTGGTGGTCGATCCGTGGTTCCTGGAAAATATCGAAGTGGTGCGCGGCCCGGCCTCAGTGCTGTACGGCCGCTCTTCGCCCGGCGGCATCGTCGCCCTCACCTCGCGTAAACCC  ?????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????  NM:i:52 AM:i:9  YI:f:64.14
gi|828959694|gb|CP011636.1|_Klebsiella_oxytoca_strain_CAV1374,_complete_genome_(reversed)-1-146 83      006_Koxy_tonB_PC        209     18      4=2X2=37I100=   =       6       -311    ATCAAGCTGTTTGCCGGGAATAACAACCAGCGCGGGGCGGCGTTTGACGTTACCGGGGCGCTGGACGATAACGATCGCGTGGCGGCGCGCTTAAGCGGCATGACCCGCTATGCAGACTCGCAGTTTGATACCTTAAAAGAGCAGC  ?????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????  NM:i:39 AM:i:18 YI:f:73.10
gi|828959694|gb|CP011636.1|_Klebsiella_oxytoca_strain_CAV1374,_complete_genome_(reversed)-1-146 163     006_Koxy_tonB_PC        6       9       87=1X2=47I1X1=1X3=2X    =       209     311     AATGATGGGCGATACCAACTCGCACAGCTCGCTGGTGGTCGATCCGTGGTTCCTGGAAAATATCGAAGTGGTGCGCGGCCCGGCCTCAGTGCTGTACGGCCGCTCTTCGCCCGGCGGCATCGTCGCCCTCACCTCGCGTAAACCC  ?????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????  NM:i:52 AM:i:9  YI:f:64.14
Thanks, Kate

Last edited by GenoMax; 05-17-2017 at 03:39 PM.
sk8bro is offline   Reply With Quote
Old 05-23-2017, 10:48 AM   #508
darthsequencer
Member
 
Location: california

Join Date: Feb 2012
Posts: 29
Default

Quote:
Originally Posted by darthsequencer View Post
Great! Yeah I'd like to know the duplicate membership for containments too.
Hi Brian - I think "renamecluster=t" in dedupe.sh will do what I want with regards to tracking which contigs are greater than "minidentity".

I run dedupe.sh like this: dedupe.sh in=file.fa out=file.fa mindentity=99

Will adding "renameclusters=t" change what dedupe is doing besides renaming the sequences to what cluster they belong?
darthsequencer is offline   Reply With Quote
Old 05-23-2017, 11:28 AM   #509
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,669
Default

No, it doesn't change the behavior otherwise. But for clustering you may want to disable duplicate removal (or you won't get accurate cluster sizes) with "am=f ac=f" and you need to add some clustering flags, "fo c pc".

So the command might look like:

dedupe.sh in=file.fa out=file.fa mindentity=99 am=f ac=f fo c pc
Brian Bushnell is offline   Reply With Quote
Old 05-23-2017, 03:46 PM   #510
darthsequencer
Member
 
Location: california

Join Date: Feb 2012
Posts: 29
Default

Quote:
Originally Posted by Brian Bushnell View Post
No, it doesn't change the behavior otherwise. But for clustering you may want to disable duplicate removal (or you won't get accurate cluster sizes) with "am=f ac=f" and you need to add some clustering flags, "fo c pc".

So the command might look like:

dedupe.sh in=file.fa out=file.fa mindentity=99 am=f ac=f fo c pc
Hi I gave that a shot and dedupe throws a bunch of errors and produces about double the clusters than dedupe does without clustering.

Here's a copy of log: https://www.dropbox.com/s/t8aqj7fgf2...er.fa.log?dl=0
darthsequencer is offline   Reply With Quote
Old 05-31-2017, 08:53 AM   #511
jweger1988
Member
 
Location: Paris, France

Join Date: Apr 2017
Posts: 34
Default

Hi Brian,

Is there a tool in the suite that will take a .gff or other annotation file and annotate a .vcf from callvariants.sh? By annotation I mean adding the result of the variant (i.e. AA change). I've tried to do it with a bunch of other tools but I'm having trouble.

I'm sure it's not currently supported, but is there any possibility that one could callvariants from a .gbk or .gff reference as opposed to a fasta to get the coding changes directly?

Thanks in advance.
jweger1988 is offline   Reply With Quote
Old 05-31-2017, 09:04 AM   #512
HESmith
Senior Member
 
Location: Washington DC

Join Date: Oct 2009
Posts: 486
Default

@jweger1988, what tools have you tried for VCF annotation and what problems have you encountered?
HESmith is offline   Reply With Quote
Old 05-31-2017, 09:10 AM   #513
jweger1988
Member
 
Location: Paris, France

Join Date: Apr 2017
Posts: 34
Default

@HESmith, thanks for the response.

I've tried Annovar, SnpEff and a few others. Do you have any experience with these or other recommendations? It's totally possibly I'm just messing this up.

I am working with a virus that is not in any of the database, I tried to build my own file in the database and was getting an error message about the file not having the correct chromosomes set (with SnpEff).

Any clues?
jweger1988 is offline   Reply With Quote
Old 05-31-2017, 09:18 AM   #514
HESmith
Senior Member
 
Location: Washington DC

Join Date: Oct 2009
Posts: 486
Default

Not trying to be a jerk, but have you tried variant calling with the test data that are included with the software? Success would eliminate installation and command errors as the source of your problem.

Also, a general rule for online troubleshooting of software problems is 1) to provide the exact command syntax that your used, and 2) report the exact error message that you received.
HESmith is offline   Reply With Quote
Old 05-31-2017, 09:20 AM   #515
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,440
Default

@jweger1988: I second @HESmith's suggestion.

Create a new thread with any errors or problems you have encountered with annotation programs.

Much as we love BBMap suite there are always going to be things that you will need to use a different program for functionality not available in BBMap suite.
GenoMax is offline   Reply With Quote
Old 05-31-2017, 10:21 AM   #516
jweger1988
Member
 
Location: Paris, France

Join Date: Apr 2017
Posts: 34
Default

@HESmith, you weren't being a jerk at all. It's a fair point.

I've created a new thread. http://seqanswers.com/forums/showthr...975#post207975. Thanks for any help you can provide.
jweger1988 is offline   Reply With Quote
Old 07-17-2017, 12:56 PM   #517
jweger1988
Member
 
Location: Paris, France

Join Date: Apr 2017
Posts: 34
Default

Hi Brian,

I've been using kmercountexact and it's been very useful to give me the kmers with their counts.

I'm wondering if any of your tools has the capability to give me a list of all of the kmers present at given length regardless of being unique or not. Basically a list that would also include the redundant kmers without counts.

Thanks in advance for your help.
jweger1988 is offline   Reply With Quote
Old 07-17-2017, 01:17 PM   #518
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,669
Default

For K=3 and the input file:

Code:
>
AAAAA
You would want an output file:

Code:
>
AAA
>
AAA
>
AAA
Is that correct? I don't have anything that will do that; sorry. What did you want to use it for?
Brian Bushnell is offline   Reply With Quote
Old 07-17-2017, 11:20 PM   #519
jweger1988
Member
 
Location: Paris, France

Join Date: Apr 2017
Posts: 34
Default

Thanks for the reply. That is correct.

I have a virus that I introduced some degenerate nucleotides in to track bottlenecks.

I suppose I could just reformat to the area of the read I'm interested in and then just convert to fasta and use that.
jweger1988 is offline   Reply With Quote
Old 07-31-2017, 07:51 PM   #520
darthsequencer
Member
 
Location: california

Join Date: Feb 2012
Posts: 29
Default

Hi - I love that bbmap and its tools can directly make bam files. I noticed that it's using samtools with 8 threads. Is there a way to increase the number of threads?

Thanks
darthsequencer is offline   Reply With Quote
Reply

Tags
bbmap, metagenomics, rna-seq aligners, short read alignment

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 09:05 AM.


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