Filtering before the linear modelling and empirical bayes steps
2
0
Entering edit mode
@shamim-sarhadi-9395
Last seen 4.6 years ago

I apologize if this is a duplicate question,I read limmas paper (filtering section) it referred that after normalizing microarray data, for downstream analysis, it is usually worthwhile to remove probes that appear not be expressed but on the other hand in this post Non-specific filtering methodogies for ExpressionSet in R/Bioconductor , it referred that, filtering is not much needed if we use limma for microarray analysis , now if you think it is better to apply filtering before the linear modelling and empirical Bayes steps, please let me know what kind of filtering I can use without conflict with distribution of variance ? 

DEG Limma Filtering • 1.0k views
ADD COMMENT
2
Entering edit mode
@ryan-c-thompson-5618
Last seen 11 months ago
Scripps Research, La Jolla, CA

I think it is reasonable to filter out probe sets that can be confidently identified as non-expressed before the eBayes step. The purpose of eBayes is to squeeze the variance of each probe set toward the average variance of all probe sets, under the assumption that this average is a reasonable representation of the overall behavior of the data. Ideally, that average should only reflect the biology of your system, not the technical noise in probes that aren't measuring anything. So removal of non-expressed probes should help ensure that. However, it might make little difference in practice.

One thing you definitely want to avoid is to filter probe sets based on their variance. Doing so would significantly bias the eBayes estimate of the average variance of the data set and compromise your entire analysis.

ADD COMMENT
0
Entering edit mode
svlachavas ▴ 780
@svlachavas-7225
Last seen 6 hours ago
Germany/Heidelberg/German Cancer Resear…

Dear Beyondmind65,

actually the approach of "non-specific intensity filtering", is beneficial prior using limma, in order to remove non-interesting probes(i.e. probes that not expressed in the majority of your samples, & also to reduce the number of multiple hypothesis tests). On the other hand, is also data-dependent and it's up to you if and how you could filter your data prior statistical inference.

You could check also some other older posts:

A: Non specific filtering prior to differential expression with limma based on vari

If you have more details about the data/technology you are using and your experimental design, would be also more helpful than some general suggestions.

Best,

Efstathios-Iason

ADD COMMENT
0
Entering edit mode

Thank you Dear svlachavas for useful comment, I read the link, I found very helpful points to consider,I am doing a DEG analysis, my study is a class comparison between to condition, normal vs cancer on Affymetrix Human Genome U133 Plus 2.0 Array and Affymetrix Human Genome U133A Array, as you said, I think removing of non-informative probs could be helpful for reducing the number of multiple hypothesis tests
I read the link and I found very intresting notes that you and Dr. Gordon Smyth have mentioned and based of that I changed my scripts, it's kind of you if let me know about any problem or advice

design <- model.matrix(~factor(test$Disease))
fit <- lmFit(test, design)
keep <- fit$Amean > median(fit$Amean)
fit <- eBayes(fit[keep,], robust=TRUE, trend=TRUE)
tab <- topTable(fit, coef=2, adjust="fdr", n=150)

 

all the best

ADD REPLY
1
Entering edit mode

Well, i personally believe that filtering on the A values from above would be too strigent, as you could possibly loose informative probesets, that are lowly expressed in only one of the two groups--except you have already performed DE analysis, and you get very few genes with FDR<5%--or, from an MDS plot or a PCA you don't "see" relative differences between your groups of comparison. Actually, as i saw that you use the affymetrix hgu133a & plus2 arrays, my suggestion (there are still many alternative ways to consider), is to use the nice R package panp (https://www.bioconductor.org/packages/release/bioc/html/panp.html), which has been created for these two specific platforms. In detail, you could create present/absent calls for each of your probesets in each sample (that is if a probeset is expressed or not in a given sample), and then you could follow very simply the following implementation:

 

# let's say for instance that each group has 10 samples

sel <- rowSums(mat=="P") >10 # where mat the matrix output of panp-see page 5 of the vignette

eset.filtered <- eSet[ind,] # where eSet your normalized expressionSet

i hope that helps,

Efstathios

