[rfilbert]

I see you aren't getting it.

Genes with very few reads, of which there can be hundreds or thousands, are unlikely to be called as differentially expressed, unless the differences are very large.

On the other hand, because of multiple testing correction, many genes of higher expression will not pass the threshold of being considered differentially expressed.

Filtering will result in a handful of genes with low expression that are differentially expressed being discarded, but will allow for many more genes of higher expression to pass the threshold of being differentially expressed.

So overall MORE genes are called differentially expressed and the overall number of false negatives is decreased.

You should take a day or two to actually read some papers on the matter. Like the linked PNAS paper, the DESeq, EdgeR, etc papers and vignettes. These all explain exactly what they do so that nothing is hidden, unlike the black boxes that you are putting your data into.

What black box am I putting my data into? BTW, I fully understand how FDR works and why a smaller starting list results in a less stringent FDR. Apparently, however, you have missed my point. Since you are so convinced on filtering out lowly expressed genes, exactly what is your threshold? I suppose your threshold is going to be sufficiently small that there is no possibility of filtering out something that's real.

[chadn737]

What black box am I putting my data into? BTW, I fully understand how FDR works and why a smaller starting list results in a less stringent FDR. Apparently, however, you have missed my point. Since you are so convinced on filtering out lowly expressed genes, exactly what is your threshold? I suppose your threshold is going to be sufficiently small that there is no possibility of filtering out something that's real.

The key to filtering is not that there is "no possibility of filtering out something that's real", but rather that your by filtering you observe more real differences than what you filter out.

Let me reiterate the basic concept which you seem intent on simply ignoring.

Genes with low read counts are less likely to be differentially expressed. So by retaining them, you will only identify a small number of differentially expressed genes that have the most significant differences.

By removing them, you do loose some real differences, but now you identify an even greater number of real differences in the retained data because of your multiple testing correction.

So its a tradeoff.

Lets say that by retaining all the values, you observe 5 real differences in those genes that would be lost to filtering. However, by filtering these out, you now observe 20 real differences that before were not called because of your multiple testing correction.

So by loosing 5, you gain 20

The big elephant in the room that you are ignoring is that you also loose real differences by NOT filtering...but you claim to understand how multiple testing correction works.

[rfilbert]

I understood all that the first time you explained it. However, you've just pulled the numbers 5 and 20 out of thin air - not from real data or even simulated data. Why not do something useful and answer my question? What criteria do you use to filter out low expressed genes? Please be specific.

[chadn737]

Let me give you a real example then.

Before filtering: 2034 differentially expressed genes (FDR adjusted p-value < 0.05)

Because I don't feel like spending a lot of time on this, as I am only trying to prove a point, I simply filtered out the bottom 25% of genes by their raw read counts. I then test for differential expression and find that I have 2134 genes that are differentially expressed.

Thats an increase in 100 genes detected and this was just a very simple and rough filtering, although the DESeq vignette provides an explanation and justification for this sort of filtering, as well as a more through approach.

Of course, if you would bother to read the linked PNAS paper or the vignettes of any of these programs, you would be more familiar with this.

[rfilbert]

First of all, I don't think it is wise or friendly of you to assume that if I don't agree with your approach, I have not read the literature or are somehow less learned than you are. I will ignore your personal attacks and just try to focus on a discussion that might be useful to others. From your post, I can only guess that you recommend filtering out the bottom 25% of genes, and I assume that is by their average expression across all samples? You still did not give any specific advice on what threshold you used.

As an example where I am concerned about such a drastic approach, consider a 2 factor design from say, gastric cancer which affects men and woman. The 2 factors are tumor status (tumor, normal) and gender (male, female). If a gene is only expressed in the tumor/female group, assuming a balanced design, it would only be expressed in only 25% of the samples. It's easy to give real world examples like this where the group comprises much less than 25% as well. If one just took the average expression across all samples, I worry this gene would be filtered out, and this is an example of the false negatives that can come from such a pre-filtering approach. My approach is to post-filter to avoid this kind of problem. I can always disgard results if I feel that there are not enough reads to be reliable, but I want to know that I did not discard something like the example I gave here.

[Jeremy]

What black box am I putting my data into? BTW, I fully understand how FDR works and why a smaller starting list results in a less stringent FDR. Apparently, however, you have missed my point. Since you are so convinced on filtering out lowly expressed genes, exactly what is your threshold? I suppose your threshold is going to be sufficiently small that there is no possibility of filtering out something that's real.

R and the packages therein (DESeq, EdgeR etc.) are free and open source, so it is not the black box that commercial software is (*cough*Partek*cough*). You can see the underlying code and know exactly what is happening with the data.

[Simon Anders]

Dear rfilbert

I've followed this thread now for a while and think that chadn737 deserves some support. I applaud him for the patience he has with you. Your discussion style is offensive and not on the level that we are used to in this forum.

