Modelling RNASeq like a car in traffic: using hybrid qualitative and quantitative factors in Sleuth to find elusive significant contrasts

tl;dr If you are using gene knockouts or knockdowns (RNAi) as experimental factors in RNASeq experiments, sensitivity to detect differential expression for ALL genes can be improved by simultaneously modelling both the qualitative (knockout -/+) and quantitative (transcripts per million) values for the gene(s) being manipulated. When adding qualitative factors, optimize the regression model by minimizing the residual sum of squares of the model fit. In the example below, we went from zero mRNA transcripts with significant contrasts for factor interactions (e.g. shRNA gene knockdown+treatment) to hundreds by using this hybrid factor modelling.

The best analogy that I can think of is to imagine that you are trying to determine several metrics of a motor vehicle (such as the fuel consumption rate, RPMs, oil temperature) of a car at any given moment, given only snapshots of the motor's speed (quantitative), and the type of road the vehicle is on (highway/city, qualitative). Either known factor by itself isn't a great indicator of the desired metrics since you could be stuck in rush hour traffic on the highway, or racing stop-and-go in the city.  Together, they provide a more predictive context for the other metrics (i.e. implicitly modelling the gear the car is in, which affects fuel economy and RPMs for given speeds). In the same way, the reaction of other genes to the levels of p53 and the shRNA are likely to be better implicitly modelled by knowing if there is a knockout/down (qualitative) AND what the level of that muting is (quantitative). 

____________

Recently our centre sequenced and analyzed a three factor RNASeq experiment done in triplicate, meaning there were 24 samples.  The first factor was two cell lines (one a p53 knockout, the other wild type, i.e. minus and plus), the second factor was a shRNA knockdown of a gene of interest or GFP, and the third was a drug treatment. At first, I did a straightforward Kallisto/Sleuth analysis of the factors and their interactions. There were a bunch of p53 knockout-specific changes, and shRNA-specific changes, but not really anything else (including no treatment responses).  It was time to break out the the sample distance matrix analysis using an orthogonal count method.  Here's the matrix:

There are a few things to point out here.  First, there is a batch effect, which is evident by the grouping (clustering tree on the left margin) of samples in bottom of the matrix by sample extraction date (D1, D2, D3). So I added the date to the model, as well as its interaction with other factors (this has worked out well previously).  So my metadata table looks like this:


sample path p53 sh treatment date
HCT116_p53_minus_shATM_DMSO_Dec_1 HCT116_p53_minus_shATM_DMSO_Dec_1.kallisto minus shATM ctl D2
HCT116_p53_minus_shATM_DMSO_Dec_2 HCT116_p53_minus_shATM_DMSO_Dec_2.kallisto minus shATM ctl D3
HCT116_p53_minus_shATM_DMSO_Nov_30 HCT116_p53_minus_shATM_DMSO_Nov_30.kallisto minus shATM ctl D1
HCT116_p53_minus_shATM_PARPi_Dec_1 HCT116_p53_minus_shATM_PARPi_Dec_1.kallisto minus shATM exp D2
HCT116_p53_minus_shATM_PARPi_Dec_2 HCT116_p53_minus_shATM_PARPi_Dec_2.kallisto minus shATM exp D3
HCT116_p53_minus_shATM_PARPi_Nov_30 HCT116_p53_minus_shATM_PARPi_Nov_30.kallisto minus shATM exp D1
HCT116_p53_minus_shGFP_DMSO_Dec_1 HCT116_p53_minus_shGFP_DMSO_Dec_1.kallisto minus GFP ctl D2
HCT116_p53_minus_shGFP_DMSO_Dec_2 HCT116_p53_minus_shGFP_DMSO_Dec_2.kallisto minus GFP ctl D3
HCT116_p53_minus_shGFP_DMSO_Nov_30 HCT116_p53_minus_shGFP_DMSO_Nov_30.kallisto minus GFP ctl D1
HCT116_p53_minus_shGFP_PARPi_Dec_1 HCT116_p53_minus_shGFP_PARPi_Dec_1.kallisto minus GFP exp D2
HCT116_p53_minus_shGFP_PARPi_Dec_2 HCT116_p53_minus_shGFP_PARPi_Dec_2.kallisto minus GFP exp D3
HCT116_p53_minus_shGFP_PARPi_Nov_30 HCT116_p53_minus_shGFP_PARPi_Nov_30.kallisto minus GFP exp D1
HCT116_p53_plus_shATM_DMSO_Dec_1 HCT116_p53_plus_shATM_DMSO_Dec_1.kallisto plus shATM ctl D2
HCT116_p53_plus_shATM_DMSO_Dec_2 HCT116_p53_plus_shATM_DMSO_Dec_2.kallisto plus shATM ctl D3
HCT116_p53_plus_shATM_DMSO_Nov_30 HCT116_p53_plus_shATM_DMSO_Nov_30.kallisto plus shATM ctl D1
HCT116_p53_plus_shATM_PARPi_Dec_1 HCT116_p53_plus_shATM_PARPi_Dec_1.kallisto plus shATM exp D2
HCT116_p53_plus_shATM_PARPi_Dec_2 HCT116_p53_plus_shATM_PARPi_Dec_2.kallisto plus shATM exp D3
HCT116_p53_plus_shATM_PARPi_Nov_30 HCT116_p53_plus_shATM_PARPi_Nov_30.kallisto plus shATM exp D1
HCT116_p53_plus_shGFP_DMSO_Dec_1 HCT116_p53_plus_shGFP_DMSO_Dec_1.kallisto plus GFP ctl D2
HCT116_p53_plus_shGFP_DMSO_Dec_2 HCT116_p53_plus_shGFP_DMSO_Dec_2.kallisto plus GFP ctl D3
HCT116_p53_plus_shGFP_DMSO_Nov_30 HCT116_p53_plus_shGFP_DMSO_Nov_30.kallisto plus GFP ctl D1
HCT116_p53_plus_shGFP_PARPi_Dec_1 HCT116_p53_plus_shGFP_PARPi_Dec_1.kallisto plus GFP exp D2
HCT116_p53_plus_shGFP_PARPi_Dec_2 HCT116_p53_plus_shGFP_PARPi_Dec_2.kallisto plus GFP exp D3
HCT116_p53_plus_shGFP_PARPi_Nov_30 HCT116_p53_plus_shGFP_PARPi_Nov_30.kallisto plus GFP exp D1

And the analysis looks like so:


library(sleuth)
meta_qual <- read.table("meta_qual.tab", header=TRUE)
meta_qual$path <- as.character(meta$path)
so_qual <- sleuth_prep(meta_qual, ~p53*sh+treatment+p53:treatment+p53:date+sh:date+sh:treatment+treatment:date+date)

reading in kallisto results

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

normalizing est_counts

29870 targets passed the filter

normalizing tpm

merging in metadata

normalizing bootstrap samples

summarizing bootstraps

> so_qual <- sleuth_fit(so_qual)

fitting measurement error models

shrinkage estimation

Adding missing grouping variables: `x_group`

computing variance of betas

I'll exclude the Sleuth messages from hereon in for brevity. After running through a typical likelihood ratio test workflow for the factors, let's see how many genes are differentially expressed in each factor contrast (3) and their interactions (also 3, for a total of 6).  We modelled the date effect but aren't interested in reporting them, just accounting for them to improve the factors of interest.

> length(lrt_p53.sig_ids)
[1] 19472
> length(lrt_sh.sig_ids)
[1] 17970
> length(lrt_treatment.sig_ids)
[1] 28
> length(lrt_p53_sh_int.sig_ids)
[1] 9579
> length(lrt_p53_treatment_int.sig_ids)
[1] 4
> length(lrt_sh_treatment_int.sig_ids)
[1] 0

Okay, this is looking somewhat better.  We do have some treatment-specific changes now, and a bunch of p53:sh interaction changes.  Inspecting these, about half are mirror images of the p53 changes (gene +X change in p53 plus, -X in shATM), which suggests that we have confounding factors we haven't modelled yet. A closer inspection of the distance matrix produced earlier shows that there are four p53 minus samples that cluster quite closely with the p53 plus samples, including a mix of treated and untreated, and the rest of the shATM knockdowns are peppered semi-randomly in the p53 minus part of the tree. What are we missing?

