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
Junior Member
 
Location: Paris, France

Join Date: Apr 2017
Posts: 9
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,553
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
Junior Member
 
Location: Paris, France

Join Date: Apr 2017
Posts: 9
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,302
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
Junior Member
 
Location: Paris, France

Join Date: Apr 2017
Posts: 9
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,302
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: 23
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,553
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: 23
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
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 05:33 PM.


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