Your orginal argument was that pre-filtering will always cause information loss and hence reduce inferential power, and that using it is a sign of a lack of grasp of statistics, or at least an unfounded fear of big data. Now, after suddenly agreeing with chadn737 that there are situations in which filtering does increase overall number of detected hits, and claiming that you had understood this from the beginning, you complain that chadn737'argument is invalid because he has not explained how to find a threshold that indeed increases rather than decreases power. Changing the question to avoid admitting that one was wrong is a tactics with which you get away in politics but not in science.

So, to reiterate the facts that chadn737 tried in vain to get you to understand: 1. Filtering out genes that are unlikely to yield in significant result can increase power. 2. Obviously, filtering has the opposite effect if the filter is too stringent. 3. Hence, one needs to use a sensible approach to chose the filtering criterion.

chadn737 suggested to use the methods suggested by Bourgon et al. in their PNAS paper. He did not bother to spell them out to you, because, as you emphasize, you of course have already read the paper carefully and are hence presumably familiar with them.

BTW, as an aside to other readers: The main point of the Bourgon et al. paper is not how to chose an effective filter (though this is covered), but rather how to decide whether a filter criterion is independent from the testing procedure and so does not invalidate the multiple-testing adjustment.

[rfilbert]

Partek is not a black box, and I am unaware of any proprietary algorithms in their software. When we were first considering different software packages for our microarray analysis, our statistician did a rigorous evaluation of Partek's statistics and concluded that they were very solid. He did this primarily by comparing the results of a variety of experimental designs to see if Partek results match those from SAS and they did. He views SAS as the gold standard for statistically accuracy (another commercial package). The problem I have with R is that anyone can contribute anything and you have no idea if it is giving the correct answer. And I don't think there are many people on this planet who can verify the correctness of a statistical algorithm simply by examining the source code. That person would have to be a theoretical statistician and an expert in computer science, and they would have to write plenty of test harnesses to examine all the corner cases.

[chadn737]

First of all, I don't think it is wise or friendly of you to assume that if I don't agree with your approach, I have not read the literature or are somehow less learned than you are. I will ignore your personal attacks and just try to focus on a discussion that might be useful to others. From your post, I can only guess that you recommend filtering out the bottom 25% of genes, and I assume that is by their average expression across all samples? You still did not give any specific advice on what threshold you used.

As an example where I am concerned about such a drastic approach, consider a 2 factor design from say, gastric cancer which affects men and woman. The 2 factors are tumor status (tumor, normal) and gender (male, female). If a gene is only expressed in the tumor/female group, assuming a balanced design, it would only be expressed in only 25% of the samples. It's easy to give real world examples like this where the group comprises much less than 25% as well. If one just took the average expression across all samples, I worry this gene would be filtered out, and this is an example of the false negatives that can come from such a pre-filtering approach. My approach is to post-filter to avoid this kind of problem. I can always disgard results if I feel that there are not enough reads to be reliable, but I want to know that I did not discard something like the example I gave here.

1) No, I do not recommend simply filtering out the bottom 25%. I made it abundantly clear in my previous post that I was taking an overly simplistic approach simply to prove a point.

2) You're still ignoring the elephant in the room. False negatives are produced by NOT filtering. There is no experimental design where false negatives are completely avoided. Pre-filtering may result in genes, like in your example being excluded, but not pre-filtering results in genes being excluded. I have explained that multiple times. Take the example from real data I gave. From the original 2034, there were differentially expressed genes excluded by pre-filtering. However, pre-filtering resulted in many more genes being identified as differentially expressed than were lost by pre-filtering.

[rfilbert]

Calling your argument an elephant in the room is just your way of insulting me. As I've already stated, I understand that by reducing the gene list, FDR becomes less stringent. You really don't need to explain it over and over again. Now, on to practical advice, how do you recommend one decides exactly what criteria is used to filter out the low expressed genes? In your example that resulted in 25% being filtered out, what was your cutoff and was it based on the mean of all samples, the max, or what?

[Jeremy]

Calling your argument an elephant in the room is just your way of insulting me. As I've already stated, I understand that by reducing the gene list, FDR becomes less stringent. You really don't need to explain it over and over again. Now, on to practical advice, how do you recommend one decides exactly what criteria is used to filter out the low expressed genes? In your example that resulted in 25% being filtered out, what was your cutoff and was it based on the mean of all samples, the max, or what?

I'm not claiming this is 'THE' way, but, using a de-novo transcriptome assembly from trinity (which is not uncommon to give in excess of 200,000 assembled sequences from an assembly of multiple samples) I chose to exclude sequences that had less than 500 reads in total from 12 samples. This threshold was chosen as it was the start of linear loss of sequences when increasing the threshold by increments of 100. The threshold resulted in about 80% of the sequences being removed while 95% of the mapped reads were in the retained sequences. So obviously the removed sequences were junk, either fragments representing transcriptional noise or fragments representing chimeric sequences from library design or some other sort of sequencing/assembly error.

My guess is this whole filter criteria argument stems from everybody having a different assumption of what the data looks like. But RNA-seq the data is highly variable, it can be very different to microarray expression data.