最近安排学徒做文献图表复现，其中一个表达量芯片和测序项目都是同样的处理和对照，所以让学徒做一下这两个表达矩阵的差异分析，比较一下不同技术是否有比较好的吻合。
其中测序是：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE64486
可以看到是8个KD病人，他们一半处理一半为处理：
代码语言：javascript复制GSM1572233 KD1 - untreated
GSM1572234 KD2 - untreated
GSM1572235 KD3 - untreated
GSM1572236 KD4 - untreated
GSM1572237 KD5 - treated
GSM1572238 KD6 - treated
GSM1572239 KD7 - treated
GSM1572240 KD8 - treated
作者给出来的每个样品一个独立的表达量矩阵txt格式的文本文件：
代码语言：javascript复制GSM1572233_8502X3_case_union_COUNTS.txt.gz 175.2 Kb
GSM1572234_10229X2_case_union_COUNTS.txt.gz 186.7 Kb
GSM1572235_10229X1_case_union_COUNTS.txt.gz 189.7 Kb
GSM1572236_9374X3_case_union_COUNTS.txt.gz 177.3 Kb
GSM1572237_8898X1_case_union_COUNTS.txt.gz 189.5 Kb
GSM1572238_8898X2_case_union_COUNTS.txt.gz 188.2 Kb
GSM1572239_8502X4_case_union_COUNTS.txt.gz 186.7 Kb
GSM1572240_10229X3_case_union_COUNTS.txt.gz 191.3 Kb
其对应的文章是：The transcriptional profile of coronary arteritis in Kawasaki disease. BMC Genomics 2015 Dec ，作者做的是这些KD病人去跟正常人的对比：
这些KD病人去跟正常人的对比
做了三次差异分析，都是病人和正常人的对比，并没有做4个处理的病人和4个未处理病人的差异，所以我安排学徒做了这个。
有意思的是，默认情况下我们会读取表达量矩阵后看看常见的管家基因是否是高表达，来确定矩阵的准确性。如下所示的代码，需要自己下载GSE64486_RAW.tar文件，然后解压哦，接下来就可以读取解压后的文件夹里面的全部的 txt文件 ：
代码语言：javascript复制fs=list.files('GSE64486_RAW/',pattern = 'case',full.names = T)
fs
library(data.table)
gid=fread(fs[1],data.table = F)[,1]
head(gid)
rawcount = do.call(cbind,
lapply(fs, function(x){
fread(x,data.table = F)[,2]
}))
rawcount[1:4,1:4]
keep <- rowSums(rawcount>0) >= floor(0.75*ncol(rawcount))
table(keep)
rawcount=rawcount[keep,]
gid=gid[keep]
library(AnnoProbe)
ids=annoGene( gid ,
ID_type = 'ENSEMBL',species = 'human'
)
head(ids)
ids=ids[!duplicated(ids$SYMBOL),]
rawcount=rawcount[match(ids$ENSEMBL , gid),]
rownames(rawcount) = ids$SYMBOL
rawcount[1:4,1:4]
rawcount['GAPDH',]
cg = log2(edgeR::cpm(rawcount) 1)['GAPDH',]
gp = rep(c('untreated','treated'),each=4)
t.test(cg ~ gp)
可以看到：
代码语言：javascript复制> rawcount[1:4,1:4]
[,1] [,2] [,3] [,4]
WASH7P 6 5 10 0
AL627309.1 27 28 13 4
AL627309.5 46 101 88 29
AL627309.4 0 6 11 5
> rawcount['GAPDH',]
[1] 720 873 2374 893 2655 6928 2687 12173
表达量矩阵本身没有问题，表达量比较高就是两个超级出名的非编码基因，以及线粒体基因等等：
合理的表达量矩阵
符合常识，所以我非常确定作者的转录组上游定量流程是ok的。但是它确实在GAPDH这个基因上面，有表达量的差异，仅仅是原始counts值就出现了肉眼可见的巨大差异，哪怕是归一化后也是如此 ：
代码语言：javascript复制> log2(edgeR::cpm(rawcount) 1)['GAPDH',]
[1] 8.176485 6.939699 7.703584 6.935478 8.039263 8.904722 8.204658 9.174328
> t.test(cg ~ gp)
Welch Two Sample t-test
data: cg by gp
t = 2.7911, df = 5.9258, p-value = 0.03195
sample estimates:
mean in group treated mean in group untreated
8.580743 7.438812
当然了，我们不能是仅仅是靠肉眼，靠单个基因的检验 ，我们还是使用常见的转录组测序表达量矩阵的差异分析方法：
代码语言：javascript复制exprSet = rawcount
group_list = gp
# 加载包
library(DESeq2)
# 第一步，构建DESeq2的DESeq对象
colData <- data.frame(row.names=colnames(exprSet),group_list=group_list)
dds <- DESeqDataSetFromMatrix(countData = exprSet,colData = colData,design = ~ group_list)
# 第二步，进行差异表达分析
dds2 <- DESeq(dds)
# 提取差异分析结果，trt组对untrt组的差异分析结果
tmp <- results(dds2,contrast=c("group_list","treated","untreated"))
DEG_DESeq2 <- as.data.frame(tmp[order(tmp$padj),])
head(DEG_DESeq2)
# 去除差异分析结果中包含NA值的行
DEG_DESeq2 = na.omit(DEG_DESeq2)
DEG_DESeq2['GAPDH',]
可以看到，这个时候它仍然是统计学显著，也就是说pvalue是 0.05以下，但是矫正后padj就不显著啦，所以它被排除到差异上下调基因之外了哦。
另外，它的log2FoldChange 已经是 1.14而已，我们如果稍微提高一下阈值，比如 1.5，它这个基因也没有统计学显著的差异。我们肉眼看到的差异其实是转录组测序本身的技术噪音，它一下子对两万多个基因进行定量，肯定是有一些哪怕是本身并不会在两个分组有明显差异表达量的基因也会表现出来差异，但是差异是否在可以接受的范围就需要严格的统计学检验了哦。
眼尖的小朋友可能会注意到，其实这个探索，忽略了一个很重要的质量控制，就是三张图。我在生信技能树的教程：《你确定你的差异基因找对了吗？》提到过，必须要对你的转录水平的全局表达矩阵做好质量控制，最好是看到标准3张图：
- 左边的热图，说明我们实验的两个分组，normal和npc的很多基因表达量是有明显差异的
- 中间的PCA图，说明我们的normal和npc两个分组非常明显的差异
- 右边的层次聚类也是如此，说明我们的normal和npc两个分组非常明显的差异
如果分组在3张图里面体现不出来，实际上后续差异分析是有风险的。这个时候需要根据你自己不合格的3张图，仔细探索哪些样本是离群点，自行查询中间过程可能的问题所在，或者检查是否有其它混杂因素，都是会影响我们的差异分析结果的生物学解释。
学徒作业
针对这个GSE64486数据集进行文章的3次差异分析，并且和作者给出来的上下调基因进行对比，看看十年前的差异分析和现在的差异分析，是否会有很严重的区别！
代码语言：javascript复制GSE64486_case_vs_control_foldgenes_broad.txt.gz 1.3 Mb
GSE64486_treated_vs_control_foldgenes_broad.txt.gz 1.1 Mb
GSE64486_untreated_vs_control_foldgenes_broad.txt.gz 1.2 Mb 
