What can you learn about your samples just by looking at normalized read counts?

Question

[sandKings]

Hi everyone,

My study involves multiple cell lineages which were isolated from the mammary gland and sorted based on surface markers. The FACS data 100% agrees with published literature and what we expected to see post-treatment.

I isolated RNA from all these cell lineages and treatment groups and sent it for sequencing.

I recently finished the RNASeq analysis where I modified a script from another lab and made it work. I showed my script and outputs to several people who know better and they all agreed that the script is working flawlessly. I followed up with a GSEA analysis which looks ok to me.

The problem is, that the DEG list and GSEA analysis doesn't show all of the genes/pathways that we expected to see.

The immediate question was if my sorted population was contaminated or and/or the samples got mixed up/flipped at the sequencing core.

The sequencing core obviously denied the chances of samples flipping or being mislabeled and I want to believe them.

I suggested that, if the purity of cell lineage is a concern, then let me run a qRTPCR for the lineage-specific genes. Looking at fold change values of the KRT5, KRT14, KRT18, ESR1 etc would tell me if they are basal or luminal population.

To this, some of the people responded that why don't I look at the normalized gene counts and see the expression level of my lineage-specific markers?

My concern is, can I really do that? Can I look at the normalized read counts and say that these samples are basal cells because they have a higher count or reads of KRT5 and these are luminal cells because they have a higher count of KRT18?

Obviously, I know very little about the statistics that goes into creating all these data files but I was under the impression that the normalized read counts don't give you enough information unless you put it through an appropriate statistical analysis (like DESeq or EdgeR)

I'd appreciate any feedback or a link to any text that might help.

[ATpoint]

I would go for a principal component analysis and see if the samples cluster in the expected fashion based on your FACS data. If you need a nice tutorial on how to get going with PCA, start with this fancy tutorial from Mr. Blighe ;-)

[Benn]

Yes, PCA would be a good approach. Another similar method would be an MDS plot from edgeR/limma package using raw read counts. If there is a batch effect in your set, you can also remove that first in edgeR and then do the MDS plot with the batch removed values.

[sandKings]

That's what I want to do but I don't understand their assumption that normalized read counts of a gene should tell me if my cells are 'pure'. Wouldn't qRTPCR be more appropriate to establish the purity.

[Benn]

So your real question is whether you can use normalized counts of some lineage specific marker genes instead of RT-qPCR values? In theory yes, the normalized values are expression values, just like you get from qPCR. However, keep in mind that RNAseq is high-through put data, which means it is not as reliable as (a well designed) qPCR. Furthermore, you can't use the normalized counts for normal statistics (such as t-test), but you need ebayes statistics for that (which is done with limma or edgeR, etc.). Last point to keep in mind, how reliable are your markers? Do they act as a switch: in one line very high expressed, while in the other 0 expression? Usually not, so you'll have to find a cut off, how? And if you have contamination, what to expect from the values? Half? Very difficult questions and in my opinion not the problem for bioinformaticians. Get the wet lab work reliable first before you start sequencing is my advice.

[sandKings]

Thank you so much! This was very helpful.

Answer 1

[Kevin Blighe]

Sincere apologies to my colleagues but I am going to butt-in and provide an answer here:

Performing PCA will indeed show you the relationships between your samples in an unbiased fashion, but it is important to remember that the PCA transformation is fundamentally based on covariance. In addition, if you are coming from RNA-seq counts, your PCA transformations should be performed on the regularised log of variance-stabilised counts. Performing it on the normalised unlogged counts will be somewhat misleading because the data distribution will be on the negative binomial.

If you do PCA and then plot the values for PC1 versus PC2, you may see a 'natural' separation between 2 groups of samples along either of these (or other) PCs. What you can do then is check which genes are contributing most to the segregation of your samples by looking at their loadings to whatever is the PC of interest. I've written more about this here:

• A: PCA in a RNA seq analysis
• A: PCA plot from read count matrix from RNA-Seq

Realistically, what your colleagues may have meant is to look at the Z-scores of your genes' expression values. By transforming your data to the Z-scale, you can generally infer whether a gene is more or less expressed in a sample by only looking at genes whose Z-scores are greater than +2 or less than -2, respective. A Z-score of 2 relates to 2 standard deviations away from the mean expression value. You can also use higher numbers as cut-offs, of course. Once you identify the top genes for each sample, you may then be able to infer to which tissue they belong. Remember that tissues are also definitively defined by what they do not express, though.

Solely by looking at the normalised counts, you can still extract some meaningful data. If you do RNA-seq on blood, haemoglobin genes will always have the highest normalised counts. Mitochondrial genes are usually always near the top, too.

MDS is also a good idea, and should show similar results as PCA.

Hope that this helps Kevin