This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Calculating the expression level of genes

Hi,

I have seq data in some time points. How I can plot the expression level of genes (individually or together) so that the gene is expression more, less or equal to mean of expression level??? I though about something like heat map that colours before and after mean are lighter and darker respectively. But I don't know how to deal with this

gene expression rna-seq r
• 4,728 views • link
written by Za Za

Some example data would help.

• 1 views • link
written by cpad0112 cpad0112 
head(a[,1:10])

Actually, I want to map the level of genes expression in each time point around the mean or central tendencies

• 1 views • link
written by Za Za

some thing like this? Rplot02

• 1 views • link
written by cpad0112 cpad0112

Oh thank you. I think yes.

• 1 views • link
written by Za Za

In above image:

  1. Each box represents time point (time 1, 2 and 3 consecutively)
  2. Each box has expression value on Y-axis and Samples on X-axis.
  3. Each dot represents an expression value. Expression value above mean (for that sample) are represented by red and below mean are cyan.
  4. The black line in each box touches all the means of samples within the time point. It is not over all mean, for all the samples in the time point. I am not sure of %CV here. I used mean as cut off for color.

In OP data, most of them are zeros or near zero. I create a data frame. If you could post following information that would be great:

  1. Is each cell and sample are same? for eg there are about 200 cells. Is each cell is the sample it self or in each cell, there are multiple samples?
  2. Are these cells same across time points?
  3. If possible, post non-zero values for few of the samples and few of the time points to get an idea.
• 1 views • link
written by cpad0112 cpad0112 
data$mpg_type <- ifelse(data$mpg_z < mean "below", "above")  # above / below avg flag
data <- data[order(data$mpg_z), ]  # sort
data$`Gene` <- factor(data$`Gene`, levels = data$`Gene`)  # convert to fa
• 1 views • link
written by Za Za

I simulated data frame with gene names DDB_G0267178 for 10 genes, time points for h10-16 and 12 samples (s1-12). Following are graphs using ggplot. Since you mentioned violin polots, I tried recreating violin plots. Unfortunately, simulated data set is small. Hence When I faceted, the violins are out of shape. Hence I included image without faceting. I included both mean and mean+se. Points are colored based on the mean. The colors can be customized. Change the code as you wish. Images are at the end:

Data preparation and per sample faceting code

## Simulate data frame with 10 gene names and create a matrix of 10rx 12 c expression data
df1=data.frame(genes = paste0("DDB_G0267",sample(c(178:200), 10, replace = F)), matrix(sample(round(rnorm( 120, 4, 4), 2)), 10, 12))
## Change the names of the columns
colnames(df1)[-1] = paste0("sample", seq(1, 12))

## create a data frame for 12 sample names and 4 time points
pheno= data.frame(samples = paste0("sample", seq(1, 12)),time = rep(c("h10","h12","h14","h16"), each = 3))

## Append expression means of each sample to the pheno data
samples=data.frame(samples=colnames(df1[-1]), sample_means=colMeans(df1[-1]), stringsAsFactors = F, time=pheno$time)
geno=data.frame(genes=df1[1],gene_means=rowMeans(df1[-1]), stringsAsFactors = F)

## Load library tidyr
library(tidyr)

## convert wide format to long format
df2=gather(df1,"samples","expression",-"genes")

## Merge the metadata with data
df3=merge(df2,samples, by="samples")
df4=merge(df3,geno, by="genes")
df5=df4 %>% group_by(genes,time) %>% mutate(gene_time_mean=mean(expression))

## Plot using ggplot2. There are 3 graphs below. First one is with faceting for each  time point
## Faceting per sample with mean

s1=ggplot(df3, aes(samples, expression))+
    geom_violin()+
    geom_jitter(aes(color=expression>sample_means), width = 0.2)+
    scale_colour_manual(values=c("darkgreen", "red"))+
    stat_summary(fun.y=mean, geom="point", aes(group=1), size=1)+
    stat_summary(color="black",  geom = "line", aes(group=1))+
    facet_grid(~time, scales = "free")+
    theme_bw()+
    theme(
        legend.position = "none",
        strip.text.x = element_text(size = 24),
        strip.background = element_blank(),
        axis.text.x = element_text(size = 14, angle=90),
        axis.text.y = element_text(size = 24),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()
        )

