BioinformaticsFMRP / TCGAbiolinks

TCGAbiolinks
http://bioconductor.org/packages/devel/bioc/vignettes/TCGAbiolinks/inst/doc/index.html
295 stars 112 forks source link

get.pair from TCGAbiolinks not working. #435

Closed NikolayAlabi closed 2 years ago

NikolayAlabi commented 3 years ago

Hello,

The following is my script.

library(ELMER) library(parallel)

distal.probes <- get.feature.probe(genome = "hg38", met.platform = "450K")

library(MultiAssayExperiment) mae <- createMAE(exp = acc.exp, met = acc.met, save = TRUE, linearize.exp = TRUE, filter.probes = distal.probes, save.filename = "mae_thca.rda",met.platform = "450K", genome = "hg38", TCGA = TRUE) mae <- mae[,!mae$is_ffpe]

tgpx.mae <- as.character(read.csv("mae.assignments.csv", header=F)$V1) tgpx.mae <- as.factor(tgpx.mae) mae$tgpx <- tgpx.mae

sig.diff_lowvsint <- get.diff.meth(data = mae, group.col = "tgpx", group1 = "low", group2 = "int", minSubgroupFrac = 0.2, sig.dif = 0.3, diff.dir = "hypo", cores = 1, dir.out = "lowvsint", pvalue = 0.01) sig.diff_lowvshig <- get.diff.meth(data = mae, group.col = "tgpx", group1 = "low", group2 = "hig", minSubgroupFrac = 0.2, sig.dif = 0.3, diff.dir = "hypo", cores = 1, dir.out = "lowvshig", pvalue = 0.01) sig.diff_intvshig <- get.diff.meth(data = mae, group.col = "tgpx", group1 = "int", group2 = "hig", minSubgroupFrac = 0.2, sig.dif = 0.3, diff.dir = "hypo", cores = 1, dir.out = "intvshig", pvalue = 0.01)

nearGenes_lowvsint <- GetNearGenes(data = mae, probes = sig.diff_lowvsint$probe, numFlankingGenes = 20) nearGenes_lowvshig <- GetNearGenes(data = mae, probes = sig.diff_lowvshig$probe, numFlankingGenes = 20) nearGenes_intvshig <- GetNearGenes(data = mae, probes = sig.diff_intvshig$probe, numFlankingGenes = 20)

pair_lowvsint <- get.pair(data = mae, group.col = "tgpx", group1 = "low", group2 = "int", nearGenes = nearGenes_lowvsint, minSubgroupFrac = 0.4,permu.dir = "lowvsint/permu", permu.size = 100, raw.pvalue = 0.05, Pe = 0.01, filter.probes = TRUE, filter.percentage = 0.05, filter.portion = 0.3, dir.out = "lowvsint", label = "hypo")

Selecting U (unmethylated) and M (methylated) groups. Each groups has 20% of samples
-------------------
* Filtering probes
-------------------
For more information see function preAssociationProbeFiltering
Making sure we have at least 5% of beta values lesser than 0.3 and 5% of beta values greater 0.3.
Removing 99265 probes out of 124472
Calculating Pp (probe - gene) for all nearby genes
|                                                    |  0%                      Error in Exps[Gene, , drop = FALSE] : subscript out of bounds
Error in if (value == n) { : argument is of length zero

Everything was working smoothly until that line. And I am having a difficult time interpreting the error. If anyone could help me solve this.

My R and TCGAbiolinks is up to date.

best, Nikolay

tiagochst commented 3 years ago

@NikolayAlabi Sorry, this is an ELMER issue. I would need to debug your code/object to understand the issue. Did you check nearGenes_lowvsint, nearGenes_lowvshig, and nearGenes_intvshig objects ? It seems there is a problem when trying to subset the gene expression data. For example, the gene expression row names inside the MAE object does not match the one in the probe/gene pairs.