上下调基因数量不平衡？再探！

dat <- read.table("STAD.rnaseqv2__illuminahiseq_rnaseqv2__unc_edu__Level_3__RSEM_genes_normalized__data.data.txt",
                  sep = "\t",
                  header = TRUE)
TCGA_id<- colnames(dat)[-1]
library(tidyverse)
library(stringr)
split_list <- str_split_fixed(TCGA_id,"\\.",7)
apply(split_list,2,table)
# [[1]] project
# TCGA 
# 450 
# [[4]] sample type
# 01A 01B 11A 
# 413   2  35
# 说明有35个normal对照样本 415个tumor样本

# participants
table(split_list[,3]) %>% table()
# 1   2 
# 386  32 
# 说明有32个参与者被收集了两种组织

EBV_groups_dat <- readxl::read_xlsx("STAD_Master_Patient_Table_20140207.xlsx")
EBV_groups_dat$group <- ifelse(EBV_groups_dat$`CIMP Category`=="GASTRIC-EBV-CIMP","EBVa","EBVn")
table(EBV_groups_dat$group)
# EBVa EBVn 
# 28  267
head(EBV_groups_dat)
# # A tibble: 6 × 47
# `TCGA barcode` `Lauren Class` `Intestinal Type Subclass` `Signet Ring` `WHO Class`     `Pathologic T`
# <chr>          <chr>          <chr>                      <chr>         <chr>           <chr>   
#   1 TCGA-B7-5816   Diffuse        NA                         0             Poorly_Cohesive T4a   
# 2 TCGA-B7-5818   Intestinal     NA                         NA            NA              T2    
# 3 TCGA-BR-4183   Intestinal     Tubular                    NA            Tubular         T4a   
# 4 TCGA-BR-4184   Intestinal     Tubular                    NA            Tubular         T4a   
# 5 TCGA-BR-4187   Diffuse        NA                         1             Poorly_Cohesive TX    
# 6 TCGA-BR-4188   Diffuse        NA                         1             Poorly_Cohesive TX   

new_TCGA_id <- paste(split_list[,1],split_list[,2],split_list[,3],sep = "-")
head(new_TCGA_id)
# "TCGA-3M-AB46" "TCGA-3M-AB47" "TCGA-B7-5816" "TCGA-B7-5818" "TCGA-B7-A5TI" "TCGA-B7-A5TJ"
colnames(dat) <- c("gene_id",new_TCGA_id)
dat <- dat[-1,]
rownames(dat) <- dat$gene_id
dat$gene_id <- NULL

# 看看2014的样本信息在不在这里
table(EBV_groups_dat$`TCGA barcode` %in% new_TCGA_id)
# FALSE  TRUE 
# 15   280 
# 看起来大部分都在
# 重点看EBVa在不在
table(EBV_groups_dat$`TCGA barcode`[EBV_groups_dat$group=="EBVa"] %in% new_TCGA_id)
# FALSE  TRUE 
# 1    27 
# 太好了大部分都在

# 提取对应数据
# 先确保是肿瘤样本的
sub_dat <- dat[,ifelse(split_list[,4]=="11A",FALSE,TRUE)]
dim(sub_dat)  # 20531   415
table(duplicated(colnames(sub_dat)))
# FALSE 
# 415
# 获取2014EBV分型样本
EBVa <- EBV_groups_dat$`TCGA barcode`[EBV_groups_dat$group=="EBVa"]
EBVn <- EBV_groups_dat$`TCGA barcode`[EBV_groups_dat$group=="EBVn"]
EBVa_dat <- sub_dat[,which(colnames(sub_dat) %in% EBVa)]  # 27
EBVn_dat <- sub_dat[,which(colnames(sub_dat) %in% EBVn)]  # 250
EBV_dat <- cbind(EBVa_dat,EBVn_dat)
# 发现读进来的表达量值为character
# numeric化
EBV_dat_num <- apply(EBV_dat, 2, as.numeric)
rownames(EBV_dat_num) <- rownames(EBV_dat)

######差异表达分析######

# 作者直接使用的R’s built-in wilcox.test function
# Wilcoxon秩和检验进行非参数检验

# 输入EBV_dat_num 和d
d <- 27
n <- length(colnames(EBV_dat_num))
deg_pvalue <- apply(EBV_dat_num, 1, function(x){
  # wilcoxon
  # print(x)
  res <- wilcox.test(x[1:d],x[(d+1):n],
                     exact = FALSE)
  return(res$p.value)
})
deg_log2fc <- apply(EBV_dat_num, 1, function(x){
  fc <- mean(x[1:d]) / mean(x[(d+1):n])
  return(log2(fc))
})
# > table( qvalue(deg_pvalue)$qvalues < 0.05)
# FALSE  TRUE 
# 12235  8039 
deg_qvalue <- qvalue(deg_pvalue)
# 放一起
dat_to_draw <- data.frame(deg_log2fc,deg_qvalue$qvalues)
# 去掉全部是0的基因
dat_to_draw <-  na.omit(dat_to_draw)

dat_to_draw$regulation <- ifelse(
  dat_to_draw$deg_log2fc<(-0.4) & dat_to_draw$deg_qvalue.qvalues<0.05, "down",
  ifelse(dat_to_draw$deg_log2fc>0.4 & dat_to_draw$deg_qvalue.qvalues<0.05,"up","normal")
)
table(dat_to_draw$regulation)
# down normal     up 
# 4416  14518   1340 

