CHGI Bioinformatics

March 1, 2017

This is not normal: dealing with outliers, correlated factors, and skewed factor effects in RNASeq (a Kallisto/Sleuth/DESeq2/U-Test/skewness mash-up)

tl;dr there are RNA experiments where despite reasonable sample selection, no significant changes will be found. Especially in patient-cenetered studies, sometimes you just can't get true "replicates".

In three stages we can correct for outliers, spurious correlations, and a quantitative factor like patient age that has a skewed (as opposed to Gaussian) dosage effect on transcripts. This process gets us from 2 significantly regulated transcripts to 285 (q-value < 0.05). This involves using Kallisto, Sleuth, DESeq2, and some home brewed code around R's built-in Wilcox-Mann-Whitney U Test and skewness measure.
__________________________

Let's start with an experiment designed with a reasonable balance of males and females in each of three age categories, under 40, over 60, and in between. I noticed on the sample submission form that the input RNA concentrations were unusually divergent, so I've included those values in case they were informative (might affect PCR duplicate levels).

library(sleuth)

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

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

meta$age_group <- as.factor(meta$age_group)

Stage 1: Identifying outlier samples

First, we need to collate the normalized TPM transcript abundance data as per my previous blog post. You could probably do it with the so info you've already got from the above analysis, but it's trickier I find. Let's load all the transcripts with expression values using DESeq2.

abundance_full <- read.table("all_abundance.tsv", header=TRUE, row.names=1)

abundance_full_integer <- apply(abundance_full, c(1,2), function(x){as.integer(x)})

library("DESeq2")

dds <- DESeqDataSetFromMatrix(countData=abundance_full_integer, colData=data.frame(meta), design=~sex*age_group+age+sex:age+conc)

dds_redux <- dds[rowSums(counts(dds)) > 1,]

Normalize the count data so that the variance is similar across the range of expression values (necessary for a good Poisson distance measure). This is a nice, easy to use feature of DESeq2. Then visualize the distance matrix.

rld <- rlog(dds_redux, blind=FALSE)

library("pheatmap")

library("RColorBrewer")

library("PoiClaClu")

poisd <- PoissonDistance(t(abundance_full))

samplePoisDistMatrix <- as.matrix( poisd$dd )
rownames(samplePoisDistMatrix) <- paste(meta$sex, meta$age, sep=":")

colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)

pheatmap(samplePoisDistMatrix, clustering_distance_rows=poisd$dd,

clustering_distance_cols=poisd$dd, col=colors)

Let's also check that in closely related samples (such as 8 & 12 in the lower right corner) show fairly even variance across the range of expression values, otherwise our nice matrix isn't very meaningful.

plot(assay(rld)[,c(8,12)], xlim=c(0,16), ylim=c(0,16), pch=16, cex=0.3)

Looks good, i.e. it's fairly evenly thick above zero, and there's no feathering in the lower left. Using just a log 2 transform on the Kallisto TPM data yields a different tree, and the scatter plot above is a lot fatter at the bottom, so it was worth using DESeq2's normalization.

Back to the outliers. In the matrix, you'll see that there are three samples, with indices 4, 5, and 14 along the bottom that are strong outliers. Not only are they strong outliers, but they are closely related to each other despite not sharing either a sex or age factor level. This gives us a strong suspicion that there's something wrong here like cross-contamination, prep batch effect, etc. but it's not obvious what. Let's exclude these samples and see how much the model of known factors improves.

> meta_no_outliers <- meta[c(1:3,6:13,15:18),]

> so_no_outliers <- sleuth_prep(meta_no_outliers, ~sex*age_group+age+sex:age+conc)

...

> so_no_outliers <- sleuth_fit(so_no_outliers)

...

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

[1] 193932.1

The 194K residual sum of squares (RSS) is a lot better than the 347K we had in the original analysis with all 18 samples! If we work through the original analysis with just the 15, we get...

> length(lrt_age.sig_ids)

[1] 8

> length(lrt_sex.sig_ids)

[1] 8

> length(lrt_int.sig_ids)

[1] 4

