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 01:40 PM
visualisation tool for paired-end, to find the mate? Fad2012 Bioinformatics 7 06-03-2013 11:06 AM
Tool for obtaining counts for paired-end strand specific RNAseq data? JenBarb RNA Sequencing 2 03-26-2013 10:59 AM
Getting insertion counts from DNA sequence data anjulka Bioinformatics 0 02-27-2013 01:14 PM
Which tool to find prokaryotic ribosome in unknown genome sequences? marmue60 RNA Sequencing 0 06-08-2012 02:34 AM

Reply
 
Thread Tools
Old 09-12-2014, 10:02 AM   #1
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default tool to find how many counts for each sequence

Hi all,

It is my first time to do gene sequencing. I am sequencing short PCR products and want to get the counts for each sequence (I know the sequences in the sample for PCR). Is there anyone knew any tool for this?

Thank you for the help.

Zheng
lyw1 is offline   Reply With Quote
Old 09-12-2014, 10:10 AM   #2
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

How long are the sequences, and do you expect 100% identity, or what? Do you already know the sequences you are expecting?
Brian Bushnell is offline   Reply With Quote
Old 09-12-2014, 10:33 AM   #3
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default 30 to 90 bp.

Yes, I expect 100% identity and know the sequences.
lyw1 is offline   Reply With Quote
Old 09-12-2014, 10:38 AM   #4
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

In that case, just map the reads to a reference fasta file you make of the sequences. For example, with BBMap, you can output a file containing the counts of reads that mapped to each reference sequence:

bbmap.sh -Xmx2g ref=sequences.fasta in=reads.fastq scafstats=stats.txt

What kind of sequencing are you doing? You may need to adapter-trim the input reads first, for example, if you have a fixed read length with NGS reads.
Brian Bushnell is offline   Reply With Quote
Old 09-12-2014, 10:44 AM   #5
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Thanks a lot. I use MiSeq.
lyw1 is offline   Reply With Quote
Old 09-12-2014, 10:52 AM   #6
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

In that case, you will need to either trim the adapters first, or use a different approach such as matching substrings, if the reference sequences are sufficiently different.

You can trim adapters like this:

bbduk.sh -Xmx1g in=reads.fq out=clean.fq ktrim=r ref=nextera.fa.gz,truseq.fa.gz k=25 mink=12 hdist=1

...where those nextera and truseq adapter files are included with the BBTools package (which contains BBMap and BBDuk) in the /resources/ directory.

Alternately, you can generate a report directly from BBDuk like this:

bbduk.sh -Xmx1g in=reads.fastq ref=sequences.fasta stats=stats.txt k=29 fbm

...if the reference sequences are sufficiently different that they do not share any 29-mers.
Brian Bushnell is offline   Reply With Quote
Old 09-12-2014, 10:55 AM   #7
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

If I add a few samples that include inserts with 30 bp unknown sequences, is it possible for me to count reads for 50 sequences that appear most frequently?
lyw1 is offline   Reply With Quote
Old 09-12-2014, 11:02 AM   #8
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

The above tools require known sequences. If some are unknown, you could assemble the adapter-trimmed reads (with, say, Velvet) which should result in a fasta file containing one copy of each sequence present in your library, then map to that assembly to get the counts, as in post 4.
Brian Bushnell is offline   Reply With Quote
Old 09-12-2014, 11:13 AM   #9
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Brian, it is very helpful. Thanks again.
lyw1 is offline   Reply With Quote
Old 09-12-2014, 11:44 AM   #10
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Oh, I should probably add, this may be difficult to assemble since the coverage is likely to be extremely high. You may need to normalize first. Also, what length are the reads, and are they paired or single?
Brian Bushnell is offline   Reply With Quote
Old 09-12-2014, 12:35 PM   #11
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

The reads are from both sides and 30 for each side.
lyw1 is offline   Reply With Quote
Old 09-12-2014, 03:09 PM   #12
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

In that case, there's an easier solution, if all of the reads of the same product are expected to start at the same location.

kmercountexact.sh in=reads.fq outk=counts.txt mincount=4 k=30

That will produce a file "counts.txt" containing the count of every 30bp string in the file, excluding strings present fewer than 4 times.
Brian Bushnell is offline   Reply With Quote
Old 09-16-2014, 02:58 PM   #13
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default It works good

until

zheng@zheng-XPS-8500:~/Desktop/bbmap/20140916ngs$ kmercountexact.sh in=ERa-05RNS160m_S6_L001_R2_001.fastq outk=ERa-05RNS160m_S6_L001_R2_001_counts16.txt mincount=3 k=16
java -ea -Xmx11628m -cp /home/zheng/Desktop/bbmap/current/ jgi.CountKmersExact in=ERa-05RNS160m_S6_L001_R2_001.fastq outk=ERa-05RNS160m_S6_L001_R2_001_counts16.txt mincount=3 k=16
Executing jgi.CountKmersExact [in=ERa-05RNS160m_S6_L001_R2_001.fastq, outk=ERa-05RNS160m_S6_L001_R2_001_counts16.txt, mincount=3, k=16]

