Mapping RNAseq counts to genomic loci

Question

[Ssri]

Hello.

I am trying to figure out a way to map count reads from RNAseq datafiles to genomic loci and display this data in an elegant way (ideally circos plots). However, most of the tools I have come across (edgeR, DEseq2) require a control sample and report differential gene expression. What I am looking for is a gene expression profile for bacteria that have sequenced transcriptomes to understand how gene expression varies across the genome. There is no "test" sample. Is there a tool that can assist in providing this information, given the sequenced genome and RNAseq counts? Thank you!

[seidel]

Yes, a genome browser can display what you're looking for, but you'd have to convert your data to a bedGraph file to draw the coordinates of the genes, and the heights according to their expression values. You could experiment with IGV from MIT for this. You mention circos plots...you could also use these if you have gene coordinates and expression values. Here's a toy example for 10 genes in a 360 base coordinate space:

library(circlize)
# make one sector called "a"
df <- data.frame(sectors=rep("a",2), x=c(1,360))

circos.par("track.height"=0.2)
circos.initialize(df$sectors, x=df$x)
# create a track
circos.track(df$sectors, ylim=c(0, 10))

# make up some gene coordinates
rpoints <- sample(1:360,20)
rpoints <- rpoints[order(rpoints)]
# make up some expression values on a scale of 1 to 10
rvalues <- runif(10,1,10)

# set our gene starts and ends by choosing
# the even or odd indexes of our sorted coordinates
starts <- rpoints[seq(1,20, by=2)]
ends <- rpoints[seq(2,20, by=2)]

y_bottom <- rep(0,length(x_left))
y_top <- rvalues

# draw genes as boxes, with size = expression value
circos.rect(starts, y_bottom, ends, y_top, col="lightblue")

# gene names
gnames <- paste0("g", seq_along(starts))
midpoints <- starts + (ends - starts)/2
circos.text(midpoints, 0.2+mm_y(-3), gnames, sector.index="a", track.index=1, cex=0.8)

[Ssri]

Thank you so much for your feedback! This looks great! I can try this out. Thank you :)

[seidel]

Could you make a genome browser coverage track? In this case, you would set the coverage of your features to be the normalized RNA Seq counts. Alternatively you could make this as a circos plot using R or python, as well as a few others I imagine. What tools/languages are you looking to use? What are you familiar with?

[Ssri]

I'm not sure I understand what a genome browser coverage track is. Is this kind of like UCSC Genome Browser? I would like to stick to tools or packages that are already being used, but if unavailable I am comfortable with R and would like to stay within R.