## without faceting per sample with mean
s2=ggplot(df3, aes(samples, expression))+
    geom_violin()+
    geom_jitter(aes(color=expression>sample_means), width = 0.2)+
    scale_colour_manual(values=c("darkgreen", "red"))+
    stat_summary(fun.y=mean, geom="point", aes(group=1), size=1)+
    stat_summary(color="black",  geom = "line", aes(group=1))+
    theme_bw()+
    theme(
        legend.position = "none",
        strip.text.x = element_text(size = 24),
        strip.background = element_blank(),
        axis.text.x = element_text(size = 14, angle=90),
        axis.text.y = element_text(size = 24),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()
    )

## without faceting per sample with mean and se
s3=ggplot(df3, aes(samples, expression))+
    geom_violin()+
    geom_jitter(aes(color=expression>sample_means), width = 0.2)+
    scale_colour_manual(values=c("darkgreen", "red"))+
    stat_summary(fun.y=mean, geom="point", aes(group=1), size=1)+
    stat_summary(color="black",  geom = "smooth", aes(group=1))+
    theme_bw()+
    theme(
        legend.position = "none",
        strip.text.x = element_text(size = 24),
        strip.background = element_blank(),
        axis.text.x = element_text(size = 14, angle=90),
        axis.text.y = element_text(size = 24),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()
    )
• 1 views • link
written by cpad0112 cpad0112

@OP: you can also use heatmap like this:

heatmap

• 0 views • link
written by cpad0112 cpad0112

Excuse me, what is the black dot in graphs?

• 1 views • link
written by Za Za

black dot is mean per each sample.

• 1 views • link
written by cpad0112 cpad0112

Thank you, so the values pf green dots even placed under block dot, are higher than mean?

• 1 views • link
written by Za Za

No. any dot above mean is above mean and any dot below mean is below. In this case, red dots are above mean and green dots are below mean. Seems I was plotting sample means for gene means. I have updated the code and images as well. Thank you for the question.

• 1 views • link
written by cpad0112 cpad0112

Sorry, in new code 

aes(color=expression<gene_means)

in old one 

aes(color=expression<means)

when I am running new one by replacing gene_means with means gives me 

Error in expression < mean : 
  comparison (3) is possible only for atomic and list types

When I am running new code without any change, gives me 

Error in FUN(X[[i]], ...) : object 'gene_means' not found
• 1 views • link
written by Za Za

Per gene faceting requires a data frame than above. In earlier code, I was highlighting per sample mean on per gene data. This confusion was/is due to my code. I updated the code. Now the code shows per sample faceting. I will add per gene faceting.

• 1 views • link
written by cpad0112 cpad0112

Sorry, ,ay you please look at this picture

I used your code for one time point and only two genes. So if I am not wrong, only in few cells (x axis) coloured by green these two genes are expressing less than mean in this time point

https://ibb.co/dcb53T

• 1 views • link
written by Za Za

.....................................................

• 1 views • link
written by Za Za

You already have time sample data with time points. You don't have to create that data frame. Since I do not have that data, i used it. Use the data you already have to create the data.

• 1 views • link
written by cpad0112 cpad0112

Thank you, you are right

• 1 views • link
written by Za Za

Thank you, I plotted as you are plotting

I just wondering, how some people could be this much kind dedicating their precious time for helping people...

Thank you again

• 1 views • link
written by Za Za

Thanks a lot for your immense kindness,

May I have something like this for heatmap?

The scale of colour in heat map likely shows the difference with mean well

• 1 views • link
written by Za Za

For heat map:

ggplot(df3, aes(samples,genes, fill = expression)) +
    geom_tile(colour = "black") +
    facet_grid(~time,scales ="free") +
    scale_fill_gradient(low="red", high="green") +
    labs(x="sample",
         y="")+
    theme_bw()+
    theme(
        strip.text.x = element_text(size = 24),
        strip.background = element_blank(),
        axis.text.x = element_text(size = 18, angle=90),
        axis.text.y = element_text(size = 18),
        axis.title.x = element_text(size = 24),
        axis.title.y = element_text(size = 24),
    )