> length(lrt_conc.sig_ids)

[1] 71

D'oh! Looks like the concentration of input RNA is explaining most of the variance. Collectively this is 84 genes, which we will write out to a file, even if the concentration effect is suspect.

Stage 2: Spurious factor correlations

Getting back to our potential sources of failure in the analysis, let's check for spurious correlation between the concentration and the age factor (the main factor of interest). Plot the data points and a trend line.

> plot(data.frame(meta$age,meta$conc))

> abline(lm(meta$age ~ meta$conc))

Ouch, the concentration as a dosage factor is probably absorbing age dosage effect. The same plot for the outlier-less metadata set eliminates the lower right dot, leaving us with two really low concentrations that are centered in the age distribution (50 & 52). That's lucky for us, because it means that we can change the concentration numeric factor into a binary one to capture just those two points without disturbing the balance of the age factor (e.g. if we'd set the low threshold to 100 it'd have all sample that are older and could still confound the age factor modelling). It's the best we can do with the card we've been dealt.

> plot(data.frame(meta$age,meta$conc))

> meta$conc_low <- meta$conc <= 55

> meta$conc_low

[1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

[13] FALSE TRUE FALSE FALSE TRUE FALSE

In case you're wondering how just excluding the outlier samples and the concentration from the model works out, you get just one significant change, only in age_group1 (between 40yo and 60yo), and it's a transcript at that showed up in the previous analyses repeatedly. You'll see shortly how the binary concentration modelling works out a lot better.

Stage 3: Rationally defining a factor level threshold for a skewed factor effect (numeric factor -> two-level factor)

Let's keep only the samples with a TPM of at least 1 on average.

> abundance <- abundance_full[rowSums(abundance_full)>length(meta_clean$age),]

Let's define all the combinations of young/old possible by setting the young/old threshold to each observed value in turn.

> oldness_treshold_combos <- t(sapply(meta_clean$age, function(x){meta_clean$age > x}))

Let's define a function to run the Wilcox-Mann-Whitney U Test, as I have in a previous post. As an aside, I also tried using the Kruskal-Wallis H test which extends the U Test notion to three or more factor levels (using the 0/1/2 age_group factor from earlier), but it's underpowered with this number of samples. As you'll see below, we'll have enough trouble with Wilcox already. We'll stick the distribution of old/young across samples in a global variable called "m" later.

> wil <- function(x){d <- data.frame(tpm=t(abundance[x,])[,1], old=m); w <- wilcox.test(tpm ~ old, data=d); w$p.value}

Now define the function that assigns a young/old combination to "m", then runs the Wilcox test and finally applies the FDR (Benjamini-Hochberg) p-value multiple testing correction to give us "q-values". The initial check for the sum combo below is a special condition for the oldest sample, which gives no age contrast if used as a threshold and buggers up the whole process.

> perm_qvals <- function(combo){if(!sum(combo)){return(rep_len(c(NA), length(abundance[,1])))}; m <<- as.factor(combo); pvals <<- c(); for(i in 1:length(abundance[,1])){pvals <<- append(pvals,wil(i))}; return(p.adjust(pvals, method="fdr"))}

Run the qval test on each combination we calculated earlier.

> combo_results <- apply(oldness_treshold_combos, 1, perm_qvals)

There were 50 or more warnings (use warnings() to see the first 50) [you can ignore these, it's complaining about ties in the rankings for the Wilcox test for certain transcripts, in which case it can't calculate an exact p-value)]

> apply(combo_results, 2, function(x){sum(as.integer(x<.05))})

[1] 285 NA 0 0 0 0 0 0 0 0 0 0 0 0 0

> res <- combo_results[,1]

> names(res) <- rownames(abundance)

It may seem strange that only the first threshold (which is old = age > 50 if you look back at the order of samples in meta.tab) has significant changes, but you have to keep in mind, as we determined in a previous post, that you need at least eight samples in each factor level to make the U Test p-value survive FDR multiple testing correction in Kallisto-based RNASeq where you'll have 18,000+ tests (i.e. expressed mRNA isoforms in a tissue). Only an 8/7 split comes close to satisfying that when we have 15 samples, like we do here.

