rnaseq_对表达矩阵进行标准化

hudi18

已于 2024-04-03 15:42:30 修改

阅读量934

点赞数 4

CC 4.0 BY-SA版权

分类专栏： rnaseq 文章标签： r语言

于 2024-04-02 21:34:25 首次发布

本文链接：https://blog.csdn.net/qq_46454013/article/details/137291227

rnaseq 专栏收录该内容

6 篇文章

订阅专栏

本文详细描述了作者在转录组分析过程中对表达矩阵进行标准化的过程，包括将count值转换为TPM（TranscriptsPerKilobaseMillion）以及使用DESeq2的VST（方差稳定变换）进行标准化。这些步骤对于整理和理解复杂的数据至关重要。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

从2022年十月份到现在，做了那么久的转录组分析，从来没有记录一点东西，导致脑子里面的东西越来越乱，所以我把我所做的工作记录一下，供自己以及同行参考。

上游分析做的比较久远，先把最近做的下游分析写一下。

表达矩阵的标准化

由上游分析得到了表达矩阵。

在后续分析中，要先把count值进行标准化，我是转换成tpm和用Deseq2进行VST标准化

count转换成tpm

###### counts转tpm  ################
counts <- rs[,7:ncol(rs)]
rownames(counts) <- rs$transcript_id
geneid_efflen <- subset(len2,select = c("X.ID","transcript_length"))
colnames(geneid_efflen) <- c("geneid","efflen")
geneid_efflen_fc <- geneid_efflen #用于之后比较
dim(geneid_efflen)
efflen <- geneid_efflen[match(rownames(counts),
                              geneid_efflen$geneid),
                        "efflen"]
TPM (Transcripts Per Kilobase Million)  每千个碱基的转录每百万映射读取的Transcripts
counts2TPM <- function(count=count, efflength=efflen){
  RPK <- count/(efflength/1000)       #每千碱基reads (“per million” scaling factor) 长度标准化
  PMSC_rpk <- sum(RPK)/1e6        #RPK的每百万缩放因子 (“per million” scaling factor ) 深度标准化
  RPK/PMSC_rpk                    
}  
tpm <- as.data.frame(apply(counts,2,counts2TPM))
rs <- tpm
q <- rs
transcript_id <- rownames(q)
rs <- cbind(transcript_id,q)
rownames(rs) <- NULL

VST标准化

########### VST标准化 ############
#导入counts数据矩阵
count <- read.csv("C:\\Users\\DI\\Desktop\\PRJNA516151_lnc\\lncRNA_Gastrocnemius_transcript_matrix.csv",header=T,row.names=1)
## 过滤在所有重复样本中小于1的基因
count <- count[rowMeans(count)>1,]
##载入样本信息
data <- read.table("C:\\Users\\DI\\Desktop\\PRJNA516151_lnc\\sample_data.txt",header = T,row.names = 1)
#变为因子数据
data[,1] <- as.factor(data$Type)
all(rownames(data) %in% colnames(count))
all(rownames(data) == colnames(count))
dds <-  DESeqDataSetFromMatrix(countData = count,colData = data,design = ~ Type)
dim(dds)
#过滤
dds <- dds[rowSums(counts(dds)) > 1,]
nrow(dds) 
#方差稳定变换
vsd <- vst(object=dds,blind=T) 
head(assay(vsd),3)
rs <- assay(vsd)
rs <- as.data.frame(rs)
q <- rs
transcript_id <- rownames(q)
rs = as.data.frame(lapply(rs,as.numeric))
rs <- cbind(transcript_id,q)
rownames(rs) <- NULL

到这里，标准化就基本完成，我后续分析用的是VST标准化。