排序

两个月前的一个学徒作业:绘图本身很简单但是获取数据很难,完成率超级低,仅仅接到了不到十个邮件,而且有3个人做的是错的!!超级尴尬,其中有一个错误很明显,就是自以为是的排序,然后比对肿瘤组织和配对的正常组织的表达量,其实呢,排序错误会导致配对失败。

下面让优秀学徒具体讲解一下:

rm(list = ls())
options(stringsAsFactors = F)
library(ggstatsplot)
load('dat_input22.Rdata')

目前我的 dat 数据是这样的,可以看到同一个病人是有肿瘤组织和配对的正常组织的表达量的,而且呢,理论上是每一行一个样品的表达量信息:

image-20201103141112594

对 pid 这一列排序后,group 这一列应该是相对应的奇数时是肿瘤组,偶数正常组。这时候就出现了问题,排列的没有规律性,如下:

image-20201103141216110

后面的数据就无法取,于是思考了一下两列的排序问题。(起初我并没有想到这一点,而是采用了其它复杂的方法完成了这个目标。但是jimmy老师点醒了我:凡是Excel能实现的数据操作,理论上R语言也可以,其实就是按照两列元素进行排序)

本来就只是一个简单的排序问题,随便搜搜就会有很好的答案,例如这样

df = dat
df = df[order(df[,4],df[,3],decreasing=TRUE),]

也就是说上面的代码呢,首先按照第4列排完序了,然后再来排一下第3列,我的数据也就得到了解决。

image-20201103141622979

排列的整整齐齐:

image-20201103141848605

并且后续的分析只需要在正常组和原位肿瘤组织中,不需要转移的肿瘤的这两个数据,应该删掉就行:

image-20201103142623820

删除了多余的转移肿瘤的数据之后其实就完美了(都是那多出来的四个数据的问题,不然,第一次按照一列排序就可以很好)

之后就可以分别取出肿瘤样本和正常样本对应的 TP53 的表达量:

d=cbind(d[seq(1,nrow(d),by=2),],d[seq(2,nrow(d),by=2),])
identical(dat[,4],dat[,8])
head(dat)
normal=dat[,6];tumor=dat[,2]

这时候我的数据就结束了。

但是,

人们有时候真是有些奇奇怪怪的要求,索性研究一下排序问题。

哈哈,其实没有了,排序哪有什么问题,没有遇到具体问题时,哪会有排序的需求。下面还是让jimmy老师给大家讲解排序的需求吧:

扩展 筛选基因

我们读取一个表达矩阵文件,这个 GSE75688_GEO_processed_Breast_Cancer_raw_TPM_matrix.txt.gz 你应该是知道如何下载的

rm(list=ls())
options(stringsAsFactors = F)
a=read.table('GSE75688_GEO_processed_Breast_Cancer_raw_TPM_matrix.txt.gz',
 header = T)
a[1:4,1:4]
tail(a[,1:4])

表达矩阵如下:

> a[1:4,1:4]
 gene_id gene_name gene_type BC01_Pooled
1 ENSG00000000003.10 TSPAN6 protein_coding 2.33
2 ENSG00000000005.5 TNMD protein_coding 0.00
3 ENSG00000000419.8 DPM1 protein_coding 60.70
4 ENSG00000000457.9 SCYL3 protein_coding 47.93
> tail(a[,1:4])
 gene_id gene_name gene_type BC01_Pooled
57910 ERCC-00168 ERCC-00168 ERCC 0.00
57911 ERCC-00170 ERCC-00170 ERCC 0.00
57912 ERCC-00171 ERCC-00171 ERCC 0.00
57913 SPIKE1 EC2 SPIKE_IN 14940.70
57914 SPIKE2 EC15 SPIKE_IN 985.82
57915 SPIKE3 EC18 SPIKE_IN 0.00

可以看到前面的3列是基因信息,后面的才是表达矩阵的各个样品表达量信息。

ids=a[,1:3]
head(ids)
a=a[,4:ncol(a)]
a[1:4,1:4]
kp=rowSums(a>1) > 20
table(kp)
a=a[kp,]
ids=ids[kp,]
rownames(a)=ids[,1]
a[1:4,1:4]

虽然说第一列是ensembl的基因ID,第二列是我们想要的基因symbol。如果这个时候,我们强行把 rownames(a)=ids[,2] 就会报错,如下:

image-20201105083313284

可以看到有大量的基因出现了多次,因为它们其实对应着不同的ensembl的基因ID,但是我们最后仍然是想要基因symbol。这个时候,我们就可以应用起来了我们的两列排序技巧:

image-20201105083650506

可以看到, 我们的ids数据框,首先是按照基因的symbol排序了,然后按照基因表达量排序了,所以可以简单的去冗余就拿到了合适的基因。全部的代码如下:

ids$median=apply(a,1,median)
ids=ids[order(ids$gene_name,ids$median,decreasing = T),]
ids=ids[!duplicated(ids$gene_name),]
a=a[ids$gene_id,]
rownames(a)=ids[,2]

这样,你的表达量矩阵,就是以唯一的基因symbol命名的啦。

虽然说两个不同的ensembl的基因ID,对应着同样的基因symbol,但是我们的挑选策略是,仅仅是保留表达量大的那个ensembl的基因ID。

如果你要问为什么两个不同的ensembl的基因ID会对应着同样的基因symbol

只能说是因为id本来就不统一,而且基因数量那么多,是超出人类认知范围的!这些知识点统称为生物信息学背景知识咯,甚至可以写一本书:

为什么要转换id?
有多少种ID?
什么id权威?
id是一一对应的吗?
ID是什么生信组织维护?
id有版本吗?
id一定正确吗?
什么情况下选择什么id?
不同数据库下载的id对应表一定一样吗?

Comments are closed.