TCGA

TCGA差异分析|1.非参数检验

2022-02-11  本文已影响0人  高大石头

无意中翻到生信大神果子老师和洲更老师神仙大家的帖子,样本足够多时非参数检验就够用了(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")

上一篇下一篇

猜你喜欢

热点阅读