Validation of RNA-Seq knockout

In this example we're going to look at a publicly avaliable dataset in which two genes (Yap1 and Taz) have been knocked out. The data comes from GEO study GSE79115.

The data was processed through the standard clusterflow fastq_hisat2 pipeline which did trimming with TrimGalore and mapping with hisat2. The quality of the data all looked fine and the species screen showed that it was indeed mouse, as expected.

I've imported the data into SeqMonk and cleaned up the sample names a bit, I've also created replicate sets for the two biological conditions (WT and cKO), so it's all ready to go.

Project Summary

Basic Project Info

Data Sets

Replicate Sets

Basic QC

Visual Inspection

The data looks pretty good on a basic inspection. The reads are nicely restricted to exons and cover the whole length of the transcripts. Oddly the library is non-directional. The paper says that the library prep used the Illumina Tru-Seq kit, which should give a directional library, so something a bit odd is going on.

There is no obvious DNA contamination, but there is maybe a suggestion of a small amount of PCR duplication since isolated reads often appear as duplets.

RNA-Seq QC Plot

This all looks good. The library is very clean with almost all reads being in exons, but it was polyA selected so this is expected. We see no rRNA, and only a small amount of MT. There is a bit of variation in the library sizes, and this does seem to link to the condition, with all of the WT samples having more, but it's still less than 2fold different overall so I'm not too concerned. The directionality is still lacking which matches what we saw in the raw data.

Duplication

New Probe Set (411632 probes)

Feature generator using mRNA split into exons duplicates removed Over feature from 0-0

Whilst there is duplication in the samples it is in line with the local data density, so it seems that this is biological in nature rather than technical. There is nothing to worry about overall in terms of the duplication rates in the samples.

Quantitation

We'll do a standard gene level quantitation. We'll use the default log2FPM (correcting for library size but not for transcript length) since we're only interested in comparing the same transcript between samples - we don't care about comapring different transcripts in the same sample.

New Probe Set (35266 probes)

Transcript features over mRNA

Probes Quantitated

RNA-Seq pipeline quantitation on merged transcripts counting reads over exons. Log transformed. Assuming a Non-strand specific library

Normalisation

The default normalisation looks pretty good, but it might improve slightly if we add in size factor normalisation, so let's try that. We'll remove all of the obviously unexpressed genes first.

New Probe List: Expressed somewhere (19140 probes)

Filter on probes in All Probes where at least 1 of cKO-1 , cKO-2 , cKO-3 , WT-1 , WT-2 , WT-3 had a value above 50.0. Quantitation was RNA-Seq pipeline quantitation on merged transcripts counting reads over exons. Log transformed. Assuming a Non-strand specific library

Probes Quantitated

RNA-Seq pipeline quantitation on merged transcripts counting reads over exons. Log transformed. Assuming a Non-strand specific library transformed by size factor normalisation using Add using the Expressed somewhere (19140) probe list

It's slightly nicer with the additional normalisation so we'll keep that.

Exploration

Clustering

We can start by using PCA to see how the samples cluster.

The samples form two obvious clusters and these split by the knockout status on PC1 so this all looks good. The cKO-1 sample seems a little divergent on PC2 so that might need looking at.

Comparison

There are some very obvious differences between the samples, but the structure of the differences looks odd. There is a pretty obvious group of genes which appear to be very consistently upregulated in the knockout sample, all changing by a log2 diff of around 2.5 over a wide range of expression levels. This doesn't look normal and suggests that there may be an additional effect on these samples beyond the knockout, possibly a copy number change? This asymmetry in the changes probably explains why the size factor normalisation produced some benefit.

There are some genes which show quite large expression decreases in the cKO samples, but nothing looks like it really moves from being highly expressed to being absent, so we should check on the knockouts.

Knockout validation

Since the cKO samples are supposed to be a double knockout of Yap1 and Taz we should look at those genes specifically.

New Probe List: Yap1 (1 probes)

Probes from Expressed somewhere whose name matches any of Yap1 after stripping probe generator suffixes

New Probe List: Taz (1 probes)

Probes from Expressed somewhere whose name matches any of Taz after stripping probe generator suffixes

That doesn't look good. Neither of the genes which are supposed to be knocked out actually appear to show any real change in the cKO group. Depending on how the knockout was done this might not be an issue. It it's very tightly targetted then it could be that the transcripts are still being produced but in a non-functional form.