ways=37
Initial:
Memory: free=187m, used=64m

Final:
Memory: free=125m, used=193m

Input: 396838 reads 11111464 bases.
Result: 0 reads (0.00%) 0 bases (0.00%)
Unique Kmers: 4345841
Exception in thread "main" java.lang.AssertionError: ACCTTTTACCTAGTTT 3

at kmer.KmerNode.dumpKmersAsText(KmerNode.java:155)
at kmer.HashForest.dumpKmersAsText(HashForest.java:211)
at kmer.HashArray.dumpKmersAsText(HashArray.java:249)
at jgi.CountKmersExact.dumpKmersAsText(CountKmersExact.java:749)
at jgi.CountKmersExact.process(CountKmersExact.java:383)
at jgi.CountKmersExact.main(CountKmersExact.java:58)

Here I tried k=25, 20, 16. It works for k=25 or 20,

Some sequences at both ends may show "N". Do this command search short sequence from different start positions? Or always start from a fixed position? thanks,

Zheng
lyw1 is offline   Reply With Quote
Old 09-16-2014, 03:28 PM   #14
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Oh, sorry about that; I left a debug assertion in. I will remove that from the next version. Until then, add the "-da" flag to your command and it will work correctly.

Any subsequence containing an 'N' will be ignored. Otherwise, it looks at all valid subsequences. So, for example:

NNTTANGNNCT

for K=2 would produce:
TT 1
TA 1
CT 1

but nothing with the G that is surrounded by Ns. Thus, the reason I suggest K=30 for length-30 reads is that you will produce exactly one substring per read, while K=20 would produce 11 substrings per read. But if a 30bp read contains even a single N, it will yield no substrings.
Brian Bushnell is offline   Reply With Quote
Old 09-16-2014, 04:20 PM   #15
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Thanks.

Seems still have a problem after with flag -da.

zheng@zheng-XPS-8500:~/Desktop/bbmap/20140916ngs$ kmercountexact.sh in=ERa-05RNS160m_S6_L001_R2_001.fastq outk=ERa-05RNS160m_S6_L001_R2_001_counts16.txt mincount=3 k=16 -da
java -da -Xmx11773m -cp /home/zheng/Desktop/bbmap/current/ jgi.CountKmersExact in=ERa-05RNS160m_S6_L001_R2_001.fastq outk=ERa-05RNS160m_S6_L001_R2_001_counts16.txt mincount=3 k=16 -da
Executing jgi.CountKmersExact [in=ERa-05RNS160m_S6_L001_R2_001.fastq, outk=ERa-05RNS160m_S6_L001_R2_001_counts16.txt, mincount=3, k=16, -da]

ways=37
Exception in thread "main" java.lang.RuntimeException: Output file ERa-05RNS160m_S6_L001_R2_001_counts16.txt already exists, and overwrite=false
at jgi.CountKmersExact.<init>(CountKmersExact.java:335)
at jgi.CountKmersExact.main(CountKmersExact.java:55)
lyw1 is offline   Reply With Quote
Old 09-16-2014, 07:11 PM   #16
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Oh - that's an intentional protection from overwriting files. Just delete the output file first or add the "overwrite" flag.
Brian Bushnell is offline   Reply With Quote
Old 09-17-2014, 09:44 AM   #17
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default high contaninants

Thanks.

Input is being processed as unpaired

Input: 385043 reads 10781204 bases.
Contaminants: 341911 reads (88.80%) 9573508 bases (88.80%)
Result: 43132 reads (11.20%) 1207696 bases (11.20%)

What is diffinition of contaminants? It looks very high.
lyw1 is offline   Reply With Quote
Old 09-17-2014, 10:02 AM   #18
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

I need to read 30 nt for sequences. Miseq read 32 nt in sequencing. Thus many sequences have NN at last 2 positions. Does this relate to high contaminant rate?
lyw1 is offline   Reply With Quote
Old 09-17-2014, 10:23 AM   #19
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Are you using bbduk.sh? That's the only one that prints anything about contaminants. Can you show your specific command line?

Anyway, if you tried filtering out adapters and you got a result like that, it means you have almost no product and mostly adapter sequence.
Brian Bushnell is offline   Reply With Quote
Old 09-17-2014, 10:28 AM   #20
lyw1
Member
 
Location: Portland, OR

Join Date: Sep 2014
Posts: 33
Default

Yes, bbduk.sh.

Input is being processed as unpaired

Input: 385043 reads 10781204 bases.
Contaminants: 341911 reads (88.80%) 9573508 bases (88.80%)
Result: 43132 reads (11.20%) 1207696 bases (11.20%)
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 09:22 PM.


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