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 at 16:13
  • $\begingroup$ Sorry I don't know how to do that $\endgroup$ – Fereshteh Jan 18 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 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 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 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 >=.

$\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 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 ]
$\endgroup$
  • $\begingroup$ Thank you, I added output of what you suggested in my post $\endgroup$ – Fereshteh Jan 21 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$ – Fereshteh Jan 21 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 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 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 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.