作者：oct
审稿：童蒙
编辑：amethyst
配对差异分析经常应用在肿瘤配对样本的差异分析中。在肿瘤研究中，通常需要进行肿瘤组织与相邻正常组织之间的差异分析，以便研究肿瘤组织的特异性。那么在差异分析时，配对和非配对的差异分析有何区别，该如何选择？经查阅文献，一篇名为Differential Expression of miRNAs in Colorectal Cancer: Comparison of Paired Tumor Tissue and Adjacent Normal Mucosa Using High-Throughput Sequencing的文章或许会给大家带来一些新的思路。
背景介绍
MicroRNA (miR) 是长度为 18-25 个核苷酸的小型非编码 RNA 分子，于 1990 年代初首次在秀丽隐杆线虫中发现。它们通过改变不同细胞过程（如分化、增殖、存活和凋亡）中的基因表达来维持体内平衡 。据估计，超过 10% 的编码人类基因的蛋白质可能受这些机制的调控 。miRBase 数据库中记录的人类 miR 数量超过一千 。研究表明，miRs 可能在不同的人类癌症中失调，因此充当肿瘤抑制基因或癌基因。它们可能是诊断或预后的潜在生物标志物，并作为癌症特异性治疗的潜在靶点。
材料与方法
01 样本选择
使用 Illumina 高通量测序技术研究 miR 表达的肿瘤特异性变化从8个患者的手术标本中收集正常粘膜和肿瘤组织，从而产生一组独特的成对样品，7个患者的肿瘤细胞含量大于60%，其中一例为非典型神经内分泌肿瘤 (NET)，其他为腺癌。
02 数据分析
1 数据处理：获得高通量测序的fastq数据，使用FASTX-Toolkit去接头。测序数据与 hg18 基因组参考比对，允许一个错误匹配。使用 miRanalyzer 进一步处理测序数据。该工具允许从 miRBase数据库中识别经过验证的 miR，并包括用于预测新 miR 的机器学习算法。
2 miR差异分析：使用edgeR与Deseq进行miR差异分析，两种工具都利用负二项式分布对每个miR的读取计数进行建模，并实现了对计数进行归一化的方法。
- 非配对差异分析：Deseq
- 配对差异分析：edgeR
主要结果
通过韦恩图可以看出，在两种方法中，有37个miR在配对和非配对分析中共同检出，有 81 个 miR 在非配对分析中未鉴定，证明非配对差异分析相比配对分析更保守，它不需考虑患者之间的基线差异。
差异分析R代码
1、加载R包和数据
library(DESeq)
library(edgeR)
setwd("/Users/Julian/Documents/Prosjekt/PLOS/Dataset/")
2、数据处理
targets <- read.delim(file="Targets.txt", stringsAsFactors=FALSE) #读取数据
cts <- readDGE(targets)
countsTable <- cts$counts #读取count
colnames(countsTable) <- c("T1", "T2", "T3", "T4", "T5", "T6", "T7", "T8", "N1", "N2", "N3", "N4", "N5", "N6","N7", "N8") #列为样本
conds <- c("Ne","Ac","Ac","Ac","Ac","Ac","Ac","Ac","No","No","No","No","No","No","No","No") #分组信息
cds <- newCountDataSet(countsTable, conds) #构建cds对象
cds <- estimateSizeFactors(cds) #归一化
cds <- estimateVarianceFunctions(cds) #方差估计
3、使用 DESeq 计算正常黏膜 (No) 与腺癌 (Ac) 中失调的 miRs
resNoAc <- nbinomTest(cds, "No", "Ac") #对负二项模型进行T检验
resNoAcSig <- resNoAc[resNoAc$padj<.1,] #设置显著阈值，padj值小于0.1
resNoAcSig <- resNoAcSig[order(resNoAcSig$padj),]
subset(resNoAcSig, select=c(1,5,6,8))
4、使用 DESeq 计算正常黏膜 (No) 与神经内分泌肿瘤 (Ac) 中失调的 miRs
resNoNe <- nbinomTest(cds, "No", "Ne")
resNoNeSig <- resNoNe[resNoNe$padj<.1,]
resNoNeSig <- resNoNeSig[order(resNoNeSig$padj),]
subset(resNoNeSig, select=c(1,5,6,8))
5、方差函数拟合
diagForT <- varianceFitDiagnostics (cds1, "T")
smoothScatter( log10(diagForT$baseMean), log10(diagForT$baseVar) )
lines( log10(fittedBaseVar) ~ log10(baseMean), diagForT[ order(diagForT$baseMean), ], col="red" )
abline(0,1,lty=2)
6、在配对的正常粘膜与腺癌中使用 edgeR 计算失调的 miR
targetsPaired <- read.delim(file="TargetsPaired.txt", stringsAsFactors=FALSE)
d <- readDGE(targetsPaired)
colnames(d) <- c("T2", "T3", "T4", "T5", "T6", "T7", "T8", "N2", "N3", "N4", "N5", "N6", "N7", "N8")
patient <- factor(c(2, 3, 4, 5, 6, 7, 8, 2, 3, 4, 5, 6, 7, 8)) #设置配对信息，相同数字为配对样本
design <- model.matrix(~patient + d$samples$group)
rownames(design) <- rownames(d$samples)
design[,8] <- c(1,1,1,1,1,1,1,0,0,0,0,0,0,0) # 设置分组信息
colnames(design)[8] <- "tumor"
d <- estimateGLMCommonDisp(d, design)
glmfit.d <- glmFit(d, design, dispersion = d$common.dispersion)
lrt.d <- glmLRT(d, glmfit.d, coef = 8)
options(digits = 4)
topTags(lrt.d, n=118)
7、上调/下调 miR 统计
sum(lrt.d$table$logFC > 0)
sum(lrt.d$table$logFC < 0)
top <- topTags(lrt.d,n=118)
sum(top$table$logFC > 0)
sum(top$table$logFC < 0)
总结
- edgeR配对分析模型中省略了回归交互项。配对差异分析时，对单个患者的比较组不感兴趣，相反，感兴趣的是一组患者的比较组的平均差异。
- 配对样本进行配对差异分析或非配对差异分析，结果大不相同，所以在选择差异分析方法时，需要考虑实验设计的目的。
- Deseq2和edgeR两种R包都可以进行配对差异分析，Deseq则无法进行配对差异分析。
参考资料
- edgeR: differential expression analysis of digital gene expression data；
- http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#can-i-use-deseq2-to-analyze-paired-samples；
- Hamfjord J, Stangeland AM, Hughes T, Skrede ML, Tveit KM, Ikdahl T, Kure EH. Differential expression of miRNAs in colorectal cancer: comparison of paired tumor tissue and adjacent normal mucosa using high-throughput sequencing. PLoS One. 2012;7(4):e34150. doi: 10.1371/journal.pone.0034150. Epub 2012 Apr 17. PMID: 22529906; PMCID: PMC3328481.；
- Roukos DH (2010) Novel clinico-genome network modeling for revolutionizing genotype-phenotype-based personalized cancer care. Expert Rev Mol Diagn 10: 33–48.