Rplot

• 1 views • link
written by cpad0112 cpad0112

Sorry, your last comment about heat map was vanished. I plotted for two genes and one time point as picture. So, only few cells (x axis) are expressing these two genes less than mean in this time point

• 1 views • link
written by Za Za

I replaced previous heatmap pic with new one with code. There was no other text other than image. Image you have attached seems to be too busy and mean for many of them is incorrect. For mean to be in the middle of the graph, it should almost equal number of points above and below it. Please look into your code and data again.

• 1 views • link
written by cpad0112 cpad0112

Sorry, you are plotting 10 genes for each sample. I just plotted one gene per sample (samples are my cells, so that I have 209 samples (cells) for time point h16). If I plot more than a gene, I am not able to see which of them in a sample is expressiong above or below the mean.

• 1 views • link
written by Za Za

I agree with you that simulated data is small and uniform. In general, real life big data (such as *omics, BFSI) is complicated and one needs several tricks to make sense of data. I can only help you with simulated data as my english is poor and cannot understand the description of data. I rather would like to see the data. I have updated the code with per gene, per sample tiling for time points. Please go through the code. Do not emulate it, rather see how a final dataframe is created for plotting. Then start plotting.

For other two points that are raised (gene expression and sample number), first issue is rather a statistical issue and second one is a non-issue.

• 1 views • link
written by cpad0112 cpad0112

Excuse me,

I have cells for one time point and the expression of one gene for these cells ,

• 1 views • link
written by Za Za

...continuation of data preparation and sample faceting code post...

## with faceting, gene expression data with mean
g1=ggplot(df5, aes(genes, expression))+
    geom_violin()+
    geom_jitter(width = 0.2, aes(color=expression>gene_time_mean))+
    scale_colour_manual(values=c("darkgreen", "red"))+
    stat_summary(fun.y=mean, geom="point", aes(group=1), size=1)+
    stat_summary(color="steelblue",  geom = "line", aes(group=1), size=1.2)+
    theme_bw()+
        theme(
        legend.position = "none",
        strip.text.x = element_text(size = 24),
        strip.background = element_blank(),
        axis.text.x = element_text(size = 14, angle=90),
        axis.text.y = element_text(size = 24),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()
    )+
    facet_grid(~time)

## without faceting, gene expression data with mean only
g2=ggplot(df5, aes(genes, expression))+
    geom_violin()+
    geom_jitter(width = 0.2, aes(color=expression>gene_means))+
    scale_colour_manual(values=c("darkgreen", "red"))+
    stat_summary(fun.y=mean, geom="point", aes(group=1), size=1)+
    stat_summary(color="steelblue",  geom = "line", aes(group=1), size=1.2)+
    theme_bw()+
    theme(
        legend.position = "none",
        strip.text.x = element_text(size = 24),
        strip.background = element_blank(),
        axis.text.x = element_text(size = 14, angle=90),
        axis.text.y = element_text(size = 24),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()
    )

## without faceting, gene expression data with mean + se

g3=ggplot(df5, aes(genes, expression))+
    geom_violin()+
    geom_jitter(width = 0.2, aes(color=expression>gene_means))+
    scale_colour_manual(values=c("darkgreen", "red"))+
    stat_summary(fun.y=mean, geom="point", aes(group=1), size=1)+
    stat_summary(color="steelblue",  geom = "smooth", aes(group=1), size=1.2)+
    theme_bw()+
    theme(
        legend.position = "none",
        strip.text.x = element_text(size = 24),
        strip.background = element_blank(),
        axis.text.x = element_text(size = 14, angle=90),
        axis.text.y = element_text(size = 24),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()
    )

library(gridExtra)
grid.arrange(s1,s2,s3,g1,g2,g3, nrow = 2, ncol=3)

Rplot01

• 1 views • link
written by cpad0112 cpad0112

Thanks a lot

This is a link of this data

209, 153, 179 and 205 samples for each time point respectively.

• 1 views • link
written by Za Za

0 answers

No answers yet.

Log in to answer this question.

More posts like this