最讨厌这样的样品命名体系

生信技能树

发布于 2025-06-09 21:13:39

1600

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

最近看到了一个转录组-时间序列的实验设计，数据分析环节开创了TrendCatcher这个分析方法，就安排我们马拉松授课的学员们试试看，复现一下。

结果学员在第一步就卡住了，但是并不是学员本身编程能力的问题，是数据集自己的bug。在GEO界面就可以下载到表达量矩阵文件，然后简单的读取即可：

data<-data.table::fread("GSE211266_raw_count.csv.gz",
                        data.table = F)

data=data[!duplicated(data$V1),]
mat<-data[,c( :ncol(data))]
rownames(mat)=data[,]
mat[:,:]
keep_feature <- rowSums (mat > ) >  ;table(keep_feature)
ensembl_matrix <- mat[keep_feature, ]  
rownames(ensembl_matrix)=rownames(mat)[keep_feature]
ensembl_matrix[:,:]

可以看到，是一个很正常的转录组测序后的表达量count矩阵：

                   JR-10S-RK-01 JR-10S-RK-010 JR-10S-RK-02 JR-10S-RK-03
ENSMUSG00000051951            0            18            0            9
ENSMUSG00000102851            0             0            0            2
ENSMUSG00000104017            0             0            0            0
ENSMUSG00000103025            0             0            0            0

可以看到，这个时候是基因id是ens开头的，不符合人类的阅读习惯啦。这个时候可以通过人工智能大模型简单的学习一下Ensembl ID的命名规则和转换方法，Ensembl ID是Ensembl数据库使用的基因、转录本、外显子、蛋白等分子特征的唯一标识符。其命名规则如下：

格式：ENS[物种符号][分子特征][独一无二的11位数字]
- G：基因（Gene）
- T：转录本（Transcript）
- E：外显子（Exon）
- P：蛋白（Protein）
- R：调控特征（Regulatory feature）
- FM：Ensembl蛋白家族（Ensembl protein family）
- GT：基因树（Gene tree）
- ENS：Homo sapiens（人类）
- ENSMUS：Mus musculus（小鼠）
- ENSRNO：Rattus norvegicus（大鼠）
- 物种符号：
- 分子特征：
版本号：ID末尾的小数点和数字表示版本号，例如ENSG00000223972.5。版本号在特征发生变化时会增加。
这里推荐使用AnnoProbe进行id转化，代码如下所示：

library(AnnoProbe)
head(rownames(ensembl_matrix))
ids=annoGene(rownames(ensembl_matrix),'ENSEMBL','mouse')
head(ids)
tail(sort(table(ids$biotypes)))
ids=ids[ids$biotypes=='protein_coding',]
ids=ids[!duplicated(ids$SYMBOL),]
ids=ids[!duplicated(ids$ENSEMBL),]
symbol_matrix= ensembl_matrix[match(ids$ENSEMBL,rownames(ensembl_matrix)),]
rownames(symbol_matrix) = ids$SYMBOL 
symbol_matrix[:,:]

然后大家就能看到了样品名字问题：

> symbol_matrix[1:4,1:4] 
       JR-10S-RK-01 JR-10S-RK-010 JR-10S-RK-02 JR-10S-RK-03
Xkr4              0            18            0            9
Rp1               1             1            0            1
Sox17          6687          9084         9311        11551
Mrpl15          365           647          417          475
> colnames(symbol_matrix)
 [1] "JR-10S-RK-01"  "JR-10S-RK-010" "JR-10S-RK-02"  "JR-10S-RK-03" 
 [5] "JR-10S-RK-04"  "JR-10S-RK-05"  "JR-10S-RK-06"  "JR-10S-RK-07" 
 [9] "JR-10S-RK-08"  "JR-10S-RK-09"  "JR-16S-RK-11"  "JR-16S-RK-16" 
[13] "JR-16S-RK-17"  "JR-16S-RK-18"  "JR-16S-RK-19"  "JR-16S-RK-22" 
[17] "JR-16S-RK-23"  "JR-16S-RK-24"  "JR-16S-RK-25"

因为上面的表达量矩阵里面的样品名字没有办法跟GEO数据库里面的记录对应，因为上面的名字并不是普通的阿拉伯数值的排序，反而是字母表的排序。这种情况下，排序的结果会与我们通常期望的数值排序不同 ：

1. 字母表排序（Lexicographical Order）

当数字被当作字符串处理时，排序会按照字母表的顺序进行。例如：

数字 10、2 和 1 会按照字母表顺序排序为：1、10、2。
数字 23、4 和 5 会按照字母表顺序排序为：23、4、5。

