无意中翻到生信大神果子老师和洲更老师神仙大家的帖子，样本足够多时非参数检验就够用了（DESeq2可以先等等）。
wx搜索：《学好统计，我10分钟就完成了别人服务器通宵才完成的分析。》
下面是自己练习的代码，用的是TCGA前列腺癌的Counts数据：
rm(list = ls())
library(DESeq2)
library(tidyverse)
library(pheatmap)
library(vioplot)
rt <- data.table::fread("F:/TCGA-counts/HTSeq_Counts_PRAD.txt",data.table = F) %>%
column_to_rownames("Tags")
metadata <- data.frame(TCGA_id=colnames(rt),
sample=factor(ifelse(str_sub(colnames(rt),14,15)<10,"tumor","normal"),levels = c("normal","tumor"))) #分组信息
################################################################################
#根据文库大小标准化
expr_mt <- rt/rep(colSums(rt),each=nrow(rt))*1e6
#过滤低表达基因
expr_mt <- expr_mt[rowSums(expr_mt>0)>(ncol(expr_mt)/3),]
#wilcox.test差异分析
cancer_sample <- metadata[metadata$sample=="tumor","TCGA_id"]
normal_sample <- metadata[metadata$sample=="normal","TCGA_id"]
cancer_df <- expr_mt[,colnames(expr_mt) %in% cancer_sample]
normal_df <- expr_mt[,colnames(expr_mt) %in% normal_sample]
#计算logFC
logFC <- log2(rowMeans(as.matrix(cancer_df))/rowMeans(as.matrix(normal_df)))
library(future.apply)
plan(multisession)
p_values <- future_lapply(seq(nrow(cancer_df)),function(n){
res <- wilcox.test(x=as.numeric(cancer_df[n,]),y=as.numeric(normal_df[n,]))
res$p.value
})
p <- unlist(p_values)
p.adj <- p.adjust(p,method = "fdr")
网友评论