Read alignment plot for RNA-Seq data

Question

[mike]

I am having RNA-Seq data with 5 replicates per condition. I want to look at the coverage (read alignment) for a particular genes for these condition. Whats the best possible way to take average of 5 replicates per condition and check the coverage? I have already checked it for a single replicate with Gviz package in R.

Thanks

[GouthamAtla]

You can make a box plot or just plot the points across all conditions for each replicate. DESeq2 has an option called plotCounts which does the same.

I use a wrapper around plotCounts to plot the boxplot ( as I have more than 20 samples per condition ) for multiple genes,

plotGeneCounts <- function (dds, genes, intgroup) {
  for (i in 1:length(genes)) { 
    x <- plotCounts(dds, gene=genes[i], intgroup=intgroup, transform=TRUE, returnData=TRUE)
    x$gene=genes[i]
    if (exists("gene_vec")){
      gene_vec <- rbind(gene_vec,x)
    } else {
      gene_vec <- x
    }
  }
  ggplot(gene_vec, aes(x=factor(intgroup), y=count))+geom_boxplot(aes(fill=gene))+facet_wrap(~gene, scales="free")
}

Use:

genes <- c("INS","IAPP","PDX1")
plotGeneCounts(dds,genes,"treat")

Where dds is a DESeq2 object.

EDIT: I understand that you would like to show the coverage across the gene structure. For this, you could simply create the track hub for ucsc browser for all your samples and show the coverage as overlay mode, separately for treated and untreated.

[mike]

Thanks but I am looking for this image. I was able to generate this with Gviz but only with 1 replicate per condition. How can I take read counts for all replicates per condition and generate the same plot?

Answer 1

[mike]