problem while retriving the genes by using p value and log2fc as conditions
3
0
Entering edit mode
@cbd9b5e7
Last seen 12 months ago
India

I am performing differential expression analysis with DESeq2. After running dds<- DESeq(dds) i wanted to retrive all the up and down regulated genes. i gave the conditionpvalue < 0.05 and log2FoldChnage > 1.5 for up regulated genes. But in the result file also include the genes which have pvalue > 0.05. iam pasting my example dataset here with the result. I wnted to know that which command will be helpfull here for retriving genes which satisfy the above mentioned two condition.

Thank you in advance.

Code should be placed in three backticks as shown below


> dds <- makeExampleDESeqDataSet()
> dds
class: DESeqDataSet 
dim: 1000 12 
metadata(1): version
assays(1): counts
rownames(1000): gene1 gene2 ... gene999 gene1000
rowData names(3): trueIntercept trueBeta trueDisp
colnames(12): sample1 sample2 ... sample11 sample12
colData names(1): condition
> dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> res <- results(dds, alpha =0.05)
> res
log2 fold change (MLE): condition B vs A 
Wald test p-value: condition B vs A 
DataFrame with 1000 rows and 6 columns
          baseMean log2FoldChange     lfcSE       stat    pvalue      padj
         <numeric>      <numeric> <numeric>  <numeric> <numeric> <numeric>
gene1     12.61910      -0.320613  0.728490  -0.440107  0.659860  0.992529
gene2      6.01265      -1.391883  0.926580  -1.502172  0.133053  0.972932
gene3     71.48063      -0.723134  0.319153  -2.265794  0.023464  0.805855
gene4      7.56355      -0.347392  0.715907  -0.485248  0.627501  0.992529
gene5     11.86764       0.275620  0.634890   0.434123  0.664199  0.992529
...            ...            ...       ...        ...       ...       ...
gene996   27.95109     -0.4044216  0.420641 -0.9614415  0.336330  0.992529
gene997   41.27173     -0.0273015  0.341122 -0.0800344  0.936210  0.996243
gene998   32.97030      0.1177037  0.357531  0.3292127  0.741995  0.992529
gene999    2.85416     -0.4862333  1.009720 -0.4815524  0.630124  0.992529
gene1000   3.24778     -0.2205322  1.003172 -0.2198348  0.826000  0.992529
> RESULT <- res[which(res$pvalue<0.05) & abs(res$log2FoldChange > 1.5),]
Warning message:
In which(res$pvalue < 0.05) & abs(res$log2FoldChange > 1.5) :
  longer object length is not a multiple of shorter object length
> write.table(RESULT, file = "upregulated_test.txt")
> RESULT
log2 fold change (MLE): condition B vs A 
Wald test p-value: condition B vs A 
DataFrame with 30 rows and 6 columns
         baseMean log2FoldChange     lfcSE      stat     pvalue      padj
        <numeric>      <numeric> <numeric> <numeric>  <numeric> <numeric>
gene20    7.02072        1.88494  0.722398   2.60928 0.00907324  0.805855
gene45    4.76153        1.70103  0.877075   1.93943 0.05244909  0.851959
gene57    0.88648        3.34660  2.213826   1.51168 0.13061522  0.972932
gene74    1.27885        2.14873  1.632143   1.31651 0.18800332  0.992529
gene75    4.16297        1.79519  1.024415   1.75240 0.07970449  0.901962
...           ...            ...       ...       ...        ...       ...
gene868   1.34768        3.38540   1.50091   2.25557  0.0240976  0.805855
gene923   1.36849        1.68479   1.60413   1.05028  0.2935898  0.992529
gene932   6.63682        1.78308   0.86071   2.07164  0.0382986  0.851959
gene970   2.86467        1.79572   1.32983   1.35034  0.1769079  0.992529
gene977   4.02625        1.50587   1.01119   1.48920  0.1364334  0.974524

Iam also pasting the attaching the file which i created using above mentioned command.

Bioconductor Transcriptomics DifferentialExpression DESeq2 • 809 views
ADD COMMENT
0
Entering edit mode
@james-w-macdonald-5106
Last seen 21 hours ago
United States

Did you see this part of your post?

> RESULT <- res[which(res$pvalue<0.05) & abs(res$log2FoldChange > 1.5),]
Warning message:
In which(res$pvalue < 0.05) & abs(res$log2FoldChange > 1.5) :
  longer object length is not a multiple of shorter object length

If you get a warning message, it's usually a good idea to try to understand why you are getting it. The warning message provides a clue, saying the longer object is not a multiple of the shorter object. That should clue you in to the fact that R is recycling a vector, and the recycled vector isn't an even fraction of the longer vector. Having R recycle vectors is almost always not what you want.

The reason those two vectors are not the same length is because wrapping a boolean vector in which, will only return the location of the TRUE values, which is not what you want. Instead, you want two boolean vectors of the same length.

> RESULT <- res[res$pvalue<0.05 & abs(res$log2FoldChange > 1.5),]

Is what you want.

ADD COMMENT
0
Entering edit mode

Also needs to wrap abs() only around the LFC.

ADD REPLY
0
Entering edit mode
@cbd9b5e7
Last seen 12 months ago
India

I tried this command but it giving some error

Error: logical subscript contains NAs

ADD COMMENT
0
Entering edit mode
@cbd9b5e7
Last seen 12 months ago
India

Thank you so much. I solved this error using is.na() and also i got the proper output using the above mentioned command. Thank you

ADD COMMENT

Login before adding your answer.

Traffic: 868 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6