Let's try the analysis again, but see what happens if we don't exclude the outliers, so we have 18 samples, which will give us multiple 8+ split options (8/10, 9/9,10/8):

> abundance <- abundance_full[rowSums(abundance_full)>length(meta$age),]

> rownames(abundance) <- rownames(abundance_full[rowSums(abundance_full)>length(meta$age),])

> oldness_treshold_combos <- t(sapply(meta$age, function(x){meta$age > x}))

> combo_results_full <- apply(oldness_treshold_combos, 1, perm_qvals)

apply(combo_results_full, 2, function(x){sum(as.integer(x<.05))})

[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 NA 0 0 0 0

It would seem in this case that the benefit of added samples for the U Test is outweighed by the outlying nature of the expression values in those 3 extra samples, as we get no significant changes. The three samples are not just overall outliers, but also specifically in their response to the experimental factor of interest (age). That is to say, they disrupt the rank of genes grouped as young/old, not just the replicate normality assumption of the Wald and LRT tests on the TPM values. Had we gotten good results in the combo_results_full vector, we might have included the outliers and proceeded with caution in comparing the betas from Sleuth with either 15 or 18 samples.

Let's note the significant IDs from the 15 sample test for later when we want to print out results to a file. I'm doing this in two steps, for clarity.

> sig_indices <- res < .05

> ids <- sort(rownames(abundance)[sig_indices])

Let's run with the >50 cutoff in the 15 samples and redo the Sleuth analysis with both the binary age factor and the quantitative age factor. The beta values are more useful than the raw fold-changes I could calculate directly from the abundance matrix variable because 1) the effect of other factors are subtracted from the data, and 2) especially for low abundance transcripts beta more accurately reflects the effect of the age factor since it properly accounts for measurement error though the bootstrap analysis Sleuth does.

> meta_clean$old <- meta_clean$age > 50

> meta_clean

sample path sex age_group age conc conc_low old

1 57_50Y_P2 57_50Y_P2.kallisto M 1 50 55 TRUE FALSE

2 63_75Y_P2 63_75Y_P2.kallisto M 2 75 99 FALSE TRUE

3 SATP100 SATP100.kallisto M 0 27 203 FALSE FALSE

6 SATP119 SATP119.kallisto M 2 74 511 FALSE TRUE

7 SATP120 SATP120.kallisto M 1 41 466 FALSE FALSE

8 SATP125 SATP125.kallisto F 0 30 697 FALSE FALSE

9 SATP137 SATP137.kallisto M 0 19 183 FALSE FALSE

10 SATP140 SATP140.kallisto F 0 39 264 FALSE FALSE

11 SATP154 SATP154.kallisto F 2 61 256 FALSE TRUE

12 SATP68 SATP68.kallisto M 0 37 303 FALSE FALSE

13 SATP75 SATP75.kallisto M 2 65 112 FALSE TRUE

15 SATP91 SATP91.kallisto F 1 57 736 FALSE TRUE

16 SATP92 SATP92.kallisto F 1 53 136 FALSE TRUE

17 SATP95 SATP95.kallisto F 1 52 53 TRUE TRUE

18 SATP98 SATP98.kallisto F 2 73 83 FALSE TRUE

Luckily for us, the age threshold of >50 puts one low concentration sample in the young category and one in the old, so both factors can be included in the final model.

> so_final <- sleuth_prep(meta_clean, ~sex*old+age:old+conc_low)

> so_final <- sleuth_fit(so_final)

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

[1] 266102.6

By the way, did our spurious correlation removal from Stage 2 help? Let's see what the RSS is without it in the model.

> so_final <- sleuth_fit(so_final, ~sex*old+age:old, "no_conc")

> sum(so_final$fits$no_conc$summary[,2])

[1] 300209.5

Nice, it reduces the noise by over 10%. It's not as good an RSS as we had with the quantitative concentration factor, but it's much less biased. Let's go on with the inclusive model and calculate the betas/fold-changes with a Wald test.

> so_final <- sleuth_wt(so_final, "sexM")

> so_final <- sleuth_wt(so_final, "oldTRUE")

> so_final <- sleuth_wt(so_final, "sexM:oldTRUE")

> so_final <- sleuth_wt(so_final, "oldTRUE:age")

> so_final <- sleuth_wt(so_final, "oldFALSE:age")

> so_final <- sleuth_wt(so_final, "conc_lowTRUE")

> wt_sex <- sleuth_results(so_final, 'sexM')

> wt_old <- sleuth_results(so_final, 'oldTRUE')

> wt_int <- sleuth_results(so_final, 'sexM:oldTRUE')

> wt_old_age <- sleuth_results(so_final, 'oldTRUE:age')

> wt_young_age <- sleuth_results(so_final, 'oldFALSE:age')

> wt_conc_low <- sleuth_results(so_final, 'conc_lowTRUE')

> table(wt_sex[,"qval"] <.05)

FALSE TRUE

28142 2

> table(wt_old[,"qval"] <.05)

FALSE TRUE

28142 2

> table(wt_int[,"qval"] <.05)

FALSE TRUE

28142 2

> table(wt_conc[,"qval"] <.05)

FALSE TRUE

28133 11

> table(wt_old_age[,"qval"] <.05)

FALSE TRUE

28142 2

> table(wt_young_age[,"qval"] <.05)

FALSE TRUE

28143 1

So we have very few significant q values, but that's to be expected, because the Wald test assumes you have a normal distribution and an accurate assessment of variance. We have neither, but we want the beta values that the test generates. These will still be our most accurate assessment of factor driven fold-change in transcript abundance.

wt_shared_sex <- wt_sex[wt_sex$target_id %in% ids,]

wt_shared_old <- wt_sex[wt_old$target_id %in% ids,]

wt_shared_int <- wt_int[wt_int$target_id %in% ids,]

wt_shared_old_age <- wt_old_age[wt_old_age$target_id %in% ids,]

wt_shared_young_age <- wt_old_age[wt_young_age$target_id %in% ids,]

wt_shared_young_age <- wt_old_age[wt_conc$target_id %in% ids,]

wt_shared_conc_low <- wt_conc_low[wt_conc_low$target_id %in% ids,]

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

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

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

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

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

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

sig_abundance <- sig_abundance[sort.list(rownames(sig_abundance)),]

If the sample values going into the fold-changes for the age factor level (old/young) are not normally distributed, is there something useful we can say about their distribution? One summary statistic we can grasp fairly easily is skewness, which when negative indicates that the mean is less than the median, and conversely that when positive the mean is greater than the median. From the Wikipedia entry:

In our data for example, a negative skew would mean that the effect of the fold-change is more intense in the older people within the old group. A positive skew would mean that the FC values for the old-designated samples pile up closer to the old cutoff (>50yo) so increased age isn't a big a factor. Let's see if the skew of the abundance of the 285 significant transcripts are quantitatively different from the general transcript abundance.

> library(moments)

> sig_abundance <- abundance[sig_indices,]

> sk <- apply(sig_abundance, 1, function(x){skewness(as.vector(x[meta_clean$old]))})

> summary(sk)

Min. 1st Qu. Median Mean 3rd Qu. Max.

-1.72200 -0.29110 0.09093 0.07834 0.44870 1.50800

> sk0 <- apply(abundance, 1, function(x){skewness(as.vector(x[meta_clean$old]))})

> summary(sk0)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA's

-2.0620 -0.1577 0.2699 0.3049 0.7294 2.2680 62

Looks like the 285 transcripts have a long tail more to the left (more negative) than the whole set of transcripts. Let's check more thoroughly. First, are both distributions (of skewness) normal? Run the Shapiro-Wilk test to find out.

> shapiro.test(sk)

...

W = 0.9934, p-value = 0.2449

> shapiro.test(sample(sk0, length(ids)))

...

W = 0.99399, p-value = 0.3266

Yes. Do they have the same variance? Run an F-test, or you could try Levene's test or Bartlett's test, etc.. they all fail for this dataset.

> var.test(sk, sk0)

F test to compare two variances

data: sk and sk0

F = 0.699, num df = 284, denom df = 19418, p-value = 7.016e-05

alternative hypothesis: true ratio of variances is not equal to 1

95 percent confidence interval:

0.5961302 0.8307902

sample estimates:

ratio of variances

0.6990001

Which makes running a t-test moot since it assumes equality of variance in each sample group. What we can say easily is that the skewness is further left and tighter (about 1/3 as much) in the significantly changes transcripts than in a random sample of 285 transcripts. That means many transcripts are affected in the really old a lot more than the moderately old, and luckily we don't have to give a specific threshold (like 60 in the original analysis) to catch that phenomenon, nor did the effect have to be normally distributed in our samples.

> sum(sk)

[1] 22.32632

> sum(sample(sk0,length(sk)), na.rm=TRUE)

[1] 78.9279

> mean(sk)

[1] 0.07833795

> mean(sk0, na.rm=TRUE)

[1] 0.3049191

For all the 285 U-test passing transcripts let's write out all the factor contrasts we've calculated so far, and the old samples skew stat adjusted for sk0 (to help account for the effect of the specific distribution of sample ages we have in the analysis).

sk_adjusted <- sk - (mean(sk0, na.rm=TRUE)-mean(sk))

> d <- data.frame(ids,avg_tpm=rowMeans(sig_abundance[ids,]),u_test_qval=res[ids],skew=sk_adjusted,wt_shared_sex[,"b"],wt_shared_sex[,"qval"],wt_shared_old[,"b"],wt_shared_old[,"qval"],wt_shared_int[,"qval"],wt_shared_int[,"b"],wt_shared_young_age[,"qval"],wt_shared_young_age[,"b"],wt_shared_old_age[,"qval"],wt_shared_old_age[,"b"])

> write.csv(d, "u_test_no_outliers_q_lt_0.05.csv")

Pondering the merits and validity of including a kurtosis measure is left as an exercise for the reader.

So, that's one way to analyze RNASeq when the expected expression values aren't necessarily normally distributed and don't have a reasonable estimate of variance. Now you have two gene lists to play with: those that appear normally distributed amongst the age predefined groups (84, mostly attributable to input RNA concentration conflated with age), and those that show a less/greater then 50yo dichtomy (285). Some transcripts on the former list are also in the latter list. While the beta (~fold change) values in the latter aren't super meaningful, taken together with the skewness value for the transcript you can focus on transcript changes that are skewed to the really old, or more middle aged, etc.

P.S.: significant changes were corroborated with consistency scores in pathway analysis (IPA) from both the 84 and 285 item lists, albeit different pathways, therefore it's probably worth following up on both.

February 16, 2017

RNASeq mic drop: U Test for detecting differential expression when factor dosage is unknown

tl;dr If you have an unknown dosage factor, make it binary, and run a Wilcox-Mann-Whitney U test to improve statistical power if you have ~20 samples or more. In a real test case of 20 samples, we went from 2 significant changes using the likelihood ratio test (LRT) for each sample group independently, to over 1000 significant using the (non-parametric) U test.

Suppose you have RNASeq samples that can be divided into two groups with 1) variable dose of a factor, and 2) no dose. Also suppose you don't know the dose for any given sample in group 1. I ran into this issue analyzing environmental sampling across 5 riparian sites (4 samples at each site) up and downstream of a wastewater source. By running a U test between upstream and downstream samples the downstream pollution dosage doesn't matter. That's because it's a non-parametric test that does not make assumptions of normality for the "replicates" in the upstream and downstream groups. Once you have about 20 samples, real differentially expressed genes start to survive Benjamini-Hochberg multiple testing correction.

