limour-blog.github.io/source/_posts/2022-07-30-【迁移】使用limma包进行差异基因分析.md

---

title: 【迁移】使用limma包进行差异基因分析 urlname: shi-yong-limma-bao-jin-hang-cha-yi-ji-yin-fen-xi date: 2022-07-30 19:13:36 index_img: https://api.limour.top/randomImg?d=2022-07-30 19:13:36

tags: ['DEG', 'limma']

虽然现在已经是高通量测序的时代，大家基本都是从counts矩阵出发，使用DESeq2进行差异表达分析，但是GEO和ArrayExpress上的仍有海量且持续更新的芯片数据，有时候也不可避免遇到一些FPKM格式乃至已经进行了z-score转换的数据，对于这些数据的分析，我们可以认为其在适当变换下(log2FPKM)，满足正态分布，那么仍可以使用limma直接进行分析。下面博主以E-MEXP-1422为例，写一份分析代码的demo。

require(limma)
f_DE_limma <- function(cts_bb, rowInfo, ControlN, TreatN, rm.NA=T, trend=T){
    # trend 表示先验方差是否与基因表达值的大小相关，False表示其为常数
    cts_b <- cts_bb[,c(ControlN, TreatN)]
    conditions <- c(rep("Control",length(ControlN)), rep("Treat",length(TreatN)))
    design <- model.matrix(~0+factor(conditions))
    colnames(design) <- levels(factor(conditions))
    rownames(design) <- colnames(cts_b)
    print(design)
    contrast.matrix <- makeContrasts('Treat-Control', levels = design)
    print(contrast.matrix)
    fit <- lmFit(cts_b, design) # 拟合线性模型
    fit2 <- contrasts.fit(fit, contrast.matrix) # 计算拟合系数和标准误差
    fit2 <- eBayes(fit2, trend=trend) # 通过经验贝叶斯方法估计统计量和logFC值
    tempOutput <- topTable(fit2, coef=1, n=Inf)
    if(!is.null(rowInfo)){tempOutput <- cbind(rowInfo[rownames(tempOutput),], tempOutput)}
    if(rm.NA){tempOutput <- na.omit(tempOutput)}
    tempOutput
}
# 经过 oligo::rma 标准化后提取出来的表达矩阵
data.exprs
# SDRF <- read.delim('E-MEXP-1422.sdrf.txt') 
# 从sdrf文件可知 AF15、AF16为PROX1 siRNA组
# AF6、AF14为GFP siRNA组
Ct1 <- c('AF6.CEL', 'AF14.CEL')
Tt1 <- c('AF15.CEL', 'AF16.CEL')
f_DE_limma(data.exprs, NULL, Ct1, Tt1, F)

蛋白质组学TCPA数据集

获取数据

• 进入TCPA的下载页面选择感兴趣的L4数据
• unzip TCGA-PRAD-L4.zip

清洗数据

f_dedup_IQR

tcpa <- read.csv('tmp/TCGA-PRAD-L4.csv')
type <- as.numeric(substr(tcpa$Sample_ID, 14, 15))
tcpa <- subset(tcpa, type < 10) # tp
rowNa <- substr(tcpa$Sample_ID,1, 12)
tcpa <- f_dedup_IQR(tcpa[-(1:4)],rowNa)
tcpa

TCGAbiolinks下载蛋白质组数据

之前通过tcpa下载过蛋白数据]，而TCGAbiolinks也有下载蛋白质组学数据的示例，后者看上去更全面一点。

下载数据

library(TCGAbiolinks)
query.rppa <- GDCquery(
    project = "TCGA-PRAD", 
    data.category = "Proteome Profiling",
    data.type = "Protein Expression Quantification"
)
GDCdownload(query.rppa) 
rppa <- GDCprepare(query.rppa)
saveRDS(rppa, 'PRAD_rppa.rds')

清洗数据

pMiss <- function(x){round(sum(is.na(x))/length(x),3)}
rppa <- rppa[apply(rppa, 1, pMiss) < 0.05, ]
rppa <- rppa[, apply(rppa, 2, pMiss) < 0.05]
sum(is.na(rppa))
rowInfo <- rppa[1:5]
rppa <- rppa[-(1:5)]
rppa <- rppa[, substr(colnames(rppa), 14, 16) == "01A"]
colnames(rppa) <- substr(colnames(rppa),1,12)
group <- readRDS('../fig5/tcga.predict.rds')
gRow <- intersect(colnames(rppa), rownames(group))
group <- group[gRow,]
rppa <- rppa[, colnames(rppa) %in% gRow]
Ct1 <- rownames(group)[group$group == 'Low Risk']
Tt1 <- rownames(group)[group$group == 'High Risk']

计算差异蛋白

r1 <- f_DE_limma(rppa, rowInfo, Ct1, Tt1, trend=F)
rownames(rppa) <- rowInfo$AGID
save(r1, rppa, file = 'PRAD_TCPA_DE.rdata')