[Brian Bushnell]

For K=3 and the input file:

>
AAAAA

You would want an output file:

>
AAA
>
AAA
>
AAA

Is that correct? I don't have anything that will do that; sorry. What did you want to use it for?

[jweger1988]

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]

@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.

[GenoMax]

@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.

[HESmith]

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.

[jweger1988]

@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?

[HESmith]

@jweger1988, what tools have you tried for VCF annotation and what problems have you encountered?

[jweger1988]

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.

[darthsequencer]

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

[Brian Bushnell]

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

[darthsequencer]

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?

[sk8bro]

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.

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...

@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

[GenoMax]

If you wanted to use dedupe with PE files you could do it this way

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:

dedupe=t optical=f
All duplicates are detected, whether optical or not.  All copies except one are removed for each duplicate.

[jweger1988]

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

[GenoMax]

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.