To gauge how well we've modelled the driving factors of the gene expression, we can look at the residual sum of squares for the regression.  Minimization is optimization for this sum. Where can we find this?  A little digging shows that the so object has a fits member, which itself has named regression models.  By default, the model we first fitted is called "full". Each regression model contains a summary data frame.  Let's check it out the first row of that data frame:

> so_qual$fits$full$summary[1,]

           x_group      rss sigma_sq sigma_q_sq mean_obs  var_obs target_id[...]

1 (-0.000966,0.01] 29.05466 2.021473   1.206822 1.246232 3.421352 NM_000195[...]

That second column, "rss" is the Residual Sum of Squares we're after. Let's look at the sum across all the transcripts modelled:

> sum(so_qual$fits$full$summary[,2])

[1] 232757.7

If we're going to improve the regression model, which leads to better statistical power to pick up the interaction terms, we'll want to get a total RSS of less than 232758. RNAi knockdown is infamous for having highly variable efficiency, so maybe we should not be using a binary on/off factor to model it, but rather model using the actual expression value since we have that from the RNAseq data itself.  My Kallisto reference file uses RefSeq MRNA models, so I look up the main NM_####### ID for ATM, which can be found by clicking the RefSeq mRNA link on the right hand sidebar of the ATM RefGene page.  It's NM_000051.  The transcripts per million information is readily available as the last column in the Kallisto abundance files in each sample's output directory.

-bash-4.2$ grep NM_000051 *.kallisto/abundance.tsv
HCT116_p53_minus_shATM_DMSO_Dec_1.kallisto/abundance.tsv:NM_000051 13147 12968 1882.6 4.93759
HCT116_p53_minus_shATM_DMSO_Dec_2.kallisto/abundance.tsv:NM_000051 13147 12968 2618.66 5.9662
HCT116_p53_minus_shATM_DMSO_Nov_30.kallisto/abundance.tsv:NM_000051 13147 12968 2392.05 4.00764
HCT116_p53_minus_shATM_PARPi_Dec_1.kallisto/abundance.tsv:NM_000051 13147 12968 2212.67 4.54057
HCT116_p53_minus_shATM_PARPi_Dec_2.kallisto/abundance.tsv:NM_000051 13147 12968 2345.31 5.07381
HCT116_p53_minus_shATM_PARPi_Nov_30.kallisto/abundance.tsv:NM_000051 13147 12968 3151.77 5.03178
HCT116_p53_minus_shGFP_DMSO_Dec_1.kallisto/abundance.tsv:NM_000051 13147 12968 3851.38 7.11942
HCT116_p53_minus_shGFP_DMSO_Dec_2.kallisto/abundance.tsv:NM_000051 13147 12968 1981.88 4.0778
HCT116_p53_minus_shGFP_DMSO_Nov_30.kallisto/abundance.tsv:NM_000051 13147 12968 2621.98 8.19509
HCT116_p53_minus_shGFP_PARPi_Dec_1.kallisto/abundance.tsv:NM_000051 13147 12968 2458.8 3.66504
HCT116_p53_minus_shGFP_PARPi_Dec_2.kallisto/abundance.tsv:NM_000051 13147 12968 2710.64 4.47078
HCT116_p53_minus_shGFP_PARPi_Nov_30.kallisto/abundance.tsv:NM_000051 13147 12968 2243.22 5.37964
HCT116_p53_plus_shATM_DMSO_Dec_1.kallisto/abundance.tsv:NM_000051 13147 12968 2441.1 4.88496
HCT116_p53_plus_shATM_DMSO_Dec_2.kallisto/abundance.tsv:NM_000051 13147 12968 3736.36 4.75133
HCT116_p53_plus_shATM_DMSO_Nov_30.kallisto/abundance.tsv:NM_000051 13147 12968 3259.82 6.86958
HCT116_p53_plus_shATM_PARPi_Dec_1.kallisto/abundance.tsv:NM_000051 13147 12968 3175.35 4.85537
HCT116_p53_plus_shATM_PARPi_Dec_2.kallisto/abundance.tsv:NM_000051 13147 12968 1790.03 4.44556
HCT116_p53_plus_shATM_PARPi_Nov_30.kallisto/abundance.tsv:NM_000051 13147 12968 3625.5 6.01011
HCT116_p53_plus_shGFP_DMSO_Dec_1.kallisto/abundance.tsv:NM_000051 13147 12968 1662.14 4.27169
HCT116_p53_plus_shGFP_DMSO_Dec_2.kallisto/abundance.tsv:NM_000051 13147 12968 2445.64 5.1909
HCT116_p53_plus_shGFP_DMSO_Nov_30.kallisto/abundance.tsv:NM_000051 13147 12968 2185.2 4.63558
HCT116_p53_plus_shGFP_PARPi_Dec_1.kallisto/abundance.tsv:NM_000051 13147 12968 1640.12 4.15573
HCT116_p53_plus_shGFP_PARPi_Dec_2.kallisto/abundance.tsv:NM_000051 13147 12968 2548.95 4.62631
HCT116_p53_plus_shGFP_PARPi_Nov_30.kallisto/abundance.tsv:NM_000051 13147 12968 1864.82 4.69346

I did the same for the p53 mRNA level (NM_000546), since there is no guarantee that p53 is evenly expressed in each of the samples. Now my sample metadata file looks like (meta_quant.tab):

sample path p53 sh treatment date
HCT116_p53_minus_shATM_DMSO_Dec_1 HCT116_p53_minus_shATM_DMSO_Dec_1.kallisto 1.07691 4.93759 ctl D2
HCT116_p53_minus_shATM_DMSO_Dec_2 HCT116_p53_minus_shATM_DMSO_Dec_2.kallisto 4.67624 5.9662 ctl D3
HCT116_p53_minus_shATM_DMSO_Nov_30 HCT116_p53_minus_shATM_DMSO_Nov_30.kallisto 1.17641 4.00764 ctl D1
HCT116_p53_minus_shATM_PARPi_Dec_1 HCT116_p53_minus_shATM_PARPi_Dec_1.kallisto 5.94978 4.54057 exp D2
HCT116_p53_minus_shATM_PARPi_Dec_2 HCT116_p53_minus_shATM_PARPi_Dec_2.kallisto 5.98334 5.07381 exp D3
HCT116_p53_minus_shATM_PARPi_Nov_30 HCT116_p53_minus_shATM_PARPi_Nov_30.kallisto 7.44437 5.03178 exp D1
HCT116_p53_minus_shGFP_DMSO_Dec_1 HCT116_p53_minus_shGFP_DMSO_Dec_1.kallisto 0.53949 7.11942 ctl D2
HCT116_p53_minus_shGFP_DMSO_Dec_2 HCT116_p53_minus_shGFP_DMSO_Dec_2.kallisto 8.92308 4.0778 ctl D3
HCT116_p53_minus_shGFP_DMSO_Nov_30 HCT116_p53_minus_shGFP_DMSO_Nov_30.kallisto 5.6365 8.19509 ctl D1
HCT116_p53_minus_shGFP_PARPi_Dec_1 HCT116_p53_minus_shGFP_PARPi_Dec_1.kallisto 10.7367 3.66504 exp D2
HCT116_p53_minus_shGFP_PARPi_Dec_2 HCT116_p53_minus_shGFP_PARPi_Dec_2.kallisto 10.1989 4.47078 exp D3
HCT116_p53_minus_shGFP_PARPi_Nov_30 HCT116_p53_minus_shGFP_PARPi_Nov_30.kallisto 16.878 5.37964 exp D1
HCT116_p53_plus_shATM_DMSO_Dec_1 HCT116_p53_plus_shATM_DMSO_Dec_1.kallisto 0.10751 4.88496 ctl D2
HCT116_p53_plus_shATM_DMSO_Dec_2 HCT116_p53_plus_shATM_DMSO_Dec_2.kallisto 0.55286 4.75133 ctl D3
HCT116_p53_plus_shATM_DMSO_Nov_30 HCT116_p53_plus_shATM_DMSO_Nov_30.kallisto 0.38625 6.86958 ctl D1
HCT116_p53_plus_shATM_PARPi_Dec_1 HCT116_p53_plus_shATM_PARPi_Dec_1.kallisto 0.77159 4.85537 exp D2
HCT116_p53_plus_shATM_PARPi_Dec_2 HCT116_p53_plus_shATM_PARPi_Dec_2.kallisto 0.82322 4.44556 exp D3
HCT116_p53_plus_shATM_PARPi_Nov_30 HCT116_p53_plus_shATM_PARPi_Nov_30.kallisto 0.59059 6.01011 exp D1
HCT116_p53_plus_shGFP_DMSO_Dec_1 HCT116_p53_plus_shGFP_DMSO_Dec_1.kallisto 0.52125 4.27169 ctl D2
HCT116_p53_plus_shGFP_DMSO_Dec_2 HCT116_p53_plus_shGFP_DMSO_Dec_2.kallisto 0.72088 5.1909 ctl D3
HCT116_p53_plus_shGFP_DMSO_Nov_30 HCT116_p53_plus_shGFP_DMSO_Nov_30.kallisto 0.39565 4.63558 ctl D1
HCT116_p53_plus_shGFP_PARPi_Dec_1 HCT116_p53_plus_shGFP_PARPi_Dec_1.kallisto 0.71229 4.15573 exp D2
HCT116_p53_plus_shGFP_PARPi_Dec_2 HCT116_p53_plus_shGFP_PARPi_Dec_2.kallisto 0.68656 4.62631 exp D3
HCT116_p53_plus_shGFP_PARPi_Nov_30 HCT116_p53_plus_shGFP_PARPi_Nov_30.kallisto 0.65503 4.69346 exp D1

