查看原文
其他

我是不会运行你的代码吗?不,我是不会导入自己的数据!

生信宝典 生信宝典 2023-06-29

常常遇到有人问起看到分享的教程导入数据的方式是data(dune)等直接调用系统的数据,而自己怎么读入自己的数据呢?

对于初学者来讲,这确实是个问题。如何准备数据、拿到正确格式的数据并导入后续的代码进行分析,是学习和应用过程中的第一个拦路虎。

为什么教程会习惯使用内置数据?

  1. 简单省事、便携可重复;这是内置数据的优势之一;

  2. 内置数据模式清晰,通常可以获得较好的结果;这是内置数据的优势之二;

  3. 别人用这个,我也用这个,这是一个偷懒的做法。

  4. 每个人常识不同。作者可能觉得这个太简单而忽略了初学者的需求。(生信学习学的是什么?常识!

但内置数据的频繁使用是导致初学者学习这个教程时经常提出上面这个问题的原因。

我不太赞成教程里面用使用内置数据,原因是:

  1. 对不会读入数据的人不友好;

  2. 不利于探索这篇教程用于实际数据时可能会遇到的问题。示例数据无脑运行,自己的数据无显著差异。

如果要使用内置数据,也需要额外提供一些信息:

  1. 详细描述内置数据的格式和生物含义,及与真实数据的对应,可以参考画一个带统计检验的PCoA分析结果

  2. 提供真实数据的格式示例和读入真实数据的代码,弥补这个“鸿沟”;
    比如写这篇文章:你的adonis用对了吗?不同因素的顺序竟然对结果有很大影响就是因为示例数据有显著差异,而自己的数据无差异。所以才从原理上其理解计算过程,并探寻解决方案。

  3. 提及可能出现的问题的解决;这也是操作了多套实际数据后,才能写出的部分。

那假如教程没有提供这么详细,自己又得用这个教程,怎么做呢?

自己如何根据教程的数据准备并读入自己的数据

1. 查看数据的结构,了解数据的构成

既然教程提供了测试数据集,不妨仔细看看测试数据集的特征,没准就找着规律了。

我们以前面文章提到的dune数据集为例,查看下其结构特征。

行名字是数字,列名字是字符串(如果我们对这些字符串不熟悉,对我们来说就没任何意义;每个字符都认识,串一起就不知道是啥了~~),中间的值是整数。除此外也看不出其它信息了。

library(vegan)
data(dune)
head(dune)

## Achimill Agrostol Airaprae Alopgeni Anthodor Bellpere Bromhord Chenalbu Cirsarve Comapalu Eleopalu Elymrepe Empenigr
## 1 1 0 0 0 0 0 0 0 0 0 0 4 0
## 2 3 0 0 2 0 3 4 0 0 0 0 4 0
## 3 0 4 0 7 0 2 0 0 0 0 0 4 0
## 4 0 8 0 2 0 2 3 0 2 0 0 4 0
## 5 2 0 0 0 4 2 2 0 0 0 0 4 0
## 6 2 0 0 0 3 0 0 0 0 0 0 0 0
## Hyporadi Juncarti Juncbufo Lolipere Planlanc Poaprat Poatriv Ranuflam Rumeacet Sagiproc Salirepe Scorautu Trifprat Trifrepe
## 1 0 0 0 7 0 4 2 0 0 0 0 0 0 0
## 2 0 0 0 5 0 4 7 0 0 0 0 5 0 5
## 3 0 0 0 6 0 5 6 0 0 0 0 2 0 2
## 4 0 0 0 5 0 4 5 0 0 5 0 2 0 1
## 5 0 0 0 2 5 2 6 0 5 0 0 3 2 2
## 6 0 0 0 6 5 3 4 0 6 0 0 3 5 5
## Vicilath Bracruta Callcusp
## 1 0 0 0
## 2 0 0 0
## 3 0 2 0
## 4 0 2 0
## 5 0 2 0
## 6 0 6 0



2. 查看数据的帮助

从数据结构和行列名字上得不到有用信息,那我们查看下帮助信息。

?dune

dune is a data frame of observations of 30 species at 20 sites. The
species names are abbreviated to 4+4 letters (see make.cepnames).

这告诉我们什么呢?这套数据包含了30个物种在20个样品的丰度信息。从dim(dune)可以看出这是一个20行X30列的矩阵;可以推测出,每一行是一个样品,每一列是一个物种
(另一个佐证是列名字长度确实为8个字符,与物种名字的4+4缩写一致)。

注:如果对数据还有疑虑,建议谷歌下数据。常见内置数据集都会有文章描述其信息,可用于佐证你的判断。

dim(dune)

## [1] 20 30

这个格式跟我们通常的OTU丰度表
(我们的表通常是每一行是一个物种,每一列是一个样品)略有不同。

3. 基本判断后,读入我们的数据,做可能的转换

如果我们有一个OTU丰度表,怎么读入并转成这个格式呢?

text <- "ID\tSamp1\tSamp2\tSamp3\tSamp4
OTU1\t2\t13\t14\t15
OTU2\t12\t13\t8\t10
OTU3\t22\t10\t14\t11"

otu_table <- read.table(text=text, sep="\t", row.names=1, header=T)

读入OTU丰度表,第一行为列名字,第一列为行名字。

otu_table <- read.table("otutable_rare",sep="\t", row.names=1, header=T)

根据上面的分析做一个转置,就可以获得可用于后续分析的输入数据了。

otu_table_t <- as.data.frame(t(otu_table))
otu_table_t

## OTU1 OTU2 OTU3
## Samp1 2 12 22
## Samp2 13 13 10
## Samp3 14 8 14
## Samp4 15 10 11


4. 示例数据中的整数代表什么意思?

这个是比较难确定的部分,只有两个判断方法:1)
教程中作者能够提及(这是最准确的方法);2)凭经验猜测。

这里涉及到另外一个经常会被问起的问题:

我这一步操作需要提供原始数据,还是标准化之后的数据?

绝大多数情况下,我们需要提供的都是标准化之后的在不同样品之间可比的数据。因为:1)我们的需求是比较不同样品的差异,数据需要在样品间可比;2)绝大部分工具是不会对数据做标准化处理的,要么直接用,要么做一些不影响数值关系的转换;3)如果某个工具自己内部会对数据做标准化,它一定会在帮助中提及,常见的比如DESeq2edgeRlimma,除了这两个半(limma算半个,因为它也可以接收标准化后的数据),一时想不起还有哪些工具是接受原始数据的。单细胞的Seurat包算是个例外,它内部调用了一些标准化算法,可以通过参数关掉。

5. 查看更多教程,总会遇到有详细描述所需数据结构的教程。

6. 跟着感觉走,不管三七二十一读进来试试,出现异常或报错再调整。学程序不是做实验,试错成本没有那么大,光看不练是假把式,大胆试才是王道。

7. 最后一步,跟教程作者沟通。我们的教程问题,欢迎在http://www.ehbio.com/Esx发帖讨论;自己努力后,带着问题和思路的讨论更容易获得解答。

生物教程还是得使用生物数据!!!

往期精品(点击图片直达文字对应教程)

机器学习

后台回复“生信宝典福利第一波”或点击阅读原文获取教程合集



您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存