Extracting features above a certain coverage

[annavilborg]

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]

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!

[jparsons]

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.

[annavilborg]

This worked very well! Thank you!

[hanshart]

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:

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

[annavilborg]

Neat! Thank you!