SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Which tool should I use to summarize counts from BWA-mem? Birdman Bioinformatics 0 03-24-2014 02:40 PM
visualisation tool for paired-end, to find the mate? Fad2012 Bioinformatics 7 06-03-2013 12:06 PM
Tool for obtaining counts for paired-end strand specific RNAseq data? JenBarb RNA Sequencing 2 03-26-2013 11:59 AM
Getting insertion counts from DNA sequence data anjulka Bioinformatics 0 02-27-2013 02:14 PM
Which tool to find prokaryotic ribosome in unknown genome sequences? marmue60 RNA Sequencing 0 06-08-2012 03:34 AM

Reply
 
Thread Tools
Old 12-23-2014, 03:23 PM   #41
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Hi Brian,

I found some sequences that don't have a known leading sequence. Those sequences are amplified in PCR but they are not my target sequence. Is it possible for me to pick all the sequences that include same leading sequences (at 5' end or 3' end)? Then I can remove the leading sequence to get my target sequences. Thanks,

Zheng
lyw1 is offline   Reply With Quote
Old 12-23-2014, 03:50 PM   #42
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Quote:
Originally Posted by lyw1 View Post
Hi Brian,

I found some sequences that don't have a known leading sequence. Those sequences are amplified in PCR but they are not my target sequence. Is it possible for me to pick all the sequences that include same leading sequences (at 5' end or 3' end)? Then I can remove the leading sequence to get my target sequences. Thanks,

Zheng
Zheng,

Yes, that is possible, if you know the sequence. You can use BBDuk for filtering a specific sequence at the beginning of the read like this:

bbduk.sh -Xmx1g in=reads.fq outm=matched.fq outu=unmatched.fq restrictleft=25 k=25 literal=AAAAACCCCCTTTTTGGGGGAAAAA

In this case, all reads starting with "AAAAACCCCCTTTTTGGGGGAAAAA" will end up in "matched.fq" and all other reads will end up in "unmatched.fq". Specifically, the command means "look for 25-mers in the leftmost 25 bp of the read", which will require an exact prefix match, though you can relax that if you want.

So you could bin all the reads with your known sequence, then look at the remaining reads to see what they have in common. You can do the same thing with the tail of the read using "restrictright" instead, though you can't use both restrictions at the same time.
Brian Bushnell is offline   Reply With Quote
Old 12-23-2014, 04:42 PM   #43
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Thanks, Brian. That is great.

Zheng
lyw1 is offline   Reply With Quote
Old 12-23-2014, 05:05 PM   #44
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Hi Brian,

There is bug here,

Exception in thread "main" java.lang.RuntimeException: Unknown parameter restrictleft=10
at jgi.BBDukF.<init>(BBDukF.java:427)
at jgi.BBDukF.main(BBDukF.java:66)
lyw1 is offline   Reply With Quote
Old 12-25-2014, 01:09 PM   #45
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Zheng,

Are you using the latest version (34.19) of BBTools? I added "restrictleft" and "restrictright" only a few weeks ago, so an old version won't support those flags...
Brian Bushnell is offline   Reply With Quote
Old 12-29-2014, 02:46 PM   #46
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Hi Brian,

I download the newest edition (34.44) and test. It showed the same problem.

Exception in thread "main" java.lang.RuntimeException: Unknown parameter restrictleft=25
at jgi.BBDukF.<init>(BBDukF.java:427)
at jgi.BBDukF.main(BBDukF.java:66)

Thanks,

Zheng
lyw1 is offline   Reply With Quote
Old 12-29-2014, 05:30 PM   #47
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Zheng,

You appear to have downloaded an old version by mistake. The latest is now 34.22, and there is no 34.44 (only a 33.44); so please download 34.22 - I have verified that it is working correctly.
Brian Bushnell is offline   Reply With Quote
Old 03-24-2015, 03:18 PM   #48
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Hi Brian,

Some sequence has quality problelm. It may have one or two "N" or one nucleotide showed wrongly. Is it possible to count this kind sequence as the correct sequence? For example, aaggctctggattacaggat is the reference sequence. Now, is it possible to count aaggctcNggattacaggat to the same group of the reference sequence?

