SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
cummeRbund 'features' error shengzizhang Bioinformatics 2 12-04-2013 05:43 AM
structural features of miRNA eng_sheringokhy General 0 04-17-2013 12:32 AM
NCBI Features bwubb Bioinformatics 2 11-09-2012 10:03 AM
BEDTools v2.9 - new tools/features quinlana Bioinformatics 4 06-22-2011 04:42 PM
Extracting base coverage and quality values from BAM files unagaswamy Bioinformatics 9 11-11-2010 06:55 AM

Reply
 
Thread Tools
Old 06-14-2013, 08:28 AM   #1
annavilborg
Junior Member
 
Location: Connecticut, US

Join Date: Feb 2013
Posts: 8
Default Extracting features above a certain coverage

Hi,

I have a bam file of all my accepted hits (output from tophat run on data from strand specific libraries/HiSeq) and an gtf file with my genes of interest for which I am trying to find potential antisense transcripts. I would like to create a list - preferably one that can be visualized in a genome browser - that shows all genes for which there are antisense reads in the accepted hits.bam file provided that there are more than a certain number of reads mapping to that gene. I can easily use bedtools intersect to get all genes that have antisense reads mapping to them, but this gives me way too many hits. I have also exported my resulting hits to excel to try to sort the file there (bioinformatics newbie that I am) but I can't get it back to a format that works for visualizing in a genome browser.

Thank you.
annavilborg is offline   Reply With Quote
Old 06-17-2013, 09:12 AM   #2
annavilborg
Junior Member
 
Location: Connecticut, US

Join Date: Feb 2013
Posts: 8
Default

Hi again,

I have come up with a possible solution, but if anyone had some input to give I would very much appreciate it. If I did bedtools coverage to get the antisense coverage for each of my genes, and then sorted the resulting file on the coverage value? And then extracted everything in the first x lines? Could that work?

I am unsure how to do the sorting and extracting. "sort" and "grep"? Any suggestions on the exact "phrasing"?

Thank you!
annavilborg is offline   Reply With Quote
Old 06-17-2013, 12:37 PM   #3
jparsons
Member
 
Location: SF Bay Area

Join Date: Feb 2012
Posts: 62
Default

Once you have a coverage file, you can sort based on the result with "sort". the -n flag will sort numbers numerically (eg: 10<20<100 instead of 10<100<11<20) and the -kX flag will allow you to choose which key/column you wish to do the actual sort on. -k 4 would choose the 4th. -r would allow you to reverse the output to get the largest numbers first instead of last.

I would suggest "head -X <filename>" to extract the first X lines, although i imagine grep could also do it relatively easily.
jparsons is offline   Reply With Quote
Old 06-21-2013, 09:27 AM   #4
annavilborg
Junior Member
 
Location: Connecticut, US

Join Date: Feb 2013
Posts: 8
Default

This worked very well! Thank you!
annavilborg is offline   Reply With Quote
Old 06-24-2013, 08:09 PM   #5
hanshart
Member
 
Location: Germany

Join Date: Nov 2011
Posts: 27
Default

Quote:
Originally Posted by annavilborg View Post
Hi,
... provided that there are more than a certain number of reads mapping to that gene...
Thank you.
Hi annavilborg,
you can additionally pipe your sorted output to awk to get only those lines (genes) with read count above your threshold. Figure out which column (field) stores the read counts and then e.g. for field 4 and Threshold 100 apply something like:
Code:
bedtools coverage ... | awk '{if ($4>100) print $0}' | sort -nr -k 4
$4 is the 4th field
100 the read count threshold
print $0 prints the content of the whole line
sort -nr -k 4 will then numerically, reverse sort all these gene lines with readcount above threshold based on the corresponding column (4).
hanshart is offline   Reply With Quote
Old 06-26-2013, 05:23 PM   #6
annavilborg
Junior Member
 
Location: Connecticut, US

Join Date: Feb 2013
Posts: 8
Default

Neat! Thank you!
annavilborg is offline   Reply With Quote
Reply

Tags
rna-seq coverage cutoff

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


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