GRanges strand identification

Question

[v.berriosfarias]

hello I have a strong headache with a task. The issue is that I wan to calculate RPM based on metatranscriptome read matches against its own Metagenome to search specific gene expression levels.

For that purpose performed an alignment with bowtie2 (aligning the metatranscriptomic reads against metagenomic contigs) and I obtained a bam file.

With that bam file I want to use bam2rpm function of "intePareto" R package to calculate the RPM for each gene of interest.

Previously I have mapped the genes of interest against the metagneomic contigs using mmseqs2 and I obtained a blastn table (output format 6: https://www.metagenomics.wiki/tools/blast/blastn-output-format-6)

Reading the documentation of bam2rpm() function, I need to build a GRanges object with the GRanges() function.

The problem is that last function: one of the needed parameters is the strand to look for when constructing the GRanges object per each contain and it should be added as input on the following manner:

gr1 <- GRanges(seqnames=Rle(c("ch1", "chMT"), c(2, 4)),
ranges=IRanges(16:21, 20), strand=rep(c("+", "-", "*"), 2))  

on the above code, the strand parameter ins asking for the strand to look for. To adores that, how do I know at which strand the program should look for per each contig? how can I retrieve that information? from the bam file? if so, how?

Thanks for your time :)

Answer 1

[benformatics]

If the you metatranscriptome data you have is stranded then you need to know what strandedness is relative to the genome you are using. This would be based on the RNA-seq protocol and is usually anti-sense or unstranded. If the protocol is unstranded (ideal for you) then I suggest you use '*' if it's not then you need to figure out the strand of your reads based on their alignments to the metatranscriptome. I would suggest you look at them in IGV to figure this out if you are not experienced in python/R.