芯片和高通量测序（HTS）：
对于基因芯片的差异表达分析而言，由于普遍认为其数据是服从正态分布，因此差异表达分析无非就是用t检验和或者方差分析应用到每一个基因上。目前在基因芯片的分析用的最多的就是limma。
高通量一次性找的基因多，于是就需要对多重试验进行矫正，控制假阳性。高通量测序(HTS)的read count普遍认为是服从泊松分布（当然有其他不同意见），不可能直接用正态分布的t检验和方差分析。
edgeR 和 DEseq2 这两个软件要求输入的文件中所有数值必须要是整数，所以只能用count数据；
limma 包做差异分析要求数据满足正态分布或近似正态分布，如基因芯片、TPM 格式的高通量测序数据。高通量测序得到的count数据不符合正态分布而服从泊松分布，对于count 数据来说，用limma包做差异分析，误差较大；
目前 DESeq2 可以处理符合泊松分布的数据，是分析高通量测序的主流。值得一提的是 edgeR 和 limma 是由一个团队开发的，算法有点过时了，DESeq2 目前使用频率较高。
1.数据处理
配对样本的lncRNA表达矩阵，供配对差异分析用
思路：选出与normal样本匹配的tumor样本，合并为新的表达矩阵
rm(list = ls())
library(stringr)
options(stringsAsFactors = F)
setwd("D:/BioInformationAnalyze/TCGAdata/CRC") # 设置工作目录（需修改）
load("TCGA-COAD_01.1.lncRNA_Expdata_count.Rdata") # 使用count数据
# 挑选出来normal的样本，然后匹配样本名
table(stringr::str_sub(colnames(exp),14,15)) #正常样本个数？
exp_nor <- exp[, str_sub(colnames(exp),14,15) == 11] # 正常样本
exp_tum <- exp[, str_sub(colnames(exp),14,15) != 11] # 肿瘤样本
patient <- str_sub(colnames(exp_nor),1,12) # 有正常样本的病人id（有配对样本的病人id）
table(str_sub(colnames(exp_tum),1,12) %in% str_sub(patient)) # 看看肿瘤样本有重复吗
# 匹配肿瘤样本
exp_tum <- exp_tum[, str_sub(colnames(exp_tum),1,12) %in% str_sub(patient)] # 找到所有配对肿瘤样本
exp_tum <- exp_tum[, str_sort(colnames(exp_tum))] # 配对肿瘤样本排序
exp_tum <- exp_tum[, !duplicated(str_sub(colnames(exp_tum),1,12))] # 排序去重，这样样本有01A和01B可选时，就会留下01A。
# 让正常样本顺序与肿瘤样本顺序一致,这样才方便写配对向量
tmp <- match(str_sub(colnames(exp_tum),1,12), str_sub(colnames(exp_nor),1,12))
exp_nor <- exp_nor[,tmp]
identical(str_sub(colnames(exp_tum),1,12), str_sub(colnames(exp_nor),1,12)) # 确认排序一致
exp_pair <- cbind(exp_nor, exp_tum) # 前一半是正常样本，后一半是肿瘤样本
# 保存
save(exp_pair, file = paste0(cancer_type, "_03.1.lncRNA_Pair.Rdata"))
2.lncRNA配对样本差异分析
library(stringr)
library(ggplotify)
library(patchwork)
library(cowplot)
library(edgeR)
rm(list = ls())
load("TCGA-COAD_01.1.lncRNA_Expdata_count.Rdata") # 表达数据的count数
load("TCGA-COAD_03.1.lncRNA_Pair.Rdata") # 配对样本数据
group_list <- rep(c("normal","tumor"), each = ncol(exp_pair)/2) %>%
factor(levels = c("normal","tumor")) # 设置分组类型
pairinfo <- rep(1:(ncol(exp_pair)/2), times = 2) # 设置配对信息
table(group_list)
# 开始用edgeR差异分析
dge <- DGEList(counts = exp_pair, group = group_list) # 创建EDGList——edgeR存储数据形式，edgeR存储数据在一个叫做DEGList的类似于list的对象中。
dge$samples$lib.size <- colSums(dge$counts)
dge <- calcNormFactors(dge) #数据标准化，calcNormFactors函数通过查找一组缩放因子来标准化文库的大小
design <- model.matrix( ~ pairinfo + group_list) # 构建分组矩阵
dge <- estimateDisp(dge, design) # 估计离散程度（Estimating dispersions）
fit <- glmQLFit(dge, design) # 广义线性模型 quasi-likelihood (QL) F-test拟合
fit2 <- glmQLFTest(fit)
result <- topTags(fit2) # 从测试对象中提取差异表达最多的基因（或序列标签），按 p 值或logFC排序。
DEG <- topTags(fit2, n = nrow(exp))
DEG <- as.data.frame(DEG) # 表中FDR的计算是根据假设检验的P-value进行校正而得到的。
#logFC_cutoff <- with(DEG, mean(abs(logFC)) + 2*sd(abs(logFC))) # 均数+2倍标准差
logFC_cutoff <- 1.5 # 设置logFC阈值
DEG$change <- as.factor(
ifelse(DEG$FDR < 0.05 & abs(DEG$logFC) > logFC_cutoff,
ifelse(DEG$logFC > logFC_cutoff ,"UP","DOWN"),"NOT")
) # 标基上下调基因
head(DEG) # 得到排序后的差异基因表
table(DEG$change)
lncRNA_DEG <- DEG
cg <- rownames(DEG)[DEG$change != "NOT"] # change gene：上调和下调基因的基因名
cgexp <- exp_pair[cg,] # 差异基因的表达量
table(head(DEG[cg, "change"], 580))
save(cg, file = paste0(cancer_type, "_03.2.lncRNA_cg.Rdata"))
画PCA图、热图、火山图
library(ggplot2) library(tinyarray) exp_pair[1:4, 1:4] # PCA图 dat <- log(exp_pair + 1) pca.plot <- draw_pca(dat, group_list); pca.plot # 热图、火山图 x <- log2(exp_pair[cg,] + 1) heatmap.plot <- draw_heatmap(x, group_list, n_cutoff = 2.6, cluster_cols = F) ; heatmap.plot volcano.plot <- draw_volcano(lncRNA_DEG, logFC_cutoff = 1.5, pkg = 2, symmetry = T, adjust = T) ; volcano.plot
