1
$\begingroup$

I conducted a differential expression analysis over several datasets, using LIMMA, each one on its own. For each dataset, I have a data frame of all the genes, with the LFC and the p-value. Those data frames are called deg1, deg2, deg3 ... and so on.

I want to use the MetaVolcanoR package, using the votecount_mv approach since it is the only approach that I can pick the p-value and LFC thresholds that I desire.

The problem is I'm not sure how to tune the parameters. There are several problems that I'm facing and I would really use some help please:

How do I set the log fold change threshold if I want only LFC > 0.5 and LFC < -0.5? Is it log fold change = 0.5 ? just making sure.

What is the metathr parameter? I read about it and still don't understand what it does. I set it as 0.01 for now.

And the main problem, the plot is showing much more gene number than there actually is. You can see in Dataset6 for example, the plot says it has over 50,000 genes while in reality, it has something like 18,000 genes. The same for Dataset8, why is this happening?

This is the plot:

enter image description here

This is my code:

totalDEG = list(Dataset1 = deg1, Dataset2 = deg2,
                Dataset3 = deg3, Dataset4 = deg4, Dataset5 = deg5, Dataset6 = deg6, 
                Dataset7 = deg7, Dataset8 = deg8 ,Dataset9 = deg9, Dataset10 = deg10,
                Dataset11 = deg11, Dataset12 = deg12, Dataset13 = deg13 , Dataset15 = deg15,
                Dataset16 = deg16, Dataset20 = deg20, Dataset21 = deg21)

totalDEG = map(totalDEG, ~ .x %>% rownames_to_column("symbol") %>% `row.names<-`(.$symbol))

meta_degs_vote <- votecount_mv(diffexp=totalDEG,
                               pcriteria='P.Value', 
                               foldchangecol='logFC',
                               genenamecol='symbol',
                               geneidcol=NULL,
                               pvalue = 0.05,
                               foldchange = 0.5, 
                               metathr=0.01,
                               collaps=FALSE,
                               jobname="MetaVolcano", 
                               outputfolder=".",
                               draw='HTML')

head(meta_degs_vote@metaresult, 50)
meta_degs_vote@degfreq

This is a small sample of the deg13 dataframe, just as an example of how it looks:

structure(list(logFC = c(3.39094029649578, 2.91185452846888, 
-4.18752593306854, 4.15866027194434, -3.46058235630518, -2.87208933878881, 
-3.08172118179278, -2.7458737595494, -4.086372553458, -2.68889593085408, 
-3.03835320252487, 2.61330670793883, -2.48123655627091, 1.49263378072949, 
1.95754205818365, -2.70230318294949, -2.9154447157702, -3.38110311321268, 
-2.35284609598761, -3.61327208834916, 2.03692346164301, -1.70931971282814, 
2.27925549238084, -3.79524834465106, 1.64515339181668, -4.30959811283984, 
3.57930610863801, -2.54394115096987, -1.98143797060853, 4.03178891878422
), AveExpr = c(6.4934411845873, 6.79507589173635, 3.87626322696796, 
5.02990511836996, 3.37154863566726, 4.54352970704667, 2.51617183960423, 
5.36760482831628, 8.74162861573995, 3.52708473144296, 1.56671205994866, 
5.87684931921495, 4.44218678637856, 5.12219572880673, 7.17539095745285, 
3.53733757517713, 3.35148844504943, 5.15160050823647, 4.75910049955374, 
2.76067393839851, 7.59404366884022, 3.26758602409262, 6.6742932845645, 
2.02865432630183, 7.35343444625846, 1.63566796485655, 6.63046957559406, 
4.0296762052982, 5.52320581646106, 5.21519413005468), t = c(7.67283765970783, 
7.48874311179818, -6.82705440075341, 6.17183091062419, -6.04487485029223, 
-5.8487763400982, -5.92000790680834, -5.64958317461715, -5.61111489307273, 
-5.58440722941211, -5.78651741270403, 5.34175048353288, -5.31088208376663, 
5.08100365503679, 5.03973023175764, -5.07773225272668, -5.04632997176375, 
-4.96736474738699, -4.82583305184561, -4.92525789154973, 4.78892270501915, 
-4.75294466621978, 4.64502680172295, -4.88600418339975, 4.62243262623093, 
-4.78866983562446, 4.56376665911121, -4.57063797766699, -4.48421852802591, 
4.47493043300733), P.Value = c(8.61557824744483e-06, 1.08572492426076e-05, 
2.57808593137873e-05, 6.41321288302935e-05, 7.70283126846018e-05, 
0.000102673342223582, 9.24388613776644e-05, 0.000138234657017388, 
0.000146500473045098, 0.000152546919410067, 0.000112607650587545, 
0.000221306832937214, 0.000232172444273951, 0.000333174926637586, 
0.000355778600664596, 0.000334909952712188, 0.000352058129044564, 
0.000399405716020234, 0.000501862041539822, 0.000427357183393412, 
0.000532900955175219, 0.000565105101377844, 0.000674574457630558, 
0.000455275577064366, 0.000700192843853878, 0.000533120436605821, 
0.000771595292301189, 0.000762850982694916, 0.00088083209108589, 
0.000894604612049976), adj.P.Val = c(0.105689892752164, 0.105689892752164, 
0.167309183326708, 0.24749466449955, 0.24749466449955, 0.24749466449955, 
0.24749466449955, 0.24749466449955, 0.24749466449955, 0.24749466449955, 
0.24749466449955, 0.322868951254968, 0.322868951254968, 0.364560714544159, 
0.364560714544159, 0.364560714544159, 0.364560714544159, 0.388801494259897, 
0.397198753468022, 0.396200809689826, 0.397198753468022, 0.397198753468022, 
0.397198753468022, 0.397198753468022, 0.397198753468022, 0.397198753468022, 
0.397198753468022, 0.397198753468022, 0.397198753468022, 0.397198753468022
), B = c(3.51320342275366, 3.35460378796296, 2.44903607778273, 
1.62363531926393, 1.57839963166078, 1.44388223679157, 1.36898233186455, 
1.24528027435772, 1.14177006290097, 1.04403303682497, 0.945243757754988, 
0.859404654981856, 0.797031901979986, 0.496397382158138, 0.450148400219772, 
0.418557695883232, 0.3417613688148, 0.324284286390148, 0.136478287041671, 
0.0982449107580958, 0.0925437909507512, -0.0280517206307875, 
-0.112494849608916, -0.129424832005619, -0.147717311976667, -0.200056933894112, 
-0.234222113297706, -0.244412573722838, -0.351160223852428, -0.386270016806261
)), row.names = c("FOS", "DUSP1", "TP53AIP1", "FOSB", "TFAP2C", 
"KIF26A", "GRHL3", "CTAGE4", "FCGBP", "FBXO27", "AC087289.3", 
"SH3BP5", "GABRE", "ADAMTS4", "DDHD1", "KCNS3", "ARHGEF16", "SERINC2", 
"OPLAH", "WNT7B", "GSN", "OVOL2", "CD163", "SVOPL", "ZFP36", 
"SPESP1", "FGD5", "FAAH", "ERBB3", "PARM1"), class = "data.frame")
$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.