Hi, Mike: Thanks a lot for the prompt response and input, which is very helpful Since some of the genes seem a bit interesting to us, and of course we love to keep. However, when I try: > resType <- results(dds, "Type_Tumor_vs_Normal",cooksCutoff=FALSE); Error in results(dds, "Type_Tumor_vs_Normal", cooksCutoff = FALSE) : unused argument (cooksCutoff = FALSE) I got error as above, I check the function results(), I did not see cooksCutoff = FALSE option there though, any advice? Also from your experience, if Cook's filtering is taken out, validation rate much worse? In reality, some genes might have large variation than others such as cancer-related genes. What do you think? Thanks and best Ming From: michaelisaiahlove@gmail.com Date: Mon, 10 Feb 2014 11:44:48 -0500 Subject: Re: DESeq2 - input data questions To: yi02@hotmail.com CC: bioconductor@r-project.org hi Ming, From the DESeq2 vignette, we have: "The results for particular genes can be set to NA, for either one of the following reasons: ... 2. If a row contains a sample with an extreme count then the p-value and adjusted p-value are set to NA. These outlier counts are detected by Cook’s distance. Customization of this outlier filtering is described in Section 3.3, along with a method for replacing outlier counts and refitting." ​Since the p-value and adjusted p-value are set to NA, and the base mean is not 0, you know that these genes were filtered based on the Cook's distances for the samples for these genes. ​See the practical and theoretical sections of the vignette as well as the manual page for ?results for more background. ​Cook's distance for a sample is the distance the model parameters would move ​if the sample were removed. From first glance, this makes sense too, in that you have shown genes with counts that have high variance: "...963 1642706 10488 1032959..." If the Cook's filtering is undesirable for your analysis you can disable it with: res <- results(dds, cooksCutoff=FALSE) Then these p-values will no longer be set to NA. So why did the genes get filtered by Cook's distance in the pooled analysis but not the paired? Cook's distance uses the residuals from the observed counts to the fitted values. When you include the subject variable, the counts have smaller residuals because more of the variance can be explained by the subject differences. In addition, DESeq2 requires 3 replicates per "cell", i.e., unique combination of covariates, in order to determine a count as an outlier based on Cook's distance. Mike On Mon, Feb 10, 2014 at 10:31 AM, Ming Yi <yi02@hotmail.com> wrote: Hi, Mike: I aksed you some questions about DESeq2 input data about a week or so as shown below. Then based on some other's experience and suggestions, I did use the rounded RSEM raw counts to run DESeq2, due to the dataset has matched tumors vs normals, I did try to use either or paired test to derive the DEGs. But I noticed a strange observation: For 7 genes of our interest, if I used paired test, I got the following: entrezID symbols gene_id baseMean log2FoldChange lfcSE 6440 6440 SFTPC SFTPC|6440 560695.51884 -4.309999 0.4811095 7476 7476 WNT7A WNT7A|7476 210.61484 -2.609732 0.3639806 8328 8328 GFI1B GFI1B|8328 11.25229 -2.445345 0.3040805 27255 27255 CNTN6 CNTN6|27255 297.51966 -5.049929 0.3276390 89780 89780 WNT3A WNT3A|89780 295.59498 -5.037872 0.2976513 256815 256815 C10orf67 C10orf67|256815 84.88613 -4.841887 0.2890447 388939 388939 C2orf71 C2orf71|388939 51.81894 -4.567162 0.3601929 pvalue padj 6440 3.292485e-19 5.363882e-18 7476 7.501115e-13 6.306398e-12 8328 8.854943e-16 1.023742e-14 27255 1.336693e-53 4.232359e-51 89780 2.922717e-64 2.206764e-61 256815 5.536370e-63 3.881589e-60 388939 7.655804e-37 6.679604e-35 However, if I used pooled test, I got the following: entrezID symbols gene_id baseMean log2FoldChange lfcSE 6440 6440 SFTPC SFTPC|6440 560695.51884 -3.654053 0.4720463 7476 7476 WNT7A WNT7A|7476 210.61484 -1.826788 0.3769339 8328 8328 GFI1B GFI1B|8328 11.25229 -1.822568 0.3423880 27255 27255 CNTN6 CNTN6|27255 297.51966 -4.103604 0.3509631 89780 89780 WNT3A WNT3A|89780 295.59498 -4.424951 0.3262553 256815 256815 C10orf67 C10orf67|256815 84.88613 -4.293802 0.3183049 388939 388939 C2orf71 C2orf71|388939 51.81894 -4.068668 0.3828477 pvalue padj 6440 NA NA 7476 NA NA 8328 NA NA 27255 NA NA 89780 NA NA 256815 NA NA 388939 NA NA Noticed that when I used paired test, the p-values of these genes look quite good, whereas if used pooled test, the pvalues are all NAs (most of other genes in dataset seem OK with valid p-values at least not NA values). almost all commands are the same for these two settings, except for the modeling step: Pooled test: dds <- DESeqDataSetFromMatrix(countData = countD,colData = colD,design = ~Type); Paired test: dds <- DESeqDataSetFromMatrix(countData = countD,colData = colD,design = ~Subject+Type); The Subject holds the sample names with matched tumor and normal samples. I did pull out the counts data for these genes, which seem not uncommon (i posted below FYI). T is tumors, N is normals,T4626_01A mactched with N4626_11A for the sample patient as 4626, and so on so forth. any idea why caused this, any issue with my data or bug(s) in the package? Here is the count data for these genes: entrezID T4626_01A N4626_11A T2668_01A N2668_11A T6145_01A N6145_11A 6440 6440 938904 1271864 963 1642706 10488 1032959 7476 7476 117 252 53 177 8 260 8328 8328 4 30 1 19 2 9 27255 27255 399 425 6 461 8 473 89780 89780 276 500 23 310 10 547 256815 256815 53 175 3 98 4 141 388939 388939 7 227 3 202 3 116 T6146_01A N6146_11A T6776_01A N6776_11A T6777_01A N6777_11A T6744_01A 6440 46915 610409 39544 855546 5535 658155 52865 7476 65 207 2 253 2 191 81 8328 3 9 2 13 2 11 3 27255 81 329 3 722 27 351 47 89780 34 330 8 555 3 465 31 256815 7 70 1 199 2 53 5 388939 1 22 1 73 18 116 9 N6744_11A T5932_01A N5932_11A T5933_01A N5933_11A T5936_01A N5936_11A 6440 667035 21596 690368 11 880241 184 858491 7476 235 29 258 1 240 55 323 8328 7 1 5 0 10 0 6 27255 501 4 545 0 676 9 674 89780 361 9 457 0 409 1 565 256815 85 2 97 0 175 6 173 388939 98 4 84 0 31 0 106 T6828_01A N6828_11A T6836_01A N6836_11A T4632_01A N4632_11A T2662_01A 6440 44239 737756 68 529280 2859 485130 13 7476 6 257 126 192 2 94 3 8328 0 14 1 9 0 10 2 27255 12 626 0 413 4 176 2 89780 11 587 8 214 0 14 3 256815 4 217 3 143 1 16 3 388939 6 94 1 57 2 0 0 N2662_11A T6147_01A N6147_11A T6743_01A N6743_11A T6829_01A N6829_11A 6440 1526175 16149 892008 9150 1138397 55 824547 7476 794 4 281 3 317 20 174 8328 19 5 16 2 5 5 8 27255 591 4 288 7 313 232 327 89780 1173 50 710 12 402 1 320 256815 254 68 330 0 114 18 80 388939 50 4 130 1 77 0 326 T6831_01A N6831_11A T4627_01A N4627_11A T4490_01A N4490_11A T6595_01A 6440 294 306899 8082 708377 115 1293737 81 7476 0 179 239 240 209 211 511 8328 1 12 6 37 3 5 2 27255 3 418 6 392 2 430 2 89780 1 333 67 488 7 298 8 256815 3 122 0 89 0 174 8 388939 0 10 0 63 1 31 0 N6595_11A T6835_01A N6835_11A T2657_01A N2657_11A T4676_01A N4676_11A 6440 1145375 41647 367263 334072 1082719 42168 959756 7476 728 283 196 10 244 18 282 8328 84 7 7 18 23 2 21 27255 975 88 348 13 448 19 606 89780 1479 31 185 145 631 23 383 256815 310 12 124 1 77 6 146 388939 82 1 6 1 108 6 31 T2661_01A N2661_11A T4625_01A N4625_11A T6745_01A N6745_11A T2665_01A 6440 413649 844142 202136 2070763 76113 1441306 82 7476 34 510 137 387 323 360 247 8328 72 32 4 44 3 8 6 27255 24 718 11 509 19 484 5 89780 15 1021 18 531 14 398 8 256815 17 121 1 152 8 192 1 388939 4 90 2 267 12 51 13 N2665_11A T2655_01A N2655_11A 6440 2200556 278151 1687880 7476 572 117 1096 8328 12 2 51 27255 895 89 921 89780 998 65 1490 256815 167 18 220 388939 123 88 142 Thanks in advance! Ming From: michaelisaiahlove@gmail.com Date: Wed, 29 Jan 2014 13:01:40 -0500 Subject: Re: DESeq2 - input data questions To: yi02@hotmail.com CC: bioconductor@r-project.org hi Ming, We do not recommend using rounded estimated values of read counts with DESeq2 (although I found an email from myself to the list one year ago contradicting this in the case that someone had no access to the raw data)​. Counts of reads which are proportionally assigned to genes and then rounded can be a bad fit for distributions like the Negative Binomial (and Poisson). For example, this procedure could generate values arising from a distribution that has variance less than the mean. As a rule, and to help avoid erroneous results, users should produce a matrix containing integer counts of reads uniquely aligned to features. Mike On Wed, Jan 29, 2014 at 11:16 AM, Ming Yi <yi02@hotmail.com> wrote: Hi, Dear Michael: Sorry to directly jump in a quick question, since I saw you respond to some DESeq questions yesterday and I know you are the author of DESeq2 package. I run into some RNAseq data issue with edgeR package (see my post in BioC yesterday), since the data I used are not integer counts of genes, but RSEM processed raw_count from TCGA RNA-seq data, in fact, are not integers, although TCGA named it as raw_count (I showed some piece of such a file from TCGA as below). file name: "unc.edu.1b4bb160-191b-4796-8759-339b11fe386d.1096727.rsem. genes.results" its content looks below: gene_id raw_count scaled_estimate transcript_id 1 ?|100130426 0.00 0.000000e+00 uc011lsn.1 2 ?|100133144 4.67 1.794813e-07 uc010unu.1,uc010uoa.1 3 ?|100134869 15.33 4.271899e-07 uc002bgz.2,uc002bic.2 4 ?|10357 218.79 1.933490e-05 uc010zzl.1 5 ?|10431 1255.00 5.033411e-05 uc001jiu.2,uc010qhg.1 6 ?|136542 0.00 0.000000e+00 uc011krn.1 I used the raw_count of this file for edgeR analysis and run into some warning messages which indicating potential issues, and had been pointed out by Dr. Smyth as you may have seen in my post yesterday. Reading the user guide from both DESeq and DESeq2, both packages mentioned the input was expected to be count data in a form of a matrix of integer values, not even (rounded) normalized counts, which would lead to nonsensical results. The TCGA data above apparently are not true raw counts but some kind of RSEM processed raw_counts. Although both RSEM authors and TCGA data staff all suggested the raw_count can be used for edgeR and DESeq (or DESeq2) analysis, because of the issues I run into in edgeR. I wonder what would be your opinion on this data for analysis in DESeq2? Any advice would be appreciated very much! Thanks and best Ming Yi NIH Maryland,USA On Wed, Jan 29, 2014 at 11:16 AM, Ming Yi <yi02@hotmail.com> wrote: Hi, Dear Michael: Sorry to directly jump in a quick question, since I saw you respond to some DESeq questions yesterday and I know you are the author of DESeq2 package. I run into some RNAseq data issue with edgeR package (see my post in BioC yesterday), since the data I used are not integer counts of genes, but RSEM processed raw_count from TCGA RNA-seq data, in fact, are not integers, although TCGA named it as raw_count (I showed some piece of such a file from TCGA as below). file name: "unc.edu.1b4bb160-191b-4796-8759-339b11fe386d.1096727.rsem. genes.results" its content looks below: gene_id raw_count scaled_estimate transcript_id 1 ?|100130426 0.00 0.000000e+00 uc011lsn.1 2 ?|100133144 4.67 1.794813e-07 uc010unu.1,uc010uoa.1 3 ?|100134869 15.33 4.271899e-07 uc002bgz.2,uc002bic.2 4 ?|10357 218.79 1.933490e-05 uc010zzl.1 5 ?|10431 1255.00 5.033411e-05 uc001jiu.2,uc010qhg.1 6 ?|136542 0.00 0.000000e+00 uc011krn.1 I used the raw_count of this file for edgeR analysis and run into some warning messages which indicating potential issues, and had been pointed out by Dr. Smyth as you may have seen in my post yesterday. Reading the user guide from both DESeq and DESeq2, both packages mentioned the input was expected to be count data in a form of a matrix of integer values, not even (rounded) normalized counts, which would lead to nonsensical results. The TCGA data above apparently are not true raw counts but some kind of RSEM processed raw_counts. Although both RSEM authors and TCGA data staff all suggested the raw_count can be used for edgeR and DESeq (or DESeq2) analysis, because of the issues I run into in edgeR. I wonder what would be your opinion on this data for analysis in DESeq2? Any advice would be appreciated very much! Thanks and best Ming Yi NIH Maryland,USA [[alternative HTML version deleted]]