Names have been changed to protect the unpublished.
____________

We'll pick up my blogged de novo assembly RNASeq analysis after the Kallisto mapping step, so we have a directory called samplename_trinity.kallisto for each sample. We make an experiment metadata file (metadata_trinity.tab) like so where we note whether the site is upstream or downstream:

sample path site downstream
15 15_trinity.kallisto site_4 1
17 17_trinity.kallisto site_4 1
21 21_trinity.kallisto site_4 1
23 23_trinity.kallisto site_4 1
29 29_trinity.kallisto site_2 0
30 30_trinity.kallisto site_2 0
34 34_trinity.kallisto site_2 0
35 35_trinity.kallisto site_2 0
40 40_trinity.kallisto site_1 0
43 43_trinity.kallisto site_1 0
4 4_trinity.kallisto site_5 1
50 50_trinity.kallisto site_1 0
57 57_trinity.kallisto site_1 0
59 59_trinity.kallisto site_3 1
62 62_trinity.kallisto site_3 1
68 68_trinity.kallisto site_3 1
6 6_trinity.kallisto site_5 1
70 70_trinity.kallisto site_3 1
7 7_trinity.kallisto site_5 1
8 8_trinity.kallisto site_5 1

If you used a reference genome instead of Trinity, just follow omit the "_trinity" from the path values.

