SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Samtools/htseq-count numbers don't add up. tirohia Bioinformatics 3 09-30-2014 12:15 AM
htseq-count low count problem gandalf886 Bioinformatics 3 08-23-2014 08:05 AM
difference between htseq-count and bedtools multicov camelbbs Bioinformatics 24 01-01-2014 08:44 PM
multiBamCov or htseq-count to count read per feature ? NicoBxl Bioinformatics 1 07-03-2012 03:05 AM
Count difference by htseq-cout and samtools DZhang Bioinformatics 2 07-03-2011 12:05 PM

Reply
 
Thread Tools
Old 09-28-2014, 09:34 PM   #1
tirohia
Member
 
Location: Auckland, NZ

Join Date: Nov 2011
Posts: 46
Default Isolating non coding RNA using bedtools/samtools/htseq-count

I am trying to isolate ncRNA from a bam file.

I have a accepted_hits.bam file that has been produced by tophat, which to the best of my knowledge contains all reads that have mapped to the reference genome, be it coding or non-coding.
I also have a bed file that I have made using getting the lengths of all the chromosomes, into a basic bed file, with only the required fields, i.e:
Quote:
Chr1 0 17971566
Chr10 0 15490421
Chr11 0 17082262
Chr12 0 13002763
Chr13 0 17372821
Chr14 0 12031045
Chr15 0 19001909
Chr16 0 11476499
Chr17 0 14675105
Chr18 0 14814751
Chr19 0 11145676
Chr2 0 15003700
Chr20 0 14375345
Chr21 0 14287830
Chr22 0 12982530
Chr23 0 20688064
Chr24 0 17702303
Chr25 0 14321294
Chr26 0 17153587
Chr27 0 10917718
Chr28 0 16406049
Chr29 0 14180785
Chr3 0 20519224
Chr4 0 13065296
Chr5 0 18825850
Chr6 0 18208465
Chr7 0 18704015
Chr8 0 21083321
Chr9 0 11598160
Unknow 0 177361715
I then used bedtools subtract with the above bed file and the gff file associated with the genome that I aligned to. My reasoning being that the gff file has all the features (genes) in it, so if I subtract those from the lengths of the various chromosomes, I'll be be left with everything outside of a coding region.
My plan was to then use
Code:
samtools view -L nonCodingRegions.bed accepted_hits.bam > nonCoding.bam
This, for some reasons gives me the same file. i.e. it outputs accepted_hits.bam into nonCoding.bam. Which isn't right. Am I using the -L option correctly or should I be using something else?

My other thought was to get the reads that are identified with the gff features using the -o option in htseq-count:
Code:
htseq-count -f bam -r name -o coding -t gene -i ID accepted_hits.bam
so that I cold then try and figure out some way of subtracting the sam file output by that from the accepted_hits.bam. This doesn't appear to work either as it only outputs a list of features that the reads have hit, with no way of identifying which reads hit that feature, like this:
Quote:
XF:Z:__no_feature
XF:Z:Achn175571
XF:Z:Achn175571
XF:Z:__no_feature
XF:Z:Achn008681
XF:Z:__no_feature
XF:Z:__no_feature
XF:Z:Achn339401
XF:Z:Achn339401
XF:Z:Achn050721
Is there something I'm missing in my attempt to get htseq-count to output the samfile of reads that map to features in the gff file?

Or alternatively, is there a better way to draw out the ncRNA?

Cheers
Ben.
tirohia is offline   Reply With Quote
Old 09-29-2014, 12:43 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Firstly, don't blindly accept anything that this method produces as being an ncRNA. Most of what you'll be getting are repeat regions that have some background expression (i.e., uninteresting noise).

Secondly, your samtools view -L ... command is correct, so it's likely a problem with your BED file. Have a look at that in IGV or another browser and make sure it seems OK.

As an aside, there are programs that have been written to find many classes of ncRNAs. You'd be strongly encouraged to look into them before trying to roll your own.
dpryan is offline   Reply With Quote
Old 09-29-2014, 10:49 PM   #3
tirohia
Member
 
Location: Auckland, NZ

Join Date: Nov 2011
Posts: 46
Default

Again, thanks for the answer.

Don't worry, I wasn't going to accept it all
Looking at my BED file, not sure how to figure out if something is/what is wrong.

These other programs, would you have any recommendations? I've come across a number of papers that haven't (that I've seen) linked to any software that they've used. There's a couple, like CPAT that I'm currently trying to get running, but if there are other options, I'd be very interested in having a look at them.

Cheers
Ben.
tirohia is offline   Reply With Quote
Old 09-30-2014, 12:13 AM   #4
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

piPipe for piRNAs, mirdeep2 for miRNAs. There are a few for lincRNAs as well, but I can't come up with the names off-hand.
dpryan 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:25 AM.


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