Thank you,

Zheng
lyw1 is offline   Reply With Quote
Old 03-24-2015, 03:23 PM   #49
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Hi Zheng,

Yes, you can do this with the "hdist=1" flag, which will allow up to 1 substitution error.

-Brian
Brian Bushnell is offline   Reply With Quote
Old 03-26-2015, 06:21 PM   #50
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Thanks, Brian.

Currently the max of hdist is 3. Is there a way to allow hdist=6?
lyw1 is offline   Reply With Quote
Old 03-26-2015, 06:45 PM   #51
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Zheng,

A high value of hdist requires exponentially more time and memory to load the reference, but you can us an arbitrarily high value if you want, by adding the -da flag, which disables assertions.

For example:

bbduk.sh in=reads.fq outm=matched.fq ref=ref.fa k=25 hdist=6 -da

However, the program will generate a factor of hdist^(3*K) additional kmers - for k=25 and hdist=6, taking 177,978,515,625 times as much time and memory to build the index, so that won't work. hdist=4 would probably still work though if you gave it enough memory.

One alternative is "meet in the middle". I added a new flag "qhdist" which stands for "query hamming distance", which mutates the query kmers instead of the reference kmers. So you would get equivalent results from the above command and this command:

bbduk.sh in=reads.fq outm=matched.fq ref=ref.fa k=25 hdist=3 qhdist=3

...but the second command would require 421,875 times less memory, so it's feasible to run. It might be very slow, though, as it would also run 421,875 times slower than normal.

There is one last alternative, though - I added support in BBDuk for degenerate bases at specific locations. You can use it like this:

bbduk.sh in=reads.fq outm=matched.fq literal=NNNNNNCCCCGGGGGTTTTTAAAAA k=25 copyundefined

With the "copyundefined" flag, a copy of each reference sequence will be made representing every valid combination of defined letter. So instead of increasing memory or time use by 6^75, it only increases them by 4^6 or 4096 which is completely reasonable, but it only allows substitutions at predefined locations. You can use the "copyundefined", "hdist", and "qhdist" flags together for a lot of flexibility - for example, hdist=2 qhdist=1 and 3 Ns in the reference would allow a hamming distance of 6 with much lower resource requirements than hdist=6. Just be sure to give BBDuk as much memory as possible.

Edit: It might be fun to write something that handles arbitrary motifs with degenerate bases and a set number of mismatches, with no time or memory penalty. Maybe I'll do that; I'll let you know.

Last edited by Brian Bushnell; 03-26-2015 at 06:51 PM.
Brian Bushnell is offline   Reply With Quote
Old 04-15-2015, 05:49 PM   #52
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Thanks, Brian. Very helpful.
lyw1 is offline   Reply With Quote
Old 04-27-2015, 07:02 PM   #53
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default comfuse about result

Hi Brian,

The total sequence length is 26 bp.

I used following command to get sequence with starting sequence "CACTTCTATAGT".

bbduk.sh -Xmx1g in=7_S7_L001_R1_001.fastq outm=7_S7_L001_R1_001matched.fq outu=7_S7_L001_R1_001unmatched.fastq restrictleft=12 k=12 literal=CACTTCTATAGT

Then I counted seqeunce frequency with following command,

kmercountexact.sh in=7_S7_L001_R1_001matched.fastq outk=7_S7_L001_R1_001_24new.txt mincount=26 k=26

From the first command, I saw every sequence include CACTTCTATAGT at leftest side when I opened 7_S7_L001_R1_001matched.fq

However, 7_S7_L001_R1_001_24new.txt result surprised me. Most of them did not include CACTTCTATAGT at its leftest side.