First, let's collate a Transcripts Per Million (TPM) file for all transcripts as produced by Kallisto for each sample. This is a normalized measure of transcript expression. Start with generating the header line:

perl -e 'print "target_id\t",join("\t",map {/(.*)\//;$1} @ARGV),"\n";' *_trinity.kallisto/abundance.tsv > all_trinity_abundance.tsv

Then append the TPM values:

paste *_trinity.kallisto/abundance.tsv | perl -ane 'print $F[0];for (1..$#F){print "\t$F[$_]" if /[49]$/}print "\n"' | tail -n +2 >> all_trinity_abundance.tsv

Alternatively if you want to use FPKM (i.e. normalized for gene length too):

paste *.kallisto/abundance.tsv | tail -n +2 | perl -ane 'print $F[0];for (1..$#F){print "\t",$F[$_]/$F[$_-2]*1000 if /[49]$/}print "\n"' >> all_abundance_fpkm.tsv

Now we can start R, load the data and get to some statistics!

meta <- read.table("metadata_trinity.tab", header=TRUE, row.names=TRUE)
abundance_full <- read.table("all_trinity_abundance.tsv", header=TRUE)

Let's only keep the transcripts that have an average of 1 TPM or more, by filtering rows with a sum of less than 20. Most technical noise bootstrap estimates in Kallisto are less than 1 in the previous analysis of this dataset.