ADD REPLY
0
Entering edit mode

ّّI'll try with panp, Thank you alot

ADD REPLY
0
Entering edit mode

Hi Efstathios, first I thank you for your nice offer about panp package,

but I have a question about using this package, in your above scripts, you said that for instance each group has 10 samples, I don't know your mean about each group... !! in my case I have 3 GSE that have 100, 120 and 130 samples ,I think again I need a cutoff to specify the number of samples !!? please let me know if there is an important note about this case

Sincerely

Shamim Sarhadi

ADD REPLY
1
Entering edit mode

Sorry to reply with delay. I just proposed a simply example of comparison groups, with equal sample size each. So, in your case, which is the sample size of cancer and control samples, respectively ? I assume you analyze each dataset, separately, right ? 

ADD REPLY
0
Entering edit mode

Hi,Thanks for your reply

Let me explain my case more clearly

I have 3 datasets

The first one>> 20 sample control+ 30 sample ductal carcinoma insitu(DCIS)+30 sample invasive ductal carcinoma(IDC)=100

The second one>>50sample(DCIS) + 70 sample(IDC)=120

The third one>> 30 sample control+ 45 (DCIS)+55 sample (IDC)=130

I want to merge these datasets with combat methode and then DEG analysis(DCIS vs control and IDC vs control)

From the previous post I found that if I use non-spesific filtering I will give a better output from my analyze, so as you suggest me, I decide to perform pa.calls function on my dataset I read its vignette but as I mentioned above I am confused with selecting the correct sample size

Thank you alot

Sincerely

ADD REPLY
0
Entering edit mode

So, you have already performed merging of your datasets and ComBat, right ? So, you have a unified dataset: if so, you could very "naively", set a cut-off number of samples by each group-however, as i see from above, the number of samples in each group varies and it is not the same. So, you could "arbitary" choose a number of samples to apply the cutoff--or alternatively (which is a bit more complicated), you could utilize the matrix of the P/A values from PANP, and keep the probesets which are present in half of the samples in at least one condition. The last, sounds more appropriate to me.

ADD REPLY
0
Entering edit mode

excuse me for these noob questions, I tried to perform pa.calls function on my merged eSet but it gave me an error 

PAPA<-pa.calls(test1)
Error in AllExprs[NSMPnames, ] : subscript out of bounds
In addition: Warning messages:
1: In if (chip == "hgu133b") { :
  the condition has length > 1 and only the first element will be used
2: In if ((chip != "hgu133a") & (chip != "hgu133atag") & (chip != "hgu133plus2")) { :
  the condition has length > 1 and only the first element will be used
3: In if ((chip == "hgu133a") | (chip == "hgu133atag")) { :
  the condition has length > 1 and only the first element will be used

also I don't know how can I deal with those esets that have more than one group, in your above codes you have supposed that each eSet has one covariate i.e. cancer or control, for example in this code

sel <- rowSums(mat=="P") >10

if the eset has the samples from case and control,this logical condition perform on both of group (case and control) I want to know there is any problem if I use this code for my dataset that has sample from case and control?

Thanks in advance

all the best

ADD REPLY
1
Entering edit mode

What platforms your individual datasets comprise of ? If they don't have the same, i suspect that this is the reason for the above error. If you type test1, it would return a slot called Annotation:

so, there you will probably have more that than one annotation platforms. Thus, if my notion is correct, then you will have to use a different kind of filterning, because PANP can't worked with merged eSets comprising of multiple annotations. In this case, you could consider something similar:

ind <- rowSums(exprs(test1)> expression.cutoff) >=N # where in the expression.cutoff you could use something like 6, and the number of samples is up to you.

 

Alternatively, if you don't have specific types of comparison you want to perform that the merging of datasets would be beneficial to increase sample size, you should also consider analyzing first each dataset separately.

Best,

Efstathios

ADD REPLY
0
Entering edit mode

Hi Efstathios

I am so grateful for your thorough answer, I tried with panp on each eset seperatly and it works,thanks for your answer

 

Best

ADD REPLY

Login before adding your answer.

Traffic: 870 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