Блог

Я пытаюсь провести анализ дифференциальных выражений в R (используя DESeq), и все работает нормально, когда я запускаю фрагменты вручную.

Однако при вязании я получаю следующую ошибку для определенного куска:

 Quitting from lines 346-445 (Formatted_M1M2_miRNA_NP_AJ-copy.rmd)  Error in h(simpleError(msg, call)) :   error in evaluating the argument 'i' in selecting a method for function '[': object 'idx' not found Calls: lt;Anonymousgt; ... withVisible -gt; eval -gt; eval -gt; [ -gt; .handleSimpleError -gt; h Execution halted  

Это тот самый кусок:

 ddsM1 = DESeqDataSetFromMatrix(countData = countsM1, colData = phenoM1_2, ~ pefull)  ddsM1 lt;- estimateSizeFactors(ddsM1) norm.ctsM1 lt;- counts(ddsM1, normalized=TRUE)  ddsM2 = DESeqDataSetFromMatrix(countData = countsM2, colData = phenoM2_2, ~ pefull)  ddsM2 lt;- estimateSizeFactors(ddsM2) norm.ctsM2 lt;- counts(ddsM1, normalized=TRUE)  mm_1 lt;- model.matrix(~ pefull, colData(ddsM1)) mm0_1 lt;- model.matrix(~ 1, colData(ddsM1)) norm.ctsM1 lt;- norm.ctsM1[rowSums(norm.ctsM1) gt; 0,] fitM1 lt;- svaseq(norm.ctsM1, mod=mm_1, mod0=mm0_1, n.sv=2) plot(fitM1$sv,pch=19,col="blue")  mm_2 lt;- model.matrix(~ pefull, colData(ddsM2)) mm0_2 lt;- model.matrix(~ 1, colData(ddsM2)) norm.ctsM2 lt;- norm.ctsM2[rowSums(norm.ctsM2) gt; 0,] fitM2 lt;- svaseq(norm.ctsM2, mod=mm_2, mod0=mm0_2, n.sv=2) plot(fitM2$sv,pch=19,col="blue")    ddsM1$SV1 lt;- fitM1$sv[,1] ddsM1$SV2 lt;- fitM1$sv[,2]  phenoM1_2$SV1 lt;- fitM1$sv[,1] phenoM1_2$SV2 lt;- fitM1$sv[,2]   ddsM2$SV1 lt;- fitM2$sv[,1] ddsM2$SV2 lt;- fitM2$sv[,2]  phenoM2_2$SV1 lt;- fitM2$sv[,1] phenoM2_2$SV2 lt;- fitM2$sv[,2]  ###################  design= ~ pefull   SV1   SV2 ddsM1_SV = DESeqDataSetFromMatrix(countData = countsM1, colData = phenoM1_2, design = design)  ddsM2_SV = DESeqDataSetFromMatrix(countData = countsM2, colData = phenoM2_2, design = design)  ddsM1_SV lt;- ddsM1_SV[idx,] ddsM1_2=DESeq(ddsM1_SV) resultsNames(ddsM1_2)  ddsM2_SV lt;- ddsM2_SV[idx,] ddsM2_2=DESeq(ddsM2_SV) resultsNames(ddsM2_2)   ddsM2_SV lt;- ddsM2[idx,] ddsM2_2=DESeq(ddsM2_SV) resultsNames(ddsM2_2)  basemean=rowMeans(counts(ddsM1_2, normalized=TRUE)) basemean=rowMeans(counts(ddsM2_2, normalized=TRUE))   res0M1 lt;- results(ddsM1_2, alpha=0.05) mcols(res0M1, use.names=TRUE) summary(res0M1) tabrM1=results(ddsM1_2, cooksCutoff=FALSE)   res0M2 lt;- results(ddsM2_2, alpha=0.05) mcols(res0M2, use.names=TRUE) summary(res0M2) tabrM2=results(ddsM2_2, cooksCutoff=FALSE)  plotMA(res0M1, ylim=c(-7,7)) #Red points have adjusted p value lt; 0.1 plotMA(res0M2, ylim=c(-7,7))  #Estimating batches by svaseq followed by providing covariates in the DESeq design is expected to perform either equally good or better than adding batch as a factor in glm.. So design matrix for deseq could be something like:  #model.matrix(as.formula(paste0("~ genotype",fit))) #Why are you doing this?  resShrinkM1 = lfcShrink(ddsM1_2, coef=2) plotMA(resShrinkM1, ylim=c(-5,5))  resShrinkM2 = lfcShrink(ddsM2_2, coef=2) plotMA(resShrinkM2, ylim=c(-5,5))  #filtering results resSigM1 lt;- subset(res0M1, pvalue lt; 0.05) resSigM1   resSigM2 lt;- subset(res0M2, pvalue lt; 0.05) resSigM2  rownames(resSigM1) rownames(resSigM2)  head(res0M1[order(res0M1$padj),], 5) # order by padjusted and print the top 5 head(res0M2[order(res0M2$padj),], 5)   resLFC1 = results(ddsM1_2, lfcThreshold = 1) #genes with a fold change of 2 (or -2) - i.e.doubling/halving. summary(resLFC1)  resLFC2 = results(ddsM2_2, lfcThreshold = 1) #genes with a fold change of 2 (or -2) - i.e.doubling/halving. summary(resLFC2)   

Что я делаю не так?