abundance <- abundance_full[rowSums(abundance_full)>20,]
rownames(abundance) <- rownames(abundance_full[rowSums(abundance_full)>20,])

Now, for each transcript let's calculate the log₂ fold change between the downstream sites and the base level of the factor, upstream.

log2fc <- function(x){log2(sum(abundance[x,meta$downstream])/sum(abundance[x,!meta$downstream]))}
l <- c()
for(i in 1:length(abundance[,1])){l <- append(l,log2fc(i))}

Now let's run our Wilcox-Mann-Whitney U test for each transcript. The test sorts the TPMs numerically and checks if there's a non-random grouping of downstream TPM values in the rankings. If all downstream sites have the lowest or highest TPMs for a transcript, that yields the smallest p-value.

wil <- function(x){df <- data.frame(tpm=t(abundance[x,])[,1], downstream=meta$downstream); w <- wilcox.test(tpm ~ downstream, data=df); w$p.value}
pvals <- c()
for(i in 1:length(abundance[,1])){pvals <- append(pvals,wil(i))}

Of course, if you have enough transcripts measured, eventually all downstream samples will rank together in a transcript by chance. Let's apply a multiple testing correction to the p-values. The False Discovery Rate (FDR) a.k.a. Benjamini-Hochberg, is quite popular as for most experiments it's neither too conservative nor too lax.

q <- p.adjust(pvals, method="fdr")

Let's decorate the abundance information with the fold change and U test stats, then note those with adjusted p-value (a.k.a. q-value) < 0.05.

abundance[,"log2fc"] <- l
abundance[,"qval"] <- q
qsig <- which(q < 0.05)
length(qsig)
[1] 1073

