ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
122 stars 30 forks source link

明明就是整数啊你怎么说不是呢! #571

Closed ixxmu closed 3 years ago

ixxmu commented 3 years ago

https://mp.weixin.qq.com/s/7xVnlsKI06B66FNmw8I-VQ

github-actions[bot] commented 3 years ago

明明就是整数啊你怎么说不是呢! by 生信星球

 今天是生信星球陪你的第772天


   大神一句话,菜鸟跑半年。我不是大神,但我可以缩短你走弯路的半年~

   就像歌儿唱的那样,如果你不知道该往哪儿走,就留在这学点生信好不好~

   这里有豆豆和花花的学习历程,从新手到进阶,生信路上有你有我!

1.问题

今天学员用counts表达矩阵去做差异分析,被Deseq2报错说不是整数,如图:

load("TCGA-HNSC_gdc.Rdata")
exp[1:4,1:4]
##                    TCGA-BB-4224-01A
## ENSG00000000003.13             2237
## ENSG00000000419.11             1606
## ENSG00000000457.12             1063
## ENSG00000000460.15             1208
##                    TCGA-H7-7774-01A
## ENSG00000000003.13             2740
## ENSG00000000419.11             1691
## ENSG00000000457.12              803
## ENSG00000000460.15              317
##                    TCGA-CV-6943-01A
## ENSG00000000003.13             2686
## ENSG00000000419.11             1649
## ENSG00000000457.12              917
## ENSG00000000460.15              402
##                    TCGA-CN-5374-01A
## ENSG00000000003.13             2086
## ENSG00000000419.11             2333
## ENSG00000000457.12             1288
## ENSG00000000460.15             1105
class(exp)
## [1] "matrix" "array"
class(exp[1,1])
## [1] "numeric"

明明就是一个看起来挺正常的counts矩阵!然而:

library(DESeq2)
colData <- data.frame(row.names =colnames(exp), 
                      condition=Group)
if(!file.exists(paste0(cancer_type,"_dd.Rdata"))){
  dds <- DESeqDataSetFromMatrix(
    countData = exp,
    colData = colData,
    design = ~ condition)
  dds <- DESeq(dds)
  save(dds,file = paste0(cancer_type,"_dd.Rdata"))
}

报错信息说有一些数值它不是整数,于是学员过来求助。

2.解决方案

我给的解决方案是apply套上as.integer,矩阵变成整数型再做差异分析

exp2 = apply(exp, 2, as.integer)
library(DESeq2)
colData <- data.frame(row.names =colnames(exp2), 
                      condition=Group)
if(!file.exists(paste0(cancer_type,"_dd.Rdata"))){
  dds <- DESeqDataSetFromMatrix(
    countData = exp2,
    colData = colData,
    design = ~ condition)
  dds <- DESeq(dds)
  save(dds,file = paste0(cancer_type,"_dd.Rdata"))
}

搞定!

3.研究一下为啥子

我猜呢,是表达矩阵里面有些小数,所以检查一下数据里是不是有小数点

library(stringr)
k = str_detect(exp,"\\.");table(k)
## k
##    FALSE     TRUE 
## 13332434    61492

确实是有的,那么提取几个出来看看

x = head(exp[k]);x
## [1] 8903  947 9114 9451 9157  876

我以为见鬼了呢!这明明就是整数,哪里有小数点啊!

x
## [1] 8903  947 9114 9451 9157  876
class(x)
## [1] "numeric"

最后的破案还是来自于尝试:

as.character(x)
## [1] "8903.00000000001" "946.999999999999"
## [3] "9113.99999999999" "9451.00000000001"
## [5] "9156.99999999999" "875.999999999999"

好吧,确实是小数,只不过太接近于整数,所以被R语言显示成了整数,真是坑你没商量哦

继续刨根问底儿,得知她的表达矩阵是从xena下载的,xena 的counts数据被取了log2(count+1),这个操作是可逆的 ,2^exp + 1即可。但做deseq2差异分析还需要加上一步exp2 = apply(exp, 2, as.integer)

补充一个不太重要的小知识

R语言中的numeric是数据类型的一种,它又分integer和double,默认都是double,如果想明确指定要integer可在写代码时加L或者用as.integer转换。

class(1)
[1] "numeric"
class(1L)
[1] "integer"
typeof(1)
[1] "double"
typeof(1L)
[1] "integer"

插个小广告!

生信零基础入门学习小组

生信入门班(四周线上直播课,长期开班)

数据挖掘班(医生/医学生首选,三周线上直播课,长期开班)

一起来学单细胞吗?

生信星球答疑公告