>32
CCCCCCCCCCCGCCACTATAGAAGTG
>30
TCCGGGGCTTCCCCACTATAGAAGTG
>72
CACTTCTATAGTGGGGAACCACCGTT
>101
CACTTCTATAGTGGGGATTTCCCCTT
>36
TGGGGAATTCCCTAACTATAGAAGTG
>40
TTGGGGAAGCCCCCACTATAGAAGTG
>33
TATGGGGAATCCCCACTATAGAAGTG
>47
TGGGGCTTTCCCATACTATAGAAGTG
>46
TAGGGGATTTCCCTACTATAGAAGTG
>27
CACTTCTATAGTGGGGAAGTACCGAT
>28
CACTTCTATAGTCGGTGGTTCCCCTT
>34
TCCCTTCTATTTCTACTATAGAAGTG
>231
TAGGTGCTTCCCCTACTATAGAAGTG
>52
TTGGGGGAAACCCCACTATAGAAGTG
>59
TCGGGATTTCCCCAACTATAGAAGTG
>43
TCGGGAAAACCCCCACTATAGAAGTG
>45
TTGGGATTTCCCCCACTATAGAAGTG
>58
CACTTCTATAGTGGGGAATCCCCCTT
>26
TCCCCCACCGTCCCACTATAGAAGTG
>61
CACTTCTATAGTAGGGAAACCACCTT
>26
TCCCCACTGGAGCTACTATAGAAGTG
>65
CGGGAATCCCCTAAACTATAGAAGTG
>46
TCGGGGATGCACCTACTATAGAAGTG

Please help. Thanks.

Zheng
lyw1 is offline   Reply With Quote
Old 04-27-2015, 08:05 PM   #54
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Hi Zheng,

There are 3 compounding issues here that make things a little confusing. First, BBDuk, by default looks for both a kmer and its reverse-complement. You can disable this behavior with "rcomp=f" (though that probably won't affect anything in this case). Second, BBDuk also ignores the middle base of a kmer in order to increase sensitivity. You can disable this with "mm=f" (which stands for maskmiddle). Lastly, and most important in this case, kmercountexact stores only a kmer or its reverse complement, whichever is alphabetically higher. I forgot to add that to the help, but this too can be disabled with "rcomp=f". This shouldn't change the counts in your case, but it will make the orientation of the output the same as the input.

It looks like all of the sequences you posted either start with "CACTTCTATAGT" or end with "ACTATAGAAGTG", which means the 26-mer was stored as its reverse-complement.

Last edited by Brian Bushnell; 04-27-2015 at 08:08 PM.
Brian Bushnell is offline   Reply With Quote
Old 04-28-2015, 02:17 PM   #55
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Thanks, Brian.

for command,
bbduk.sh -Xmx1g in=7_S7_L001_R1_001matched.fq out=7_S7_L001_R1_001matched_trimmed.fq ftr=2
The out file has nothing. ftr=10 and it works. Some sequences has a couple N at right side because of poor quality. How can I remove 2 bp at right side of the sequence?

Zheng
lyw1 is offline   Reply With Quote
Old 04-28-2015, 02:22 PM   #56
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Hi Zheng,

"ftr=2" will remove all bases after position 2 in the read (zero-based), meaning all reads will end up 3bp long, which will make them get discarded due to the default minlength=10. To remove the last 2 bases, you need "ftr2=2" instead. And you can set "minlength=0" if you don't want any reads to be discarded.

-Brian
Brian Bushnell is offline   Reply With Quote
Old 04-28-2015, 02:43 PM   #57
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Hi Brian,

Executing jgi.BBDukF [-Xmx1g, in=7_S7_L001_R1_001matched.fq, out=7_S7_L001_R1_001matched_trimmed.fq, ftr2=2]

Exception in thread "main" java.lang.RuntimeException: Unknown parameter ftr2=2
at jgi.BBDukF.<init>(BBDukF.java:392)
at jgi.BBDukF.main(BBDukF.java:62)

Is anything wrong here?

It seems it works to remove last 2 bp from 26 sequence by ftr=23?

Thanks,
lyw1 is offline   Reply With Quote
Old 04-28-2015, 02:46 PM   #58
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Zheng,

ftr=23 would work. The reason ftr2=2 failed is because you have an old version of BBTools; I added it fairly recently.

-Brian
Brian Bushnell is offline   Reply With Quote
Old 05-01-2015, 07:01 PM   #59
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Thanks, Brian. I test new version of BB Tools and it works.
lyw1 is offline   Reply With Quote
Reply

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:16 PM.


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