#####绘图######

# 火山图
library(EnhancedVolcano)
EnhancedVolcano(dat_to_draw,x='deg_log2fc',y='deg_qvalue.qvalues',
                lab=rownames(dat_to_draw),
                pCutoff = 0.05,
                FCcutoff = 0.4)

volcano_p <- ggplot(data=dat_to_draw, aes(x=deg_log2fc, y=-log10(deg_qvalue.qvalues),color=regulation)) +  # 按regulated分组颜色
  geom_point(alpha=0.5, size=1.8) +  # alpha控制点的透明度 size大小
  theme_set(theme_set(theme_bw(base_size=20))) + 
  xlab("log2FC") + ylab("-log10(Pvalue)") +
  scale_colour_manual(values = c('blue','black','red'))+theme_bw()
volcano_p

> length(EBV_dat_num[1,][d+1:n])
[1] 277
> length(EBV_dat_num[1,][(d+1):n])
[1] 250

# 整理到函数里
degs_analysis_draw <- function(mat,d){
  print(dim(mat))
  # 输入EBV_dat_num 和d
  n <- length(colnames(mat))
  deg_pvalue <- apply(mat, 1, function(x){
    # wilcoxon
    # print(x)
    res <- wilcox.test(x[1:d],x[(d+1):n],
                       exact = FALSE)
    return(res$p.value)
  })
  deg_log2fc <- apply(mat, 1, function(x){
    fc <- mean(x[1:d]) / mean(x[(d+1):n])
    return(log2(fc))
  })
  # > table( qvalue(deg_pvalue)$qvalues < 0.05)
  # FALSE  TRUE 
  # 12235  8039 
  deg_qvalue <- qvalue(deg_pvalue)
  # 放一起
  dat_to_draw <- data.frame(deg_log2fc,deg_qvalue$qvalues)
  # 去掉全部是0的基因
  dat_to_draw <-  na.omit(dat_to_draw)
  
  dat_to_draw$regulation <- ifelse(
    dat_to_draw$deg_log2fc<(-0.4) & dat_to_draw$deg_qvalue.qvalues<0.05, "down",
    ifelse(dat_to_draw$deg_log2fc>0.4 & dat_to_draw$deg_qvalue.qvalues<0.05,"up","normal")
  )
  print(table(dat_to_draw$regulation))
  
  #####绘图######
  
  # 火山图
  library(EnhancedVolcano)
  EnhancedVolcano(dat_to_draw,x='deg_log2fc',y='deg_qvalue.qvalues',
                  lab=rownames(dat_to_draw),
                  pCutoff = 0.05,
                  FCcutoff = 0.4)
  
  volcano_p <- ggplot(data=dat_to_draw, aes(x=deg_log2fc, y=-log10(deg_qvalue.qvalues),color=regulation)) +  # 按regulated分组颜色
    geom_point(alpha=0.5, size=1.8) +  # alpha控制点的透明度 size大小
    theme_set(theme_set(theme_bw(base_size=20))) + 
    xlab("log2FC") + ylab("-log10(Pvalue)") +
    scale_colour_manual(values = c('blue','black','red'))+theme_bw()
  volcano_p
  
  print("----------end-------------")
  
}

#########抽样测试########
# 先看看函数能不能用
degs_analysis_draw(mat = EBV_dat_num,d = 27)
# 27 vs 250
# sample size
index <- seq(28,277)
for (sample_size in seq(50,250,by=50)){
  # EBVn 50 100 .. 200 250
  print(sample_size)
  cur_mat <- EBV_dat_num[,c(seq(1,27),sample(index, sample_size))]
  degs_analysis_draw(cur_mat, d)
  # break
}

[1] 50
[1] 20531    77

  down normal     up 
  3823  15037   1357 
[1] "----------end-------------"
[1] 100
[1] 20531   127

  down normal     up 
  4005  14826   1411 
[1] "----------end-------------"
[1] 150
[1] 20531   177

  down normal     up 
  4400  14501   1355 
[1] "----------end-------------"
[1] 200
[1] 20531   227

  down normal     up 
  4320  14553   1387 
[1] "----------end-------------"
[1] 250
[1] 20531   277

  down normal     up 
  4416  14518   1340 
[1] "----------end-------------"

# 散点图
num_degs <- data.frame(
  down = c(3823,4005,4400,4320,4416),
  up = c(1357,1411,1355,1387,1340)
)
# basic scatterplot
scatter_plot <- ggplot(num_degs, aes(x=down, y=up)) + 
  geom_point()
scatter_plot

# sample distribution
index <- seq(28,277)
for (test in seq(1,5)){
  # EBVn 50 50 .. 50 50
  # 选小一点，这样可以把样本异质性凸显出来
  # print(sample_size)
  sample_size <- 50
  cur_mat <- EBV_dat_num[,c(seq(1,27),sample(index, sample_size))]
  degs_analysis_draw(cur_mat, d)
  # break
}
# 散点图
num_degs <- data.frame(
  down = c(2939,2932,2996,3609,3396),
  up = c(1199,1197,1217,1329,1234)
)
# basic scatterplot
scatter_plot <- ggplot(num_degs, aes(x=down, y=up)) + 
  geom_point()
scatter_plot