Hey, that's a lot more than we got using the likelihood ratio test for each individual site, where only two transcripts significant across all the downstream sites! Strength in replicate sample numbers. The q-values won't blow you away with a 8:12 split, they're in the (0.01,0.05) range, which is why I suggest this technique when you have 20 or more samples. You can try with less, but you're most likely just spinning your wheels.

The log2 fold change isn't super useful, as we have no assumption of normality amongst replicates (since the downstream sites have different and unknown pollution exposure levels). BUT we can use it to see the set union of transcripts with a fairly strong effect at one site or consistent expression changes across sites. This could be useful for pollution biomarker discovery.

length(which(abs(abundance[qsig,"log2fc"]) > 1))
[1] 118

That's still quite a bit to work with! Let's write out all the qval passed transcripts to a file:

write.csv(abundance[qsig,], file="up_downstream_wilcox_q_lt_0.05.csv")

Not all of these are going to be real, but a subsequent analysis in IPA shows a statistically significant number of observations in the RNASeq that are consistent with the literature, so we're probably in the right track. To be extra sure, let's randomly label the samples as upstream or downstream and redo the analysis. Hopefully we don't get any significant q-values.

meta$downstream <- sapply(1:length(meta$downstream), function(x) sample(c(0,1),1))
abundance <- abundance_full[rowSums(abundance_full[2:21])>20,2:21]

pvals <- c()
for(i in 1:length(abundance[,1])){pvals <- append(pvals,wil(i))}
q_random <- p.adjust(pvals, method="fdr")
qsig_random <- which(q_random < 0.05)
length(qsig_random)
[1] 0

*mic drop*

Hold on, that's a bit premature. Let's pick up the mic again because we may not be doing a fair random comparison. In the real configuration, the downstream status is consistent with site status for each samples. How many random ways could we do this with 2 upstream and 3 downstream site designations?

library(combinat)

downstream_options <- unique(permn(c(0,0,1,1,1)))

Which gives us the ten unique permutations of the sites assigned binary downstream values.

downstream_options

[[1]]

[1] 0 0 1 1 1

[[2]]

[1] 0 1 0 1 1

[[3]]

[1] 1 0 0 1 1

[[4]]

[1] 0 1 1 0 1

[[5]]

[1] 1 0 1 0 1

[[6]]

[1] 1 1 0 0 1

[[7]]

[1] 1 0 1 1 0

[[8]]

[1] 1 1 0 1 0

[[9]]

[1] 0 1 1 1 0

[[10]]

[1] 1 1 1 0 0

The first combination is the one we used already (the real one, where site_1 and site_2 are upstream). Let's run the analysis again, but for each of the 10 possibilities. Ideally only the first (real) combination gives significant results.

meta$site <- substr(meta$site, 5, 5) # strip site name to just number, serves double duty as index into downstream_options combinations

combo_qvals <- function(combo){meta$downstream <<- apply(meta, 1, function(x){combo[as.numeric(x[3])]}); pvals <<- c(); for(i in 1:length(abundance[,1])){pvals <<- append(pvals,wil(i))}; q_random <- p.adjust(pvals, method="fdr"); return(length(which(q_random < 0.05)))}

result <- sapply(downstream_options, combo_qvals)

result

[[1]]

[1] 1073

[[2]]

[1] 0

[[3]]

[1] 0

[[4]]

[1] 0

[[5]]

[1] 0

[[6]]

[1] 0

[[7]]

[1] 0

[[8]]

[1] 0

[[9]]

[1] 0

[[10]]

[1] 0

*mic drop* [for real this time]

February 13, 2017

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.

Blog Archive

March 1, 2017

This is not normal: dealing with outliers, correlated factors, and skewed factor effects in RNASeq (a Kallisto/Sleuth/DESeq2/U-Test/skewness mash-up)

Stage 1: Identifying outlier samples

Stage 2: Spurious factor correlations

Stage 3: Rationally defining a factor level threshold for a skewed factor effect (numeric factor -> two-level factor)

February 16, 2017

RNASeq mic drop: U Test for detecting differential expression when factor dosage is unknown

February 13, 2017

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