count转化为TPM

library(biomaRt)
#查看基因组参数
mart = useMart('ensembl')
listDatasets(mart)

#你需要哪个基因组,就复制它在dataset列里的词,放在下面这行的`dataset = `参数里
#此处以人类为例,植物参考注一
bmart <- biomaRt::useMart(biomart = "ENSEMBL_MART_ENSEMBL", 
                          dataset = "hsapiens_gene_ensembl",
                          host = "www.ensembl.org")

# 从输入数据里提取基因名
feature_ids <- rownames(expMatrix)##就是你需要ENSEMBL的例如ENSG00000000003


attributes = c(
  "ensembl_gene_id",
  #"hgnc_symbol",
  "chromosome_name",
  "start_position",
  "end_position"
)
filters = "ensembl_gene_id"

feature_info <- biomaRt::getBM(attributes = attributes, 
                               filters = filters, 
                               values = feature_ids, mart = bmart)
mm <- match(feature_ids, feature_info[[filters]])
feature_info_full <- feature_info[mm, ]
rownames(feature_info_full) <- feature_ids

# 计算基因的有效长度
eff_length <- abs(feature_info_full$end_position - feature_info_full$start_position)
feature_info_full <- cbind(feature_info_full,eff_length)

eff_length2 <- feature_info_full[,c(1,5)]

x <- expMatrix / eff_length2$eff_length

tpm = t( t(x) / colSums(x) ) * 1e6 

### 将 GEO 数据库中的 Count 数据转换TPM #### 工作流程概述 为了将 GEO 数据库中的 Count 数据转换TPM (Transcripts Per Million),需要完成以下几个核心操作:获取基因长度信息、合并计数数据与基因长度矩阵以及应用标准化公式。 --- #### 1. 准备阶段 在开始之前,需设置工作环境并加载必要的 R 包。以下是常用工具包及其功能说明: - `edgeR` 或 `DESeq2`: 处理 RNA-seq 计数数据。 - `tximport`: 提供从文件导入转录本级别的表达量的功能。 - `dplyr`: 方便处理和整理数据框。 代码如下: ```r library(edgeR) library(tximport) library(dplyr) ``` --- #### 2. 获取 Counts 数据 假设已经下载了 GEO 数据库中的 counts 文件(通常是一个表格),可以使用以下命令将其读入内存中: ```r counts <- read.table("GSEXXXXXX_counts.txt", header=TRUE, row.names=1, sep="\t") ``` 此步骤会创建一个名为 `counts` 的对象,其中每一行为基因 ID,每列为样本对应的 reads count [^3]。 --- #### 3. 下载基因长度信息 TPM 的计算依赖于基因的长度信息。可以通过外部资源(如 GENCODE 或 RefSeq)或者直接利用 GEO 中附带的相关元数据来提取这些。如果已知基因长度存储在一个单独的 CSV/TXT 文件里,则可通过下面的方式加载它: ```r gene_lengths <- read.table("gene_length_info.txt", header=TRUE, row.names=1, sep="\t") ``` 确保该表的第一列包含唯一的基因标识符,并且第二列表达的是对应基因的实际碱基对数目。 --- #### 4. 合并 Counts 和 Gene Lengths 为了让后续计算更加简便,应先将上述两部分结合起来形成一个新的数据结构——即扩展后的矩阵形式,在这里我们采用 dplyr 来实现这一目标: ```r merged_data <- inner_join(as.data.frame(t(counts)), gene_lengths, by="row.names") rownames(merged_data) <- merged_data$Row.names merged_data$Row.names <- NULL ``` 此时得到的新变量 `merged_data` 应当具有这样的布局:行代表各个基因;除了最后的一列记录着它们各自的物理尺寸外,其余各列则分别表示不同样品下的原始测序覆盖度统计结果。 --- #### 5. 定义 TPM 转换函数 基于定义好的数学模型,编写一段简单的脚本来执行批量转化过程: ```r convert_to_tpm <- function(count_matrix){ effective_count_per_kb <- sweep(count_matrix[, -ncol(count_matrix)], MARGIN=1, FUN="/", STATS=count_matrix[, ncol(count_matrix)] / 1e3) tpm_values <- sweep(effective_count_per_kb , MARGIN=2 , FUN="/",STATS=(apply(effective_count_per_kb ,2,sum)/1e6)) return(tpm_values) } tpms <- convert_to_tpm(merged_data) ``` 这段逻辑首先调整每个基因单位千碱基数目的相对丰度比例关系,接着再除以整个实验体系内的总有效映射片段数量从而获得最终标准化后的 TPM 数组。 --- #### 6. 输出结果至文件 一旦完成了全部数变换之后,就可以把所得成果导出到本地磁盘上以便进一步分析或分享给其他研究者们查看啦! ```r write.csv(as.data.frame(t(tpms)), file="GSEXXXXXX_TPM_expression.csv", quote=FALSE) ``` 这样就成功实现了从 raw counts 到 TPM 表达谱之间的转变全过程[^2]。 --- ###
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值