0
$\begingroup$

I have two normalized gene expression values (log2 of cpm)

> head(biomarker[,1:5])
               A2        A3        A4        A6        A7
A2M    12.4071618 12.631601 12.889748 11.842511 11.574134
ABCB11  0.5151957  6.176788  5.414905  7.134915  1.247590
ABCC2   6.7244303  6.487794 10.178615  6.384132  6.089462
ABCC6   6.3977122  5.823627  7.370761  7.397091  6.071587
ABCF1   9.1609847  9.893258 10.116638 10.225520 10.838486
ABCG2   7.6841874  7.738293  6.833152  6.888041  6.030032
>

And

> head(immune[,1:5])
              A2        A3        A4        A6        A7
A2M    10.904748 11.388404  9.910614 11.513439 12.963609
ABCB11  5.011380  6.359443  7.145992  8.451947  7.722605
ABCC2   5.040461  6.477014  3.873996  6.409777  9.133971
ABCC6   7.798441  7.601848  9.948072 11.628533 12.701460
ABCF1   8.553597  8.615425 11.145903 10.289098 11.444140
ABCG2   6.224294  5.629375  8.293416  7.979859  8.603793
> 

> dim(biomarker)
[1] 719  56
> dim(immune)
[1] 719  56
> 

I have matched samples in these experiments and same gene too. I could roughly say in both experiments everything has been constant. However, I want to do Student's t-Test to find genes that they had inconsistent gene expression patterns between the two data sets.

By this function I am getting a list but I don't know how to extract the names of genes with p-value > 0.05 or whatever cut-off

> f <- function(x,y){
+     test <- t.test(x,y, paired=TRUE)
+     out <- data.frame(stat = test$statistic,
+                       df   = test$parameter,
+                       pval = test$p.value,
+                       conl = test$conf.int[1],
+                       conh = test$conf.int[2]
+     )
+     return(out)
+ }
> t=sapply(seq(ncol(t(biomarker))), function(x) f(t(biomarker)[,x], t(immune)[,x]))

Any help please?

  str(t)

    [list output truncated]
     - attr(*, "dim")= int [1:2] 5 719
     - attr(*, "dimnames")=List of 2
      ..$ : chr [1:5] "stat" "df" "pval" "conl" ...
  ..$ : NULL
    >


> sapply(t, "[[", "pval")
Error in FUN(X[[i]], ...) : subscript out of bounds
>  


> pvalueGenes <- sapply(colnames(biomarker),function(i) t.test(biomarker[i, ],immune[i, ], paired = TRUE)$p.value)
Error in if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") : 
  missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In mean.default(x) : argument is not numeric or logical: returning NA
2: In if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") :

 Show Traceback

 Rerun with Debug
 Error in if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") : 
  missing value where TRUE/FALSE needed
$\endgroup$
  • 1
    $\begingroup$ Did you try iterate by colnames instead? Have you searched the internet about how to get names of the variables you iterate over them? BTW you might be interested in the broom package, function tidy. $\endgroup$ – llrs Jan 18 '19 at 16:13
  • $\begingroup$ Sorry I don't know how to do that $\endgroup$ – Exhausted Jan 18 '19 at 20:32
  • $\begingroup$ A simple google search would give you answers. You could use tibble::rownames_to_column(biomarker, col_name). You could even do biomarker$col_name <- colnames(biomarker) $\endgroup$ – Ram RS Jan 20 '19 at 18:07
  • 2
    $\begingroup$ If you're going to downvote, add a reason why. Don't be passive-aggressive. It's bad for the community. Thanks! $\endgroup$ – Alex Reynolds Jan 21 '19 at 20:30
  • 1
    $\begingroup$ @AlexReynolds It is not required to disclose, but I think that the downvote button already explain this: "does not show research effort, it is unclear or not useful". Asking for an explanation from downvoters is also passive-aggressive... $\endgroup$ – llrs Jan 22 '19 at 8:17
2
$\begingroup$

Starting from the definition of t:

...
> t <- sapply(seq(ncol(t(biomarker))), function(x) f(t(biomarker)[,x], t(immune)[,x]))

Set the column names to the gene names from one or the other of your conditions:

> colnames(t) <- rownames(immune)

Set the cutoff:

> cutoff <- 0.05

Filter for rows which meet the p-value cutoff criterium, and get their column names (gene names):

> colnames(t)[which(unlist(lapply(t["pval",], function(x) x < cutoff)))]
[1] "ABCC6" 
...

Only one gene name is reported here, using your demo datasets. Adjust your cutoff and criterium logic, as needed. For example, I filter here for gene names where the p-value is less than the cutoff; however, from your question, you may be asking for the inverse, for whatever reason, so you could change < to >=.

| improve this answer | |
$\endgroup$
  • 3
    $\begingroup$ I might suggest changing the variable name t to something else, as t() takes the transpose of the argument, if of matrix or data frame type. Using non-ambiguous variable names will make it easier for you to make changes or fix problems down the road. $\endgroup$ – Alex Reynolds Jan 21 '19 at 22:19
2
$\begingroup$

From your example we are comparing samples not genes, I think we need to run the test row wise, something like this:

pvalueGenes <- sapply(colnames(biomarker),
                      function(i) t.test(biomarker[i, ],
                                         immune[i, ], paired = TRUE)$p.value)

names(pvalueGenes)[ pvalueGenes < 0.05 ]
| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you, I added output of what you suggested in my post $\endgroup$ – Exhausted Jan 21 '19 at 13:02
  • $\begingroup$ Sorry, I have 2 matched matrices with matched samples and genes I want to find genes with p-value < 0.05 between 2 matrices but I could not do that $\endgroup$ – Exhausted Jan 21 '19 at 13:18
  • $\begingroup$ @FereshTeh forget your function f, there is no need, just run the commands in my answer. And check what you see in pvalueGenes object. $\endgroup$ – zx8754 Jan 21 '19 at 16:58
  • 1
    $\begingroup$ @FereshTeh The reason for error is statistics, read about it more here. Suggestion, when you get an error try to use the search, 99% of the time the error has been searched before and you will find the solution/explanation. $\endgroup$ – zx8754 Jan 21 '19 at 19:12
  • 1
    $\begingroup$ @RamRS Not sure what you mean by "logical unit". Yes both approaches work, but your version is about 1.14 times slower, tested on 100K vector. $\endgroup$ – zx8754 Jan 23 '19 at 11:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.