Блог

Я пытаюсь запустить анализ GO в R (я никогда не проводил этот анализ, поэтому я пробую разные пакеты), и я изо всех сил пытаюсь найти проблему с моим кодом в пакете goseq.

Я начинаю с этого кода, который создает список дифференциально экспрессируемых названий генов:

  de.genes <- rownames(res)[ which(res$padj < fdr.threshold amp; !is.na(res$padj)) ]
 

Затем я пытаюсь запустить этот код (основанный на странице 7 виньетки (https://bioconductor.org/packages/devel/bioc/vignettes/goseq/inst/doc/goseq.pdf )

  pwf <- nullp(de.genes, "hg38","geneSymbol")
 

но я получаю следующую ошибку:

  Can't find hg38/geneSymbol length data in genLenDataBase...
 Found the annotation package, TxDb.Hsapiens.UCSC.hg38.knownGene
 Trying to get the gene lengths from it.
 Error in if (matched_frac == 0) { : missing value where TRUE/FALSE needed
 In addition: Warning message:
 In grep(txdbPattern, installedPackages):argument 'pattern' has length > 1 and only the first element will be used
 

Я нашел этот форум: https://support.bioconductor.org/p/38580 / это говорит о том, что мне нужна «индикаторная переменная», но я не знаю, что это такое.

Мы будем очень признательны за любую помощь в устранении этой ошибки, или если вы знаете какие-либо другие пакеты GO, которые легко освоить. Спасибо!

Вы можете проверить поддерживаемые базы данных, hg38 не входит в их число:

 library(org.Hs.eg.db)
library(goseq)

supported[grep("hg38|hg19",supported$Genome),]
   Genome         Id  Id Description Lengths in geneLeneDataBase
4    hg19  knownGene  Entrez Gene ID                        TRUE
36   hg19    ensGene Ensembl gene ID                        TRUE
81   hg19 geneSymbol     Gene Symbol                        TRUE
98   hg38                                                  FALSE
   GO Annotation Available
4                     TRUE
36                    TRUE
81                    TRUE
98                    TRUE
 

Вы можете получить приблизительное представление о том, как это выглядит, используя hg19, у вас будут некоторые недостающие или несоответствующие значения should be ok. У вас должен быть двоичный вектор, и он должен быть назван, например:

 set.seed(111)
allgenes = keys(org.Hs.eg.db,keytype="SYMBOL")
de.genes = rbinom(100,1,0.3)
names(de.genes) = sample(allgenes,100)
 

Это выглядит так:

   GALNT5        TPRKB         CD48       OR52R1 LOC105372708 LOC112163649 
       0            1            0            0            0            0 
 

LOC105369203 LOC110121115 LOC105377654 LOC105371502 LOC101929964 HPC14
0 0 0 0 0 0
IGHD4-17 LOC101927993 HINT1 BCC3 RPL18P3 LOC108281192
0 0 0 0 0 1
RNU6-793P ИЮНЬ
0 0

Все будет хорошо:

 res = nullp(de.genes,"hg19","geneSymbol")