Note that the p53 expression values are all over the place in the nominally p53 "minus" samples, which gives us some hope that quantitative modelling will help. Let's rerun the analysis, except this time the regression model will treat p53 and sh as dosage effects, which is the default for any numeric factor column in R.

> library(sleuth)
> meta_quant <- read.table("meta_quant.tab", header=TRUE)
> meta_quant$path <- as.character(meta$path)

> so_quant <- sleuth_prep(meta_quant, 
~p53*sh+treatment+p53:treatment+sh:treatment+p53:date+sh:date+treatment:date+date)

> so_quant <- sleuth_fit(so_quant)

Let's check if we improved the regression model by using the quantitative factors instead of qualitative.

> sum(so_quant$fits$full$summary[,2])

[1] 325008.8

Ouch, that's a lot more than our original 232757.7 RSS!  At this point I tried a few mathematical manipulations of the expression values to see if that'd help.  After all, maybe the dosage response is logarithmic rather than linear?  No need to run sleuth_prep() again, we can just try as many extra models as we like using sleuth_fit().

> so_quant <- sleuth_fit(so_quant, ~I(log(p53))*sh+treatment+I(log(p53)):treatment+I(log(p53)):date+sh:date+sh:treatment+treatment:date+date, "full_log_p53")

You'll note that I replaced p53 with I(log(p53)).  The I is a function that Inhibits Interpretation, allowing the logarithm function to be passed to the model rather than evaluated in the command call itself.  Did that help?

> sum(so_quant$fits$full_log_p53$summary[,2])

[1] 300397.4

A bit, but not close to the original yet.  Since their are a lot of values for p53 range from 0.108 to 1, we're getting a much big spread (-2.2, 0) in the logarithm of small values than in large values, so let's round up to 1, which will effectively make samples with <1 have a dosage value of 0 (because log(1) = 0). Let's also invert the sh quantity, which gives greater weight to small values, i.e. better knockdown.

> so_quant <- sleuth_fit(so, ~I(log(ceiling(p53)))*I(1/sh)+treatment+I(log(ceiling(p53))):treatment+I(log(ceiling(p53))):date+I(1/sh):date+I(1/sh):treatment+treatment:date+date, "full_p53_log_and_sh_reciprocal")

> sum(so$fits$full_p53_log_and_sh_reciprocal$summary[,2])

[1] 274487.9

Better, but still not as good as the original qualitative-only model.  This is reflected in the LRT results if we bother to do the normal downstream workflow: 

> length(lrt_p53.sig_ids)

[1] 13282

> length(lrt_sh.sig_ids)

[1] 14224

> length(lrt_treatment.sig_ids)

[1] 4

> length(lrt_treatment_sh_int.sig_ids)

[1] 2

> length(lrt_treatment_p53_int.sig_ids)

[1] 0

> length(lrt_sh_p53_int.sig_ids)

[1] 1

Various other manipulation didn't help much.  It then struck me that modelling both the qualitative and quantitative values simultaneously could help.  Let's combine the metadata and see, first modelling just the main p53 effect seen in the distance matrix.

> meta_hybrid <- meta_qual

> meta_hybrid$p53_quant <- meta_quant$p53

> meta_hybrid$sh_quant <- meta_quant$sh

> so_hybrid <- sleuth_prep(meta_hybrid, ~p53*sh+treatment+p53:treatment+p53:date+sh:date+sh:treatment+treatment:date+date+p53_quant)

> so_hybrid <- sleuth_fit(so_hybrid)

> sum(so_hybrid$fits$full$summary[,2])

[1] 193880.4

Woohoo!  That around a 20% drop in the residual sum of squares, and already we're seeing that the additional p53_quant factor is removing enough systematic bias to start making the interaction effects of the 3 factors significant.

> so_hybrid <- sleuth_fit(so_hybrid, ~p53*sh+treatment+p53:treatment+p53:date+sh:date+sh:treatment+treatment:date+date+I(log(ceiling(p53_quant))), "hybrid_log_p53")
[a bunch of boring code goes here, see below...]
> length(lrt_p53.sig_ids)
[1] 22346
> length(lrt_sh.sig_ids)
[1] 2634
> length(lrt_treatment.sig_ids)
[1] 9
> length(lrt_sh_treatment_int.sig_ids)
[1] 64
> length(lrt_p53_treatment_int.sig_ids)
[1] 17
> length(lrt_sh_p53_int.sig_ids)
[1] 2511

Let's add in the sh quantitative effect.

> so_hybrid <- sleuth_fit(so_hybrid, ~p53*sh+treatment+p53:treatment+p53:date+sh:date+sh:treatment+treatment:date+date+p53_quant+sh_quant, "full_sh")

> sum(so_hybrid$fits$full_sh$summary[,2])

[1] 150439.9

Excellent! We've managed to reduce the residuals by almost a third (~23K to ~15K).  Further mathematical manipulations of the quantitative factors didn't help, so let's assume this is as good as it gets continue with a standard LRT differential expression workflow from here.

> so_hybrid <- sleuth_fit(so_hybrid, ~sh+treatment+sh:date+sh:treatment+treatment:date+date+sh_quant, "no_p53")

> so_hybrid <- sleuth_fit(so_hybrid, ~p53+treatment+p53:date+p53:treatment+treatment:date+date+p53_quant, "no_sh")

> so_hybrid <- sleuth_fit(so_hybrid, ~p53*sh+p53:date+sh:date+date+p53_quant+sh_quant, "no_treatment")

> so_hybrid <- sleuth_fit(so_hybrid, ~p53*sh+treatment+p53:date+sh:date+p53:treatment+treatment:date+date+p53_quant+sh_quant, "no_sh_treatment_int")

> so_hybrid <- sleuth_fit(so_hybrid, ~p53*sh+treatment+p53:date+sh:date+sh:treatment+treatment:date+date+p53_quant+sh_quant, "no_p53_treatment_int")

> so_hybrid <- sleuth_fit(so_hybrid, ~p53+sh+treatment+p53:date+sh:date+p53:treatment+sh:treatment+treatment:date+date+p53_quant+sh_quant, "no_sh_p53_int")

> so_hybrid <- sleuth_lrt(so_hybrid, 'no_p53', 'full_sh')

> so_hybrid <- sleuth_lrt(so_hybrid, 'no_sh', 'full_sh')

> so_hybrid <- sleuth_lrt(so_hybrid, 'no_treatment', 'full_sh')

> so_hybrid <- sleuth_lrt(so_hybrid, 'no_p53_treatment_int' , 'full_sh')

> so_hybrid <- sleuth_lrt(so_hybrid, 'no_sh_treatment_int', 'full_sh')