I went and looked up the paper for the original article where the Yap1/Taz knockout was created ((PMC3605093) and in both cases they used a Cre/Lox system to completely excise the first complete coding exon (exon 2). It's possible that after deletion the transcript splices around the missing exon and goes on, but we should check.

The second exon of Yap1 certainly appears to still be there, and at a similar level to the first exon. The total number of reads is fewer, but that's in line with the lower coverage overall in the cKO samples (which we saw in the QC report).

Likewise the second Taz exon also appears to be present at levels which are consistent with the exons around it. Given these observations I would be very cautious of proceeding with the analysis without going back to validate the exact knockouts in more detail, because on the face of it it looks like the expected effects are absent from these samples.

Differences examination

Despite the lack of expected knockouts in the samples we did definitely see some consistent expression differences between the two groups, but we also saw the oddly consistently regulated set of genes. To try to understand what's going on with these we can see if we can find any obvious connection between the genes in this set.

To isolate these genes we can do a crude extraction of them based on their fold change. I can just manually select these in a scatterplot.

New Probe List: Oddly similar changes (511 probes)

Filter on probes in Expressed somewhere where average difference when comparing WT to cKO had a value between 1.8 and 3.0. Quantitation was RNA-Seq pipeline quantitation on merged transcripts counting reads over exons. Log transformed. Assuming a Non-strand specific library transformed by size factor normalisation using Add using the Expressed somewhere (19140) probe list

The most obvious thing to connect these points would be phyiscal position in the genome, so let's see where they are:

Meh - they're all over the place so bang goes that theory. Let's have a look at what the actual genes are instead. We'll just take the ones which are more highly expressed to limit how many we have to deal with.

New Probe List: Highly expressed (108 probes)

Filter on probes in Oddly similar changes where exactly 1 of cKO had a value above 3.0. Quantitation was RNA-Seq pipeline quantitation on merged transcripts counting reads over exons. Log transformed. Assuming a Non-strand specific library transformed by size factor normalisation using Add using the Expressed somewhere (19140) probe list

This looks pretty telling. The genes we're looking at here all appear to be strongly related to muscles, so actin, myosin and troponin as well as mitochondrial genes which would be increased in muscle cells. I put the list through gprofiler and get enormously significant hits to muscle pathways (eg Reactome muscle contraction at p=3x10^-37.

Given that these samples are supposed to be Schwann cells, it doesn't make sense for them to be expressing these kinds of muscle markers at this level. It would seem to be much more plausible that what has actually happened is that the samples have been contaminated to varying degrees with muscle cells from the surrounding tissue, and that the amount of contamination varies between the two groups.

If this supposition is correct then it's also likely that the amount of contamination will vary within each group since it will not relate to the genotype but will be a function of the dissection or sample preparation. We can try to assess this by looking at whether the muscle related genes are unusually variable within the WT or cKO groups. We can do this with a variation plot, plotting STDEV against expression level for each of the two groups.

It is clear that in both the WT and cKO groups that these muscle markers are unusually variable, suggesting that an external contamination is a viable explanation. We can also clearly see the difference in expression level between the two groups which led us to look at these genes in the first place.

Conclusion

An examination of these samples would ideally have shown us that they contained good quality data and exhibited constant effects of the type we expect given the experimental description. Unfortunately an examination of the data instead raises serious concerns about whether these samples are indeed valid and we should probably not proceed with this analysis until the issues raised are sorted out.

Specifically there are 3 major concerns:

1. The nature of the read data looks wrong, in that this was supposed to be an Illumina Tru-Seq library, but the reads coming from it are non-directional.

2. The cKO samples are supposed to have a knockout of the Taz and Yap1 genes and the data shows no evidence of this change. Since it's a conditional knockout then maybe there is more subtelty to the design than a simple reading of the conditions suggests, but this needs to be determined since a lack of the expected effect would mean that any further work on the data would be pointless.

3. It would appear that all of the samples may be contaminated with muscle cells, and unfortuanately that the amount of contamination is not consistent between the two groups, but is systematically higher in cKO than WT. Given the conflation of the contamination and the condition it is going to be very difficult to specifically attribute any changes seen to the knockouts rather than the contamination.