This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Suggestion for Differential Expression analysis

Good morning,

My task is to perform a Differential Expression Analysis on 31 tissues of an individual to see variation among gene expression. I'm not interesting in all genes expressed but actually I have to focus my analysis only on 2 particular genes. I trimmed my raw reads and check the quality, I mapped the reads with HISAT2, I assembled my reads with StringTie and I generated two CSV files containing the count matrices for genes and transcripts with the prepDE.py script following this pipeline (https://ccb.jhu.edu/software/stringtie/index.shtml?t=manual). As I just said now I have to perform DEA for 2 genes on my 31 different tissues.

How can I proceed? Could you suggest me a pipeline to follow possibly using DESeq2?

Thank you in advance for your time and suggestion!

rna-seq
• 268 views • link
updated by Carlo Yague Carlo Yague • written by Lorenzo • Lorenzo

Could you suggest me a pipeline to follow possibly using DESeq2?

Yes, the DESeq2 manual at Bioconductor, this is all you need.

• 0 views • link
written by ATpoint ATpoint

Could you help me please? I'm a beginner in R and I'm encountering some problems.

> library(DESeq2)
> countData <- as.matrix(read.csv("gene_count_matrix.csv", row.names="gene_id"))
> colData <- read.csv("Pheno_data.csv", row.names=1)
> all(rownames(colData) %in% colnames(countData))
[1] TRUE
> countData <- countData[, rownames(colData)]
> all(rownames(colData) == colnames(countData))
[1] TRUE
> dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~ Tissue)
Note: levels of factors in the design contain characters other than
letters, numbers, '_' and '.'. It is recommended (but not required) to use
only letters, numbers, and delimiters '_' or '.', as these are safe characters
for column names in R. [This is a message, not a warning or an error]

Warning message: In DESeqDataSet(se, design = design, ignoreRank) : some variables in design formula are characters, converting to factors dds <- DESeq(dds) estimating size factors Note: levels of factors in the design contain characters other than letters, numbers, '_' and '.'. It is recommended (but not required) to use only letters, numbers, and delimiters '_' or '.', as these are safe characters for column names in R. [This is a message, not a warning or an error] estimating dispersions Error in checkForExperimentalReplicates(object, modelMatrix) : 

   The design matrix has the same number of samples and coefficients to fit,
   so estimation of dispersion is not possible. Treating samples
  as replicates was deprecated in v1.20 and no longer supported since v1.22.
> head(colData)
                                   Tissue                   Sex
 Adrenal_Cortex         Adrenal Cortex      Female
 Cerebellum_Vermis  Cerebellum vermis Female
 Cornea                      Cornea                   Female
 Spleen                       Spleen                   Female
 > head(countData)
                                        Adrenal_Cortex Cerebellum_Vermis Cornea Spleen
 DNALI1|DNALI1                           0                   0                          0          0
 MSTRG.852|RAB3B                     0                  0                          64        4
 MSTRG.10382                           3933             2809                     2078   3393
 PRAM1|PRAM1                             0                  0                          0          0

How can I solve my problem??

• 0 views • link
written by Lorenzo • Lorenzo

Can you show the full coldata ?

the same number of samples and coefficients to fit

looks like you do not have replicates, which is triggering the error on the design matrix. If you do not have replicates, then you can not use DESeq2 for your analysis.

• 0 views • link
written by Carlo Yague Carlo Yague

I posted the full colData because now I'm not using all tissues because I want to set my pipeline and see if it works on only 4 tissues. My first problem is here:

> dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~ Tissue)

Note: levels of factors in the design contain characters other than letters, numbers, '_' and '.'. It is recommended (but not required) to use only letters, numbers, and delimiters '_' or '.', as these are safe characters for column names in R. [This is a message, not a warning or an error] Warning message: In DESeqDataSet(se, design = design, ignoreRank) : some variables in design formula are characters, converting to factors

• 0 views • link
written by Lorenzo • Lorenzo

No, this is just a warning, it is not really the issue here. The issue is that you are asking DESeq2 to calculate within group variability (=dispersion) on groups of 1 samples since you have only one replicate by group (=by tissue). DESeq2 obviously can't do that that is why it throws the error The design matrix has the same number of samples and coefficients to fit, so estimation of dispersion is not possible. Unless you have replicates there is nothing you can do here to solve that issue. It is impossible to make sound differential expression analysis without replicates.

• 0 views • link
written by Carlo Yague Carlo Yague

Ok thank you for your explanation. So to perform DEA it is necessary to have at least 2 replicates for each tissue? Is there a way to consider the different tissues like different replicates? I mean I have so many tissues and in my mind I'm not able to understand why it is not possible to perform DEA on them.

• 0 views • link
written by Lorenzo • Lorenzo

You could put the tissues into groups: "brain", "immune", "whatever" you could compare groups to each other. You can't meaningfully compare a single sample to another single sample. You can technically do it, just put your library normalized counts into Excel and make ratios, but that would be exploratory only. You'd need to go back and do replicates for anyone to believe that your values are biologically relevant.

• 0 views • link
written by swbarnes2 swbarnes2

Ok thank you actually I have the same type of tissues from another individual so I can do it to determine the biological relevance. Now I have to improve my skills on R to do it. Thank you again!

• 0 views • link
written by Lorenzo • Lorenzo

You only care about 2 genes? So you sequenced the whole transcriptome? That seems like a waste of time. qPCR would have been a whole lot easier. But as others have said, if you have no replicates, you really can't get a good answer. There is no algorithmic trickery that can make up for the fact that you have zero idea of how variable expression of your genes are within your different tissues, which makes it impossible to be sure if any differences you see are real.

• 0 views • link
written by swbarnes2 swbarnes2

No it is a part of a big project and considering the fact that we have the reads coming from the transcriptome we decided to do that.

• 0 views • link
written by Lorenzo • Lorenzo

1 answer

Yes, it is an appropriate option to use DESeq2 in your case. Just note that DESeq2 will take into account all the genes to (1) estimate dispersion and (2) normalize your libraries using the median of ratios method. So even if in the end you only look at two genes, all genes are used in the analysis, which is good because it make it more robust/powerful.

There are quite a few DESeq2 examples/pipelines out there. For instance, this one is well explained and complete.

• 0 views • link
written by Carlo Yague Carlo Yague

Log in to answer this question.

More posts like this