################################################################################
# 1. TCGA 表達數(shù)據(jù)與導(dǎo)入
## https://xenabrowser.net/ 首先從xena上選好癌癥類型，下載表達數(shù)據(jù)， download as "expr.tsv", 去除無表達值的樣品
################################################################################

rm(list = ls())
library (stringr)
library (tidyverse)

## 讀取表達信息
expr <- read.table("expr.tsv",header = T, stringsAsFactors = F, quote = "",check.names = F)

##查看數(shù)據(jù)結(jié)構(gòu)
expr[1:3,1:4]

## 查看腫瘤、正常樣品： FALSE=Normal, TRUE=Tumor
table(str_sub(expr$sample,14,15)<10)

################################################################################
# 2. 從TCGA表達矩陣?yán)锊鸪雠鋵颖?################################################################################
##配對樣本的共同點是病人ID（前12位）一致，根據(jù)這個來匹配即可。

##先拆分成tumor和normal兩個矩陣，根據(jù)ID的14和15位來拆， 14-15>10: Normal, 414-15<=10: Tumor 
expr_nor = expr[str_sub(expr$sample,14,15) > 10,]
expr_tum = expr[str_sub(expr$sample,14,15)<=10,]

## 有normal樣本的病人的ID，是normal組樣本ID的前12位，也是expr_nor的行名。

## 然后將tumor矩陣中與patient相匹配的樣本名選出來，能匹配到的就會被挑出來。%in% or match 都是R語言里的神技能：match 更好，排序一致。
k<- match(str_sub(expr_nor$sample,1,12),str_sub(expr_tum$sample,1,12))
expr_tum_paired <- expr_tum[k,]

##可以把表達矩陣拼回去：
expr_paired <- rbind(expr_nor,expr_tum_paired)
expr_paired$Type <-ifelse(str_sub(expr_paired$sample,14,15) > 10, "Normal", "Tumor")

#dim(expr_paired)
save(expr_paired,file="expr_paired.Rdata")

################################################################################
# 3. 成對樣本表達差異分析：ggpaired 函數(shù)
################################################################################

## 3.1 單基因分析
rm(list = ls())
load("expr_paired.Rdata")

library(ggpubr)
library(ggstatsplot)
library(ggsci)
gene="YTHDF1" ## 選擇要查看的基因
ggpaired(expr_paired, x="Type", y=gene, color = "Type", line.color = "lightblue", point.size = 1.5,
         font.label = list(size = 16, color = "black"),
         ggtheme=theme_bw(),repel = TRUE,
         Tile=gene,
         ylab= "Expression(Log2(nCount+1))",
         xlab=" Tissue Type",
         line.size = 0.5, palette = "aaas")+ stat_compare_means(paired = TRUE)

ggsave(file=paste0(gene,".png"),width = 300,height = 450)

## 3.2 多基因批量基因分析、組圖

library(ggplotify)

p=list()

gene_nb <- length(colnames(expr_paired))-2
for ( i in 1:gene_nb){
  gene<-colnames(expr_paired)[i+1]
  p[[i]] <-ggpaired(expr_paired, x="Type", y=gene, color = "Type", 
                    line.color = "lightblue", point.size = 1.5,
                    font.label = list(size = 16, color = "black"),
                    ggtheme=theme_bw(), repel = TRUE,
                    title=gene,
                    ylab= "Expression(Log2(nCount+1))",
                    xlab=" Tissue Type",
                    line.size = 0.5, palette = "aaas")
    + stat_compare_means(paired = TRUE)
  ggsave( p[[i]], file=paste0(gene,".png"),width = 300,height = 450)
  
}

### 組圖
cowplot::plot_grid(plotlist = p,
                   labels =LETTERS[1:gene_nb],
                   ncol=2)

## 3.3 單基因多癌癥種-成對分析（待補充）
### 同上下載不同癌癥的成對樣品所有基因的表達數(shù)據(jù)，整理好數(shù)據(jù)
### 批量做圖、組圖