> so_hybrid <- sleuth_lrt(so_hybrid, 'no_sh_p53_int', 'full_sh')
> lrt_p53 <- sleuth_results(so_hybrid, 'no_p53:full_sh', test_type = 'lrt')

> lrt_sh <- sleuth_results(so_hybrid, 'no_sh:full_sh', test_type = 'lrt')

> lrt_treatment <- sleuth_results(so_hybrid, 'no_treatment:full_sh', test_type = 'lrt')

> lrt_p53_treatment_int <- sleuth_results(so_hybrid, 'no_p53_treatment_int:full_sh', test_type = 'lrt')

> lrt_sh_treatment_int <- sleuth_results(so_hybrid, 'no_sh_treatment_int:full_sh', test_type = 'lrt')

> lrt_sh_p53_int <- sleuth_results(so_hybrid, 'no_sh_p53_int:full_sh', test_type = 'lrt')

> lrt_p53.sig_ids <- lrt_p53$target_id[which(lrt_p53$qval < 0.05)]

> lrt_sh.sig_ids <- lrt_sh$target_id[which(lrt_sh$qval < 0.05)]

> lrt_treatment.sig_ids <- lrt_treatment$target_id[which(lrt_treatment$qval < 0.05)]

> lrt_sh_treatment_int.sig_ids <- lrt_sh_treatment_int$target_id[which(lrt_sh_treatment_int$qval < 0.05)]

> lrt_p53_treatment_int.sig_ids <- lrt_p53_treatment_int$target_id[which(lrt_p53_treatment_int$qval < 0.05)]

> lrt_sh_p53_int.sig_ids <- lrt_sh_p53_int$target_id[which(lrt_sh_p53_int$qval < 0.05)]

Now, let's look how this improved model has affected the factor contrast lists under the LRT test.

> length(lrt_p53.sig_ids)

[1] 23099

> length(lrt_sh.sig_ids)

[1] 13483

> length(lrt_treatment.sig_ids)

[1] 164

> length(lrt_sh_treatment_int.sig_ids)

[1] 291

> length(lrt_p53_treatment_int.sig_ids)

[1] 161

> length(lrt_sh_p53_int.sig_ids)

[1] 5515

This is a lot more in line with what we know from previous experiments and Western blots, etc.  Important to note is that we got significant changes in our original qualitative factors (p53 minus/plus and shGFP/ATM) by reducing the confounding expression "noise" with the quantitative models. Note that if you run an LRT test on a quantitative factor in Sleuth, it makes an assumption of normally distributed values for abundance in each sample, which is unlikely to be true unless you've very carefully chosen your samples. Let's go on and run the Wald tests so we can get the betas (~natural logarithm fold change), which IS something we can model for the quantitative traits on a per-quantitative-unit basis if we picked the right quantitative effect (linear, log, etc.).

> so_hybrid <- sleuth_wt(so_hybrid, 'p53plus')
> so_hybrid <- sleuth_wt(so_hybrid, 'shshATM')
> so_hybrid <- sleuth_wt(so_hybrid, 'treatmentexp')
> so_hybrid <- sleuth_wt(so_hybrid, 'p53_quant')
> so_hybrid <- sleuth_wt(so_hybrid, 'sh_quant')
> so_hybrid <- sleuth_wt(so_hybrid, 'shshATM:treatmentexp')
> so_hybrid <- sleuth_wt(so_hybrid, 'p53plus:shshATM')
> so_hybrid <- sleuth_wt(so_hybrid, 'p53plus:treatmentexp')
> wt_p53 <- sleuth_results(so_hybrid, 'p53plus')
> wt_sh <- sleuth_results(so_hybrid, 'shshATM')
> wt_treatment <- sleuth_results(so_hybrid, 'treatmentexp')
> wt_p53_quant <- sleuth_results(so_hybrid, 'p53_quant')
> wt_sh_quant <- sleuth_results(so_hybrid, 'sh_quant')
> wt_sh_treatment_int <- sleuth_results(so_hybrid, 'shshATM:treatmentexp')
> wt_p53_sh_int <- sleuth_results(so_hybrid, 'p53plus:shshATM')
> wt_p53_treatment_int <- sleuth_results(so_hybrid, 'p53plus:treatmentexp')

Now let's keep the Wald test results that are also in the LRT results, and reorder them so everyone's got the their results for the same targets (RefSeq mRNA IDs) alphabetically, then print the important result bits (ID, WT and LRT q-values, fold changes for each factor & interaction) to a file.

> lrt_ids <- unique(lrt_p53.sig_ids, lrt_sh.sig_ids, lrt_treatment.sig_ids, lrt_sh_treatment_int.sig_ids, lrt_p53_treatment_int.sig_ids, lrt_sh_p53_int.sig_ids)

> lrt_shared_p53 <- lrt_p53[lrt_p53$target_id %in% lrt_ids,]

> lrt_shared_sh <- lrt_sh[lrt_sh$target_id %in% lrt_ids,]

> lrt_shared_treatment <- lrt_treatment[lrt_treatment$target_id %in% lrt_ids,]

> lrt_shared_p53_quant <- lrt_treatment[lrt_p53_quant$target_id %in% lrt_ids,]

> lrt_shared_sh_quant <- lrt_treatment[lrt_sh_quant$target_id %in% lrt_ids,]

> lrt_shared_sh_treatment_int <- lrt_sh_treatment_int[lrt_sh_treatment_int$target_id %in% lrt_ids,]

> lrt_shared_p53_treatment_int <- lrt_sh_treatment_int[lrt_p53_treatment_int$target_id %in% lrt_ids,]

> lrt_shared_sh_p53_int <- lrt_p53_sh_int[lrt_p53_sh_int$target_id %in% lrt_ids,]

> lrt_shared_p53 <- lrt_shared_p53[sort.list(lrt_shared_p53[,1]),]

> lrt_shared_sh <- lrt_shared_sh[sort.list(lrt_shared_sh[,1]),]

> lrt_shared_treatment <- lrt_shared_treatment[sort.list(lrt_shared_treatment[,1]),]

> lrt_shared_p53_quant <- lrt_shared_p53_quant[sort.list(lrt_shared_p53_quant[,1]),]

> lrt_shared_sh_quant <- lrt_shared_sh_quant[sort.list(lrt_shared_sh_quant[,1]),]

> lrt_shared_sh_treatment_int <- lrt_shared_sh_treatment_int[sort.list(lrt_shared_sh_treatment_int[,1]),]

> lrt_shared_p53_treatment_int <- lrt_shared_p53_treatment_int[sort.list(lrt_shared_p53_treatment_int[,1]),]

> lrt_shared_sh_p53_int <- lrt_shared_sh_p53_int[sort.list(lrt_shared_sh_p53_int[,1]),]

> wt_shared_p53 <- wt_p53[wt_p53$target_id %in% lrt_ids,]

> wt_shared_sh <- wt_sh[wt_sh$target_id %in% lrt_ids,]

> wt_shared_treatment <- wt_treatment[wt_treatment$target_id %in% lrt_ids,]

> wt_shared_p53_quant <- wt_treatment[wt_p53_quant$target_id %in% lrt_ids,]

> wt_shared_sh_quant <- wt_treatment[wt_sh_quant$target_id %in% lrt_ids,]

> wt_shared_sh_treatment_int <- wt_sh_treatment_int[wt_sh_treatment_int$target_id %in% lrt_ids,]

> wt_shared_p53_treatment_int <- wt_sh_treatment_int[wt_p53_treatment_int$target_id %in% lrt_ids,]

> wt_shared_sh_p53_int <- wt_p53_sh_int[wt_p53_sh_int$target_id %in% lrt_ids,]

> wt_shared_p53 <- wt_shared_p53[sort.list(wt_shared_p53[,1]),]

> wt_shared_sh <- wt_shared_sh[sort.list(wt_shared_sh[,1]),]

> wt_shared_treatment <- wt_shared_treatment[sort.list(wt_shared_treatment[,1]),]

> wt_shared_p53_quant <- wt_shared_p53_quant[sort.list(wt_shared_p53_quant[,1]),]

> wt_shared_sh_quant <- wt_shared_sh_quant[sort.list(wt_shared_sh_quant[,1]),]

> wt_shared_sh_treatment_int <- wt_shared_sh_treatment_int[sort.list(wt_shared_sh_treatment_int[,1]),]

