Chromosome information for miRNA

Question

Apologies if this question is previously asked. I have a BAM file alignment to hg19 and miRNA gff3 file from mirbase for annotation of BAM file. I want to extract information in the following format (e.g.):

miRNA        chr    chr_start    chr_end    expression_count    sequence
hsa-mir-155  chr2       ...         ....         1020            ATGC....

I know that I can use htseq-count for BAM annotation and for expression counts but how can I get their chromosomal location. The reason for doing this specific way is that I want to show them in a CIRCOS plot.

Just counting is not the goal, their chromosomal location also needs to be extracted. So just subread/htseq-count is not sufficient for extracting genomic location information.

Previously, I had the same idea about processing the sam records. If I want to annotate chromosomal information with miRNA using .gff file, I have to check for each line in .sam records whether their chr_start and chr_end lie in the bandwidth of start and end of any miRNA in .gff file. This process can be very heavy and slow because .sam records are of more than 5 GB. That's why I want to know if there is any extremely fast and generalized way to do it.

I want to avoid processing sam records in python because my experience with extremely heavy data like fastq or sam in python is not good. I would rather prefer awk/bash or any low-level language.

Thanks.

Answer 1

you can do samtools view pathtobamfile, this will convert the file to .sam format and send it to stdout, the columns can then be parsed out which contain the information you are after.

Answer 2

The gff should have the coordinates for each feature. I'd use featureCounts from subread to count up how many reads hit each feature.