2. 数值排序（Numerical Order）

当数字被当作数值处理时，排序会按照数值大小进行。例如：

数字 10、2 和 1 会按照数值大小排序为：1、2、10。
数字 23、4 和 5 会按照数值大小排序为：4、5、23。

  library(AnnoProbe)
  library(GEOquery)
  gset = getGEO("GSE211266", destdir = '.', getGPL = F) 
  pd = pData(gset[[]]) 
  colnames(pd)

这样的话，就完全没有机会保证准确无误的对应起来了：

>   pd[,c(1,40,41)]
                  title time:ch1 treatment:ch1
GSM6460125   mEC-C-YMK1        0           PBS
GSM6460126   mEC-C-YMK2        0           PBS
GSM6460127   mEC-C-YMK3        0           PBS
GSM6460128   mEC-C-YMK4        0           PBS
GSM6460129   mEC-C-YMK5        0           PBS
GSM6460130 mEC-TW1-YMK1    week1         Tumor
GSM6460131 mEC-TW1-YMK2    week1         Tumor
GSM6460132 mEC-TW1-YMK3    week1         Tumor
GSM6460133 mEC-TW1-YMK4    week1         Tumor
GSM6460134 mEC-TW1-YMK5    week1         Tumor
GSM6460135 mEC-TW1-YMK6    week1         Tumor
GSM6460136 mEC-TW2-YMK3    week2         Tumor
GSM6460137 mEC-TW2-YMK4    week2         Tumor
GSM6460138 mEC-TW2-YMK5    week2         Tumor
GSM6460139 mEC-TW2-YMK6    week2         Tumor
GSM6460140 mEC-TW3-YMK3    week3         Tumor
GSM6460141 mEC-TW3-YMK4    week3         Tumor
GSM6460142 mEC-TW3-YMK5    week3         Tumor
GSM6460143 mEC-TW3-YMK6    week3         Tumor

其实有一个隐藏的办法是认真的阅读文章里面的生物学背景，比如提到的差异基因列表：

cg = 'Hey1, CCn3, CCn4, Smad6, Acta1, Myoz1,Myoz2'
cg = trimws(strsplit(cg,',')[[]])
cg
dat=log(edgeR::cpm(symbol_matrix)+)
cg=cg[cg %in% rownames(dat)]
n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化
n[n>]= 
n[n< -]= -
n[:,:]
pheatmap::pheatmap(t(n) )

可视化后其实是可以看到基本上的层次聚类还是很清晰的：

层次聚类还是很清晰的热图

但是，实际上也很难对应，文章里面的热图如下所示：

文章里面的热图

这种情况，就只能说是去下载这个项目的fastq文件，然后自己走转录组测序的定量流程啦

在R中，可以使用多种方法生成类似01, 02, 03, ...这样的序列，其中每个数字都以两位数的格式表示。以下是几种常见的方法：

方法1：使用`sprintf()`函数

sprintf()函数可以格式化字符串，非常适合生成这种两位数的序列。

# 生成1到10的序列
numbers <- :

# 使用sprintf格式化为两位数
formatted_numbers <- sprintf("%02d", numbers)

# 查看结果
print(formatted_numbers)

方法2：使用`paste0()`函数

paste0()函数可以将字符串拼接在一起，结合ifelse()函数可以实现类似的效果。

# 生成1到10的序列
numbers <- :

# 使用paste0格式化为两位数
formatted_numbers <- ifelse(numbers < , paste0("0", numbers), numbers)

# 查看结果
print(formatted_numbers)

方法3：使用`formatC()`函数

formatC()函数也可以用于格式化数字，生成两位数的序列。

# 生成1到10的序列
numbers <- :

# 使用formatC格式化为两位数
formatted_numbers <- formatC(numbers, width = , format = "d", flag = "0")

# 查看结果
print(formatted_numbers)

方法4：使用`str_pad()`函数（来自`stringr`包）

stringr包提供了str_pad()函数，可以方便地对字符串进行填充。

# 安装并加载stringr包
if (!require(stringr)) {
  install.packages("stringr")
}
library(stringr)

# 生成1到10的序列
numbers <- :

# 使用str_pad格式化为两位数
formatted_numbers <- str_pad(numbers, width = , pad = "0")

# 查看结果
print(formatted_numbers)

方法5：使用`paste()`函数

paste()函数也可以用于拼接字符串，结合ifelse()函数可以实现类似的效果。

# 生成1到10的序列
numbers <- :

# 使用paste格式化为两位数
formatted_numbers <- ifelse(numbers < , paste("0", numbers, sep = ""), numbers)

# 查看结果
print(formatted_numbers)