> wt_shared_p53_treatment_int <- wt_shared_p53_treatment_int[sort.list(wt_shared_p53_treatment_int[,1]),]

> wt_shared_sh_p53_int <- wt_shared_sh_p53_int[sort.list(wt_shared_sh_p53_int[,1]),]

> combined_result <- data.frame(lrt_shared_p53[,"target_id"],lrt_shared_p53[,"qval"],wt_shared_p53[,"b"],wt_shared_p53[,"qval"],wt_shared_p53_quant[,"b"],wt_shared_p53_quant[,"qval"],lrt_shared_sh[,"qval"],wt_shared_sh[,"b"],wt_shared_sh[,"qval"],wt_shared_sh_quant[,"b"],wt_shared_sh_quant[,"qval"],lrt_shared_treatment[,"qval"],wt_shared_treatment[,"b"],wt_shared_treatment[,"qval"],lrt_shared_sh_p53_int[,"qval"],wt_shared_sh_p53_int[,"b"],wt_shared_sh_p53_int[,"qval"],lrt_shared_p53_treatment_int[,"qval"],wt_shared_p53_treatment_int[,"b"],wt_shared_p53_treatment_int[,"qval"],lrt_shared_sh_treatment_int[,"qval"],wt_shared_sh_treatment_int[,"b"],wt_shared_sh_treatment_int[,"qval"])

> write.csv(combined_result, file="hct_hybrid_dge_lrt_q_0.05.csv")

The results are encouraging. Most of the sh:p53 interactions that were mirroring the p53 changes have gone away, and most of the LRT-significant changes in the interaction terms aren't mirrors either. We've dealt with the batch effects, and the variability of gene expression in knockdowns and knockouts pretty effectively.  You'll note that the q-values for the LRT tests are better than the Wald test in situations like this with so many factors in the regression model, because the LRT test has greater power.  Not all the genes listed as significant in the LRT are going to be real, but it's certainly a good starting point for pathway enrichment analysis, etc., and the Wald test results with qval < 0.05 can be taken as more of a sure thing.

Debugging RNASeq: sample swaps and batch effects

tl;dr If a standard differential expression analysis yields nothing for an experiment, look at the distance matrix between samples in an RNASeq experiment and see if anything looks out of place. Sometimes you can validate sample swaps when different cell lines are involved, and even strong batch effects can be modelled effectively to reveal the underlying factor effects you're expecting.

______________

There's a particular sinking feeling you get when you do a differential gene expression analysis of a carefully planned RNASeq experiment and you get no significant changes. I got that recently, when analyzing a two-factor Latin Square experiment done in triplicate.  The experiment involved two different cell lines (factor 1, with levels E545K and EV), and a treatment (factor 2, with level plus and minus).  We had good reason to believe that there should be an interaction effect (factor 3) between the one of the cell lines and the treatment.  I ran a Kallisto + Sleuth analysis with a linear model of ~cell_line+treatment+cell_line:treatment, and an experiment metadata file like so:

sample                path                           cell_line treatment 
E545K_minus_CP_Dec_1  E545K_minus_CP_Dec_1.kallisto  K         ctl      
E545K_minus_CP_Dec_2  E545K_minus_CP_Dec_2.kallisto  K         ctl       
E545K_minus_CP_Nov_30 E545K_minus_CP_Nov_30.kallisto K         ctl       
E545K_plus_CP_Dec_1   E545K_plus_CP_Dec_1.kallisto   K         exp      
E545K_plus_CP_Dec_2   E545K_plus_CP_Dec_2.kallisto   K         exp       
E545K_plus_CP_Nov_30  E545K_plus_CP_Nov_30.kallisto  K         exp      
EV_minus_CP_Dec_1     EV_minus_CP_Dec_1.kallisto     EV        ctl      
EV_minus_CP_Dec_2     EV_minus_CP_Dec_2.kallisto     EV        ctl       
EV_minus_CP_Nov_30    EV_minus_CP_Nov_30.kallisto    EV        ctl      
EV_plus_CP_Dec_1      EV_plus_CP_Dec_1.kallisto      EV        exp      
EV_plus_CP_Dec_2      EV_plus_CP_Dec_2.kallisto      EV        exp      
EV_plus_CP_Nov_30     EV_plus_CP_Nov_30.kallisto     EV        exp      

Note that I called the E545K samples "K", because I'm lazy and want to use the default behaviour of R which is to make the alphabetically first factor value the base value, so reported ratios will be K:EV. As I said before, I got nothing significant...so I decided to run a bwa + DESeq2 analysis instead with the same model. That analysis generates a distance matrix like so:

[Note: you can glean similar information from the Sleuth PCA plots in the sleuth_live interface, but I often find them too cluttered] Probably the most important thing to look at on this chart is the cluster tree on the left hand side. First, most of the EV cell line samples cluster together at the bottom, with the one exception of EV_minus_CP_Dec2 way up with the first sample, E545K_plus_CP_Dec2. A single E545K sample, E545K_minus_CP_Dec2, clusters with the EV samples.  Maybe these two are swapped?

For supporting evidence, we can look in the mapped reads for the G>A variant that's specific to cell line E545K.  There's something you can't do with a gene expression microarray! (Talk about kicking a platform when it's down).
for b in E*.bam; do
  echo $b;
  samtools mpileup -r chr3:178936091-178936091 $b 2> /dev/null;
done

E545K_minus_CP_Dec_1.bam
chr3 178936091 N 64 ggggggggggggggggggggggaggggggggggggggggggggggggggggggggggggg^Wg^Wg^Wg^Wg AEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEE/EEEEEEE/EEEEEEEEEEEEEEEEEEEE

E545K_minus_CP_Dec_2.bam
chr3 178936091 N 20 gggggggggggggggggggg EEEAEEEEEEEEEE/EE<EE

E545K_minus_CP_Nov_30.bam
chr3 178936091 N 77 ggggggggggggggggggggggggggggggggggggggggggggggggggggggaggggggggggggggggggggg^Wg AAAEEEEEEAEEA/EEEEEEE/EEEAEAEEEEEEEEEEEEAAEEEE/EEEEAAEEEEEEEEEEEEEE/E<<EAEEEE

E545K_plus_CP_Dec_1.bam
chr3 178936091 N 33 gggggggggaggggggggggggagggggggggg AAAEEEEEEE<EEEE<E/EEEEEEEEAEE/EEE

E545K_plus_CP_Dec_2.bam
chr3 178936091 N 56 ggggggggAgggggggggggggggggggggggggggggggggggggggaggggggg AEEEEEEEEEAEAEE/EEEEEEEEEEEEEAEAEEEEE/AE/EE<E/EEAEEEEE<E

E545K_plus_CP_Nov_30.bam
chr3 178936091 N 36 ggagggggggaggggggggggggggggggggggggg E6EEEE/EEEE/AEEEEEAEEEE/EEAAEEEEEE//

EV_minus_CP_Dec_1.bam
chr3 178936091 N 30 g$ggggggggggggggggggggggggggggg AAAAEEEEEEEEEE/EEEEEAEEE/EEEA<

EV_minus_CP_Dec_2.bam
chr3 178936091 N 47 gggggggggagggggggaggggggggggaggggggggggggggggg^Wg AEAEEEEEEEEEEAEEEEEEEEEE/EEAEEEEEEEEAE6EEEE/AEE

EV_minus_CP_Nov_30.bam
chr3 178936091 N 52 g$ggggggggggggggggggggggggggggggggggggggggggggggggggg AAEEEAEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEE6EEEEAEEEE<EA/

EV_plus_CP_Dec_1.bam
chr3 178936091 N 31 gggggggggggggggggggggggggggggg^Wg AEEEEEE/EEEEEEEAEAEEEEEEE<A/EEE

EV_plus_CP_Dec_2.bam
chr3 178936091 N 30 gggggggggggggggggggggggggggggg AEEEEAEEEEEEEEAEEEEEEEEEE6<AAE

EV_plus_CP_Nov_30.bam
chr3 178936091 N 34 ggggggggGggggggggggggggggggggggggg AAAEEEE/EEEEEEEEA/EE/EEEEEEEEEEEA<

