数据是TPM，想做差异只能转换？

叉叉滴同学 叉叉滴同学的生信笔记 2023-04-21 16:00 发表于广西

目前最常用的三个差异分析方法：limma，edgeR，DEseq

公众号之前的推文也有对其进行介绍👉三包差异分析

利用这三种包进行差异分析，数据都得使用Counts；

然而，后续的许多分析使用的数据一般都是log2(TPM+1)；

后续想要补一个差异分析，只能另外收集Counts数据或者对TPM进行转换。

虽然

有上述两个简单明了的解决思路，但是新的问题不免会让同学们陷入纠结

👉Counts数据收集之后还需要对照着已经处理好的TPM数据的标准再处理一边。麻烦不说，处理时的代码报错真的是想让人砸电脑......

👉直接利用公式，将TPM数据转换为Counts。这听上去似乎可靠可行，但需要实现从公式到代码......也是让人头疼

难道

就没有一个不用另外收集数据或者另外转换的方法吗

👇👇👇

还真有，

Wilcoxon秩和检验

老样子，概念这方面需要同学们自行去了解

下面直接进行代码演示👇

①首先需要准备TPM数据这是肯定的

具体方法可查看推文👉新版TCGA数据提取

②数据前期处理并标准化

#读取TPM矩阵readCount<-read.table(file="combined_RNAseq_TPM.txt", header = T, row.names = 1, stringsAsFactors = F,check.names = F)# edger 标准化并删除低表达基因library(edgeR)library(limma)group1=sapply(strsplit(colnames(readCount),"\\-"), "[", 4)group1=sapply(strsplit(group1,""), "[", 1)group1=gsub("2", "1", group1)y <- DGEList(counts=readCount,group=group1)#删除过低表达量基因keep <- filterByExpr(y)y <- y[keep,keep.lib.sizes=FALSE]##进行TMM标准化并转移到CPM（百万计数）y <- calcNormFactors(y,method="TMM")count_norm=cpm(y)count_norm<-as.data.frame(count_norm)

③对每个基因进行Wilcoxon秩和检验

矩阵的基因动不动就成千上万，这里用多线程帮同学们提了提速

（6万个基因运算不超过10秒）

# 对每个基因进行Wilcoxon秩和检验library(future.apply)plan(multisession)pvalues <- future_lapply(1:nrow(count_norm),function(i){  data<-cbind.data.frame(gene=as.numeric(t(count_norm[i,])),group1)  p=wilcox.test(gene~group1, data)$p.value  return(p)})fdr=p.adjust(pvalues,method = "fdr")

④计算每个基因的差异倍数（logFC）

# 计算每个基因的倍数变化group2=factor(group1)conditionsLevel<-levels(group2)dataCon1=count_norm[,c(which(group1==conditionsLevel[1]))] #肿瘤dataCon2=count_norm[,c(which(group1==conditionsLevel[2]))] #正常#肿瘤比正常，logFC大于0则为“基因在肿瘤上调”；若为正常比肿瘤，则logFC大于0表示“基因在正常中上调”foldChanges=log2(rowMeans(dataCon1)/rowMeans(dataCon2))

⑤基于FDR阈值的输出结果

pvalues0=t(as.data.frame(pvalues))outRst<-data.frame(log2foldChange=foldChanges, pValues=pvalues0, FDR=fdr)rownames(outRst)=rownames(count_norm)outRst=na.omit(outRst)fdrThres=0.05

最后导出文件

write.table(outRst[with(outRst, ((log2foldChange> 1 | log2foldChange< (-1)) & FDR < 0.05 )),], file="wilcoxout.tsv",sep="\t", quote=F,row.names = T,col.names = T)

最后得到差异结果文件

图片

是不是非常亲切图片

今天就跟同学们分享到这，需要代码及示例文件的同学可在下方扫码获取

👇感谢关注与支持，欢迎进群交流👇

图片

喜欢此内容的人还喜欢

微信扫一扫
关注该公众号