I've added in the DNA base highlights manually to show which samples have the variant. We'd expect it to be the first 6 normally (all the E545K... samples). The two that we suspect are swapped (also highlighted) defy the expectation: the EV sample has G>A, the E545K sample doesn't. Had the swap been between sample in the same cell line we wouldn't be able to use this evidence, but we're lucky I guess.

Rerunning the analysis with the samples swapped yields almost nothing.

The second thing to notice about the distance matrix is that the E545K samples (given the sample swap) pair off in the tree closely based on the date the sample was prepared, but the EV samples do not (treatment vs non is the primary clustering factor for those).  This suggests a strong batch effect that's specific to E545K cells, and is worse at later dates.  Let's swap the samples in the metadata file, and model the batch effect in a new Sleuth analysis in order to recover the real differential expression signal for the factors of interest: cell line, treatment and cell line/treatment interaction.  This means that our new metadata file (e_meta_swapped.tab) should look like this:

sample                path                           cell_line treatment date
E545K_minus_CP_Dec_1  E545K_minus_CP_Dec_1.kallisto  K         ctl       D1
E545K_minus_CP_Dec_2  EV_minus_CP_Dec_2.kallisto     K         ctl       D2
E545K_minus_CP_Nov_30 E545K_minus_CP_Nov_30.kallisto K         ctl       D0
E545K_plus_CP_Dec_1   E545K_plus_CP_Dec_1.kallisto   K         exp       D1
E545K_plus_CP_Dec_2   E545K_plus_CP_Dec_2.kallisto   K         exp       D2
E545K_plus_CP_Nov_30  E545K_plus_CP_Nov_30.kallisto  K         exp       D0
EV_minus_CP_Dec_1     EV_minus_CP_Dec_1.kallisto     EV        ctl       D1
EV_minus_CP_Dec_2     E545K_minus_CP_Dec_2.kallisto  EV        ctl       D2
EV_minus_CP_Nov_30    EV_minus_CP_Nov_30.kallisto    EV        ctl       D0
EV_plus_CP_Dec_1      EV_plus_CP_Dec_1.kallisto      EV        exp       D1
EV_plus_CP_Dec_2      EV_plus_CP_Dec_2.kallisto      EV        exp       D2
EV_plus_CP_Nov_30     EV_plus_CP_Nov_30.kallisto     EV        exp       D0

As I mentioned earlier, the batch effect seems worse (cluster tree separation is higher) at later dates, so we set Nov 30 as the base level for the date factor by giving it the alphabetically first factor level name, D0. Don't use just 0, 1, 2 unless you want to treat the date as a dosage effect (i.e. the hypotheical effect at Dec 2 should be exactly twice as much as at Dec 1).

I tried modelling just a straight date batch effect by adding date to the model.  This yielded 1000+ significant changes, but they were almost mirrors of the E545K effect (e.g. when E545K cells had a change of X for a transcript, the date factor usually had a change of -X).  This indicates that the E545K changes are artifactual, so we need to model the interaction of the date and cell line (date:cell_line). On with the a new (correct) Sleuth analysis!

library(sleuth)

meta <- read.table("e_meta_swapped.tab", header=TRUE)

meta$path <- as.character(meta$path)

so <- sleuth_prep(meta, ~cell_line*treatment+date:cell_line+date)

so <- sleuth_fit(so)

Note that cell_line*treatment is R syntactic sugar for cell_line+treatment+cell_line:treatment. This gives us a model where gene expression is affected by 3 factors and two factor interactions.*

Now, one by one, we need to exclude each of the factors of the experiment from the linear model to test what transcripts are affected by the factor's exclusion.

so <- sleuth_fit(so, ~treatment+date, "no_cell_line")

so <- sleuth_fit(so, ~cell_line*treatment, "no_date")

so <- sleuth_fit(so, ~cell_line+date:cell_line+date, "no_treatment")

so <- sleuth_fit(so, ~cell_line+treatment+date:cell_line+date, "no_interaction")

so <- sleuth_fit(so, ~cell_line*treatment+date, "no_date_interaction")

The significance test of choice here is the Likelihood Ratio Test (LRT). Run it.

so <- sleuth_lrt(so, 'no_cell_line', 'full')

so <- sleuth_lrt(so, 'no_treatment', 'full')

so <- sleuth_lrt(so, 'no_interaction', 'full')

so <- sleuth_lrt(so, 'no_date', 'full')

so <- sleuth_lrt(so, 'no_date_interaction', 'full')

Grab the results, and filter to those with a qval (multiple-testing corrected p-value) of <.05.

results_table_cell_line_lrt <- sleuth_results(so, 'no_cell_line:full', test_type = 'lrt')

results_table_treatment_lrt <- sleuth_results(so, 'no_treatment:full', test_type = 'lrt')

results_table_interaction_lrt <- sleuth_results(so, 'no_interaction:full', test_type = 'lrt')

results_table_date_lrt <- sleuth_results(so, 'no_date:full', test_type = 'lrt')

results_table_date_interaction_lrt <- sleuth_results(so, 'no_date_interaction:full', test_type = 'lrt')

cell_line.lrt.sig_ids <- results_table_cell_line_lrt$target_id[which(results_table_cell_line_lrt$qval < 0.05)]

treatment.lrt.sig_ids <- results_table_treatment_lrt$target_id[which(results_table_treatment_lrt$qval < 0.05)]

interaction.lrt.sig_ids <- results_table_interaction_lrt$target_id[which(results_table_interaction_lrt$qval < 0.05)]

date.lrt.sig_ids <- results_table_date_lrt$target_id[which(results_table_date_lrt$qval < 0.05)]
date_interaction.lrt.sig_ids <- results_table_date_lrt$target_id[which(results_table_date_lrt$qval < 0.05)]

Now, let's make a (non-redundant) list of all of the target IDs that pass the LRT significance filter. I'm excluding the LRT test results for date and date:cell_line because they're huge, and I don't really care to report the specifics of those batch effects unless they are in genes that are also affected by the original factors of interest (cell line, treatment, and their interaction).

lrt_ids <- unique(c(cell_line.lrt.sig_ids, treatment.lrt.sig_ids, interaction.lrt.sig_ids))

The LRT test doesn't generate fold-change estimates for transcripts since it only considers the factors' presence, not their levels. We need to run a Wald test for each factor level (7 total since the date factor has 2 levels to contrast with base level D0).

so <- sleuth_wt(so, 'dateD1')

so <- sleuth_wt(so, 'dateD2')

so <- sleuth_wt(so, 'cell_lineK')

so <- sleuth_wt(so, 'cell_lineK:treatmentexp')

so <- sleuth_wt(so, 'treatmentexp')

so <- sleuth_wt(so, 'cell_lineK:dateD1')

so <- sleuth_wt(so, 'cell_lineK:dateD2')

Grab the results.

results_table_date2_wt <- sleuth_results(so, 'dateD2')

results_table_date1_wt <- sleuth_results(so, 'dateD1')

results_table_date1_int_wt <- sleuth_results(so, 'cell_lineK:dateD1')

results_table_date2_int_wt <- sleuth_results(so, 'cell_lineK:dateD2')

results_table_cell_line_wt <- sleuth_results(so, 'cell_lineK')

results_table_interaction_wt <- sleuth_results(so, 'cell_lineK:treatmentexp')

results_table_treatment_wt <- sleuth_results(so, 'treatmentexp')

Let's report out only the Wald test results of each factor contrast that correspond to the LRT-passing targets (transcripts) identified earlier in the list "ids".

shared_results_cell_line <- results_table_cell_line_wt[results_table_cell_line_wt$target_id %in% lrt_ids,]

shared_results_treatment <- results_table_treatment_wt[results_table_treatment_wt$target_id %in% lrt_ids,]

shared_results_interaction <- results_table_interaction_wt[results_table_interaction_wt$target_id %in% lrt_ids,]

shared_results_date2_int <- results_table_date2_int_wt[results_table_date2_int_wt$target_id %in% lrt_ids,]

shared_results_date1_int <- results_table_date1_int_wt[results_table_date1_int_wt$target_id %in% lrt_ids,]

shared_results_date1 <- results_table_date1_wt[results_table_date1_wt$target_id %in% lrt_ids,]

shared_results_date2 <- results_table_date2_wt[results_table_date2_wt$target_id %in% lrt_ids,]

Since we want to report out all the factors together (i.e. influence of each factor for a transcript are reported on the same row), we need to sort each factor's Wald test results alphabetically by ID for consistency.

shared_results_cell_line <- shared_results_cell_line[sort.list(shared_results_cell_line[,1]),]

shared_results_treatment <- shared_results_treatment[sort.list(shared_results_treatment[,1]),]

shared_results_interaction <- shared_results_interaction[sort.list(shared_results_interaction[,1]),]

shared_results_date2_int <- shared_results_date2_int[sort.list(shared_results_date2_int[,1]),]

shared_results_date1_int <- shared_results_date1_int[sort.list(shared_results_date1_int[,1]),]

shared_results_date1 <- shared_results_date1[sort.list(shared_results_date1[,1]),]

shared_results_date2 <- shared_results_date2[sort.list(shared_results_date2[,1]),]

Write out all the LRT-significant data to a file. I could get fancy and splice out specific columns to keep, but that's a lot of syntax and it's easier to do that in Excel or with UNIX commands like cut.

write.csv(c(shared_results_cell_line,shared_results_treatment,shared_results_interaction,shared_results_date2_int,shared_results_date1_int,shared_results_date1,shared_results_date2), "e_kallisto_swapped_lrt_q0.05_max.csv")

This looks a lot better! Now we have ~330 genes that show an interaction effect for the cell line and treatment whereas before the swap and cell line specific batch effect modelling we had nothing.  It took some work for sure, but a lot less than redoing the experiment and resequencing everything. This was a quick run to get some directions for grant applications based on pathway analysis etc., so overall patterns are what matter. They'll want to run follow-up experiments and qPCR on salient transcripts in the existing samples to verify the results.

_____________

*If you tried to include the treatment:date interaction and other factor combinations in the model with these just 12 samples, you're likely to get a nasty message like:

Error in solve.default(t(X) %*% X) : 

  Lapack routine dgesv: system is exactly singular: U[7,7] = 0

Which in essence means that you only have one sample for each combination of factors. That'd leave zero degrees of freedom for statistical tests. That's no good, to put it lightly. There may be an interaction effect of the treatment and date, but as per the distance matrix shown earlier, any such effect is negligible compared to the cell line batch effect, so we'll make do without modelling it.

Transcriptome analysis without an annotated genome: pretty fast and slightly dirty

tl;dr If you're doing RNASeq for a species without a well-annotated, contiguous genome, and want useful differential expression results "quickly", run Trinity in genome guided mode, then Kallisto+Sleuth. Annotated functionally using Diamond and a few Perl one liners (below). It'll take less than a week, mostly hands off in Trinity run time.

_________________

An interesting RNA sequencing project came across my desk recently, where a group was looking at a species of fish (with a published genome) they'd corralled and raised at multiple sites upstream and downstream of a riparian urban center. The idea is to look at fish gene signatures for fauna water quality stress.

These days there are a lot of genomes out there that have been generated by groups that basically got funding to run a few HiSeq lanes, and they push out assemblies with thousands (if not more) contigs to get a paper out of it. I don't blame them, as getting a nice genome of contiguous chromosomes is a lot of work, though becoming easier with long reads from PacBio and Oxford Nanopore instruments. The academic rewards these days are not that high for such an endeavour, especially if your bailiwick is molecular biology, ecology, physiology, etc..  Even if you get a nice genome, gene structure and functional annotation in eukaryotes (except maybe fungal genomes, which I was involved in automating) takes a lot of manpower unless you can convince Ensembl your assembly is ripe for them to take on. The fish genome in question is in this dead space: no gene models, lots of contigs.

That being said, as soon as a "genome" sequence is available, researchers often assume you can easily do RNASeq experiments now with good functional analysis and the whole nine yards. I've done this before with Tophat + Cufflinks + Cuffdiff in de novo gene modelling mode for other organisms, which is okay but has issues with genes split across contigs, missassemblies, isoform detection, chimeric genes, etc. Since those tools were in vogue there have been a few advances in RNASeq processing, so I thought it's worth seeing what alternatives can do.

First off, don't map back to NCBI's UniGene. I tried this since it's easy and got nothing, which is due to the several factors such as poor coverage of genes expressed in the tissue of interest (liver), and the fact that UniGene gives you the best representative sequence for each gene, not a consensus transcript. The representative is the longest good quality read in dbEST for a gene, and because these mostly come from (poly A-tail amplified) ESTs, they tend to have a 3' mRNA bias which leads to a whole host of issues for differential expression analysis especially in fish, which often have 8 copies of genes (they went through an extra whole genome duplication, so ohnologs that have 4 copies in most vertebrates have 8 in teleosts, i.e. most fish).

Given the volume of tissue specific transcript information you get these days for cheap, why not do a de novo transcript assembly, then run a quick differential expression workflow like I blogged in Using Kallisto & Sleuth for RNASeq analysis?  In my case, this yielded about 50 genes in logical functional categories that were candidates for further study in water quality stress response. Here are the details of how to proceed.

I downloaded the reference genome from the NCBI for the fish in question.

wget  ftp://ftp.ncbi.nlm.nih.gov/sra/wgs_aux/JN/CD/JNCD01/JNCD01.1.fsa_nt.gz 

Did the same for FastA parts 2 and 3. Concatenated the files to get a final FastA genome.

cat JNCD01.*.fsa_nt.gz | gzip -cd - > genome.fna

BWA (version 0.7+) indexed the genome. If you need samtools, bwa, Trinity, etc., you may find it easiest to just install BioBuilds. Then...

bwa index genome.fna

Ran BWA against the genome for all the samples, in this case X = {site0, site1, site2, site3, site4}, and Y = {1,2,3,4} since there were four fish sampled from each site.

bwa mem -t 10 genome.fna <(gzip -cd siteX_sampleY.R1.fastq.gz) <(gzip -cd siteX_sampleY.R2.fastq.gz) | samtools sort -@ 4 -T siteX_sampleY -O bam > siteX_sampleY.bam

Merged all the BAMs into one big one (using 16 threads):

samtools merge -r -@ 16 -O BAM all_samples.bam *.bam

Ran Trinity on reference-genome-guided mode with these data to generate consensus transcripts.  NOTE: unlike most genome-guided methods, Trinity only uses the genome to pre-cluster the reads for the de novo assembly, so it doesn't inherit all the gaps and flaws of the "genome". 

I have a lot of threads available on my machine, so I let it use 80.  YMMV.

Trinity --genome_guided_bam all_samples.bam   --genome_guided_max_intron 10000 --max_memory 450G --CPU 80 

If you didn't have a reference genome, you would have just skipped the bwa bit, the samtools bit and previous line, going straight to this instead:

Trinity --seqType fq --max_memory 800G --single (<gzip -cd myreads.fastq.gz) --CPU 80 --trimmomatic

At this point, book a last minute, long weekend getaway. Trinity finished after 5 days, generating 510738 transcript contigs in trinity_output_dir. Using a normal Trinity run instead on this volume of data would normally take 2 weeks or more. #winning

I ran these against the NCBI non-redundant protein database with Diamond (blastx mode), 

Download the non-redundant protein dataset from the NCBI.

wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz

Uncompress it, then index it with Diamond.

diamond index -i nr nr

Run Diamond on the Trinity results (FastA).  It's hard to overstate how awesomely fast Diamond is at finding good protein homology for DNA sequences.

diamond blastx -d nr -q Trinity-GG.fasta -p 80 -a fathead.nr

Then run MEGAN6's taxonomic classification algorithm (after being sure to install thre "tools" during the setup script, and downloading the latest NCBI acc2tax file from the same site):

/export/common/programs/megan6/tools/blast2lca -i fathead.nr.daa -f DAA -m BlastX -o fathead.nr.diamond.tax -a2t /export/common/dbs/megan6/prot_acc2tax-Oct2017X1.abin

This yielded 142539 contigs with protein matches (the remaining 370K are mostly rRNA, short nonsense fragments/assembly chimera, low quality, etc.).  Of these 136764 match teleost proteins (the rest are Burkholderia, other contaminants, or unclassified). In my case, the fish was evolutionarily fairly close to zebrafish, so most of the top scores had e-values<10e-35. But, you might want to pimp up this one liner with a e-value or % ID filter if you've got an organism with no reasonable model organism in the evolutionary 'hood.

perl -F\; -ane 'print "$F[0]\n" if /Actinopteri/' fathead.nr.diamond.tax > fathead.nr.diamond.tax.fish_matches.txt

perl -ne 'BEGIN{%keep=split /(\s)/s, `cat fathead.nr.diamond.tax.fish_matches.txt`; $/=">"}print if /^>?(\S+)/ and $keep{$1}' Trinity-GG.fasta > Trinity-GG.fish_keep.fna

Index the kept contigs with kallisto:

kallisto index -i fishname Trinity-GG.fish_keep.fna

Run kallisto against the fishy (that's a good thing) contigs for each sample.

kallisto quant -i fishname.kallisto_index -b 100 -l 180 -s 20 -o siteX_sampleY.kallisto <(gzip -cd siteX_sample1.R1.fastq.gz) <(gzip -cd siteX_sampleY.R2.fastq.gz)

Kallisto finished running the reads against the Trinity assembled transcripts.  

I then mapped the Trinity transcript assembly IDs to the best Diamond matches in the NCBI (RefSeq NP_ preferred, then XP_, then whatever):

perl -ane 'if(/(NP_.*?)\s/ and not $p{$F[0]}++){$f{$F[0]} = "$1\t$F[10]"}elsif(/(XP_.*?)\s/ and not $p{$F[0]} and not $x{$F[0]}++){$f{$F[0]} = "$1\t$F[10]"}elsif(not $p{$F[0]} and not $x{$F[0]}){$f{$F[0]} = "$F[1]\t$F[10]"}END{for(keys %f){print "$_\t$f{$_}\n"}}' fathead.nr.diamond.tab > trinity2refid

Because it's quite likely that we only have partial coverage of genes, I want to consolidate our assembled transcripts into gene symbol level as much as possible before differential expression analysis. To this end, I'll map the refids to gene names in IPA (you could do this with DAVID too for example if you don't have an IPA license) and save to ref2symbol_ipa.txt. I then mapped the IDs using the Entrez Gene database.

wget ftp://ftp.ncbi.nih.gov/gene/DATA/gene2refseq.gz

perl -F\\t mv -ane 'print "$F[5]\t$F[15]"' gene2refseq > refid2symbol_entrez.txt

I also mapped the RefSeq IDs to UniProtKB here, then extracted the UniProtKB IDs to gene names (saved together as refid2symbol_uniprot.txt).

perl -F\\t -ane 'print "$F[0]\t$1\n" if $F[6] =~ /^(\S+)/' uniprot-yourlist.tab > refid2symbol_uniprot.txt

In the end I had a file of Trinity IDs and their corresponding gene (if available from IPA, Entrez Gene or UniProt) or RefSeq transcript by running this:

perl -ane 'BEGIN{%r2i=split /\s/s,`cat refid2symbol_ipa.txt`; %r2e=split /\s/s,`cat refid2symbol_ncbi.txt`; %r2u=split /\s/s,`cat refid2symbol_uniprot.txt`} print "$F[0]\t", ($r2i{$F[1]} || $r2e{$F[1]} || $r2u{$F[1]} || $F[1]), "\n"' trinity2refid > trinity2symbol.txt

I ran the following Sleuth (R) analysis which generates a table of the expression values for all sites where there is a significant log-ratio test qval (multiple testing corrected p-value) for the factor.  Most of the code from here down is a variation on my earlier post Using Kallisto & Sleuth for RNASeq analysis.

library(sleuth)

meta <- read.table("metadata.tab", header=TRUE)
trinity2symbol <- read.table("trinity2symbol.txt")

meta$path <- as.character(meta$path)

so <- sleuth_prep(meta, ~site, target_mapping=trinity2symbol, use_extra_bootstraps=TRUE)

so <- sleuth_fit(so)

so <- sleuth_fit(so, ~1, "reduced")

so <- sleuth_lrt(so, 'reduced', 'full')

results_table_lrt <- sleuth_results(so, 'reduced:full', test_type = 'lrt’)

results_table_lrt <- results_table_lrt[sort.list(results_table_lrt[,1]),]

lrt.sig_ids <- results_table_lrt$target_id[which(results_table_lrt$qval < 0.05)]

Now we want to perform the Wald test for each site vs. site0 (our factor base level) since this will give us the fold-changes and site-specific significance of the differential expression.

so <- sleuth_wt(so, 'site1')

so <- sleuth_wt(so, 'site2')

so <- sleuth_wt(so, 'site3')

so <- sleuth_wt(so, 'site4')

results_table_wt_1 <- sleuth_results(so, 'site1')

results_table_wt_2 <- sleuth_results(so, 'site2')

results_table_wt_3 <- sleuth_results(so, 'site3')

results_table_wt_4 <- sleuth_results(so, 'site4')

We'll want to report the beta (~fold change) and wald test q-value for each sample for each ID that passed the original log-ratio-test.

shared_results_1 <- results_table_wt_1[results_table_wt_1$target_id %in% lrt.sig_ids,]

shared_results_2 <- results_table_wt_2[results_table_wt_2$target_id %in% lrt.sig_ids,]

shared_results_3 <- results_table_wt_3[results_table_wt_3$target_id %in% lrt.sig_ids,]

shared_results_4 <- results_table_wt_5[results_table_wt_4$target_id %in% lrt.sig_ids,]

It's important at this point to make sure we have all the rows in each sample in the same index position for collating the output table, and the easiest way to do this is to sort all the results by target_id (transcript name) alphabetically, like I snuck in earlier for the LRT test results in case you weren't paying attention.

shared_results_1 <- shared_results_1[sort.list(shared_results_1[,1]),]

shared_results_2 <- shared_results_2[sort.list(shared_results_2[,1]),]

shared_results_3 <- shared_results_3[sort.list(shared_results_3[,1]),]

shared_results_4 <- shared_results_4[sort.list(shared_results_4[,1]),]

Collate the results and write to a file. The first column is the significant (q-value < .05) log-ratio test score information for the factor for each gene across all the sites, while the subsequent columns are the Wald scores and log fold-changes (kinda) for each specific site.  Note that because we have multiple values for the factor, a gene can pass the LRT test, but have insufficient replicates at any given site to have a good (qval <.05) Wald score. This is the opposite of the binary factors people usually are testing with RNASeq (e.g., treatment vs. control). Sequencing more replicates is left as an exercise to the reader.

combined_result <- data.frame(results_table_lrt[results_table_lrt$target_id %in% lrt.sig_ids,], site1_qval=shared_results_1[,"qval"], site1_b=shared_results_1

[,"b"], site2_qval=shared_results_2[,"qval"], site2_b=shared_results_2[,"b"], site3_qval=shared_results_3[,"qval"], site3_b=shared_results_3[,"b"], 

site4_qval=shared_results_4[,"qval"], site4_b=shared_results_4[,"b"])

write.table(combined_result, "all_sites_lrt_passed_qval.05.txt", sep="\t")

Now that we have the diff exp table, let's annotate the transcripts that had protein hits (note the literal ctrl-A in the command):

perl -ne 'BEGIN{%t=reverse(split /\s/s, `cut -f 2,1 trinity2refid`); $/=">"}print "$1\t$2\n" if /^(\S+)\s+(.*?)^A/ and exists $t{$1}' nr > trinity_descs.txt

sort tids | perl -ane 'BEGIN{%d=split /[\t\n]/s, `cut -f 1,2 id_desc.txt`}print $d{$F[0]},"\n"' | perl -ane 'BEGIN{%d=split /[\t\n]/s, `cut -f 2,3 trinity_descs.txt`}print $d{$F[0]},"\n"' > descs

Lets also append the p-value of the protein hit used for the annotation:

perl -F\\t -lane 'BEGIN{%e=split /\s/s, `cut -f 1,3 trinity2refid`}print $_,"\t",$e{$F[0]},"\n"' all_sites_lrt_passed_qval_desc.05.txt > lrt_any_factor_qval_lt_0.05_desc_diamond_evalues.txt

That's it!  Now you can go down the rabid whole of functional analysis, a blog post I'll get to writing eventually...