测序上游分析系列：
mRNA-seq转录组二代测序从raw reads到表达矩阵：上中游分析pipeline
miRNA-seq小RNA高通量测序pipeline：从raw reads，鉴定已知miRNA-预测新miRNA，到表达矩阵【一】
miRNA-seq小RNA高通量测序pipeline：从raw reads，鉴定已知miRNA-预测新miRNA，到表达矩阵【二】
其他文章系列：ggplot2作图篇：
（1）基于ggplot2的RNA-seq转录组可视化：总述和分文目录
（2）测序结果概览：基因表达量rank瀑布图，高密度表达相关性散点图，注释特定基因及errorbar的表达相关性散点图绘制
（3）单/多个基因在组间同时展示的多种选择：分组小提琴图、分组/分面柱状图、单基因蜂群点图拼图的绘制
（4）配对样本基因表达趋势：优化前后的散点连线图+拼图绘制
上次说到了TCGA高通量转录组测序数据的下载，整合与临床数据的清洗。接下来的必要流程必定和肯定是癌组织与癌旁组织的差异表达分析了。
R语言中有两个封装R包广泛用于差异基因的筛选：DESeq2和EdgeR。由于测序数据属于负二项分布，而且不同测序批次间总体深度会有区别，因此他们不能直接使用t检验进行差异分析！封装包帮我们解决了这个问题：本教程使用DESeq2进行筛选，不需要知道原理，操作过程非常简洁，就可以得到结果。
本例仍使用之前教程选择的的TCGA-LUSC白色人种肺鳞癌表达谱数据。下载整合操作见前文：从TCGA数据库下载并整合清洗高通量肿瘤表达谱-临床性状数据
绝对的TIPS：DESeq2以及EdgeR包要求的输入文件为【未经标准化的原始read counts表达矩阵】！DESeq2会对原始counts进行标准化（该方法优于CPM/TPM/FPKM/RPKM等方法），去除不同测序批次导致的批次效应（batch effect），因此不能使用FPKM等已经标准化过的矩阵。
所以在TCGA-portal中需要下载的文件是HTSeq-counts.txt。另外，上篇教程提到过的，counts文件在镜像操作的同时，需要在获得最初的整合count matrix时去除末尾5行总结性（mei）注（yong）释数据。
在经过整合筛选后，我们从之前的教程中获得（或者说DESeq2需要）：
（1）condition_table: containing three columns named 'TCGA_IDs', 'sample_type', and 'submitter_id'. It will be used by DESeq2 since it contains grouping information of each sample (condition=normal/cancer).
（2）mRNA_exprSet: the filtered read counts matrix whose columns represent samples and whose rows represent genes. The colnames are TCGA_IDs of samples identical with that in condition_table, and rownames are gene symbols.
注：从上个教程到这里，别忘了先把mRNA_exprSet的第一列'gene_name'并入rownames！否则后续分析columns数对不上。
rownames(mRNA_exprSet) <- mRNA_exprSet$gene_name
mRNA_exprSet <- mRNA_exprSet[,-1]
原本392个样本保留386个，其中包含344个癌组织样本，以及42个癌旁组织样本。
分析部分
- 下载并加载DESeq2包，作图需ggplot2包，热图所需pheatmap包，小提琴图所需vioplot包和管道符所在的包dplyr。
install.packages("DESeq2")
install.packages('ggplot2')
install.packages('pheatmap')
install.packages('vioplot')
install.packages('dplyr')
library("DESeq2")
library('ggplot2')
library('pheatmap')
library('vioplot')
library('dplyr')
- 将condition_table的分组信息进行因子化，并设置reference。
DESeq2要求condition_table的格式为：至少含一列样本的分组信息，且rownames需要是可以与表达矩阵colnames匹配的样本名称。此外grouping information需要是factor格式。
我们有必要让R知道我们设置组别的意义。如果不手动设置分组的reference（差异分析中的分母，也就是cancer与normal比较中的normal），DESeq2默认按照分组的首字母顺序设置reference，显然是不准确的。
condition_table <- condition_table[, c('TCGA_IDs', 'sample_type')]
rownames(condition_table) <- condition_table$TCGA_IDs
condition_table[,2] <- as.factor(condition_table[,2])
condition_table$sample_type <- relevel(condition_table$sample_type, ref = 'normal')
- 过滤counts matrix中read counts较低的基因。
如果某基因在样本间的整体表达水平很低，则小幅度的read counts扰动即可造成较大的fold change差异，但这种差异属于noise，并不具有生物学意义。所以需要设置一定的threshold将这些基因去除，减少假阳性。
过滤的原则：原理类似，具体操作各自不同。有分析使用如下原则：若某基因在所有样本中的read counts之和小于样本总数（样本均read counts<1）则去除。
本例使用相对更苛刻的原则：若某基因在大于80%的样本中read counts数<10，则将其滤除。
qualified_genes <- c()
for (genes_in_sheet in rownames(mRNA_exprSet)) {
qualification <- mRNA_exprSet[genes_in_sheet,] <= 10
if (sum(qualification) < 0.8*length(mRNA_exprSet)) {
qualified_genes <- append(qualified_genes, genes_in_sheet)
}
}
mRNA_expr_for_DESeq <- mRNA_exprSet[qualified_genes,]
- 将已经调整好的输入文件输入DESeq2，一步自动分析。
封装函数帮我们一步完成了以下步骤：
（1）estimation of size factor of each sample for normalization.
（2）estimation of dispersion.
（3）negative binomial GLM fitting and Wald tests.
由于TCGA样本量比较大，所以这一步会需要较长时间。作为参考，macbook Rstudio跑392个样本用了20多分钟。
#colData parameter needs our condition_table to recognize the samples, and design parameter needs the grouping factors to identify the classification.
dds <- DESeqDataSetFromMatrix(mRNA_expr_for_DESeq, colData = condition_table, design = ~ sample_type)
dds_DE <- DESeq(dds)
- 以表格的形式，按照自己的标准分别输出所有、上调和下调的差异基因。
dds_DE对象中包含了差异分析的所有结果。我们需要自己设定cut-off值进行输出：
（1）首先在result函数中设定alpha第一类错误率（即p值threshold）, 以及再次申明比较类型（contrast参数设置如下：c(condition_table中分组变量名，分子组名，分母组名)）并按校正后p值（padj，p.adjusted, q-value, 也是False Discovery Rate, FDR, 通常可以根据差异基因数量设置0.1，0.05或0.01）大小排序。
（2）subset函数设定校正后p值，以及log2FoldChange阈值（相对表达差异的绝对值超过该阈值视为差异表达）。
本例中样本较多，为控制差异基因数量，使用padj < 0.05和|log2FoldChange| > 2（FoldChange大于4）作为cut-off值。
TIPS：测序数据差异分析需要用校正p值FDR而不是p值作为cut-off！用于控制假阳性率。
注：最后本地保存差异基因表我选择以txt格式保存！因为csv等格式常默认用excel打开，随后某些基因名会被Excel自动转换为年月...血泪教训啊！
res_DE <- results(dds_DE, alpha=0.05,contrast=c("sample_type","cancer","normal"))
resOrdered <- res_DE[order(res_DE$padj),]
resSig <- subset(resOrdered,padj < 0.05)
resSigAll <- subset(resSig, (log2FoldChange < (-2)|log2FoldChange > 2))
resSigUp <- subset(resSig,log2FoldChange >2)
resSigDown <- subset(resSig,log2FoldChange <(-2))
#export the datasheets.
dir.create('/Users/zx/Desktop/TCGA/TCGA-LUSC-white/counts/tutorial/Sig_genes')
setwd('/Users/zx/Desktop/TCGA/TCGA-LUSC-white/counts/tutorial/Sig_genes')
write.table(resSigAll, 'DE_genes_all.txt', quote = F, sep = '\t')
write.table(resSigUp, 'DE_genes_up.txt', quote = F, sep = '\t')
write.table(resSigDown, 'DE_genes_down.txt', quote = F, sep = '\t')
所以我们就获得了代表所有/上调/下调差异基因的三个csv表格！我们只需要关注log2FoldChange, padj和表达趋势即可。
可视化部分
- 绘制各样本标准化后read counts的箱线图。
用于可视化表达数据标准化的程度。若各样本normalized read counts的分布水平相对稳定，则可以认为标准化的操作以及后续的分析是可信有效的。
使用vst函数对数据进行variance stablizing transformation，并获得转换后的表达矩阵。vst转换可以使表达矩阵适用于可视化、聚类以及机器学习等操作分析。
注：vst比rlog函数适用于大样本（如本例n=392）。rlog转换速度非常慢！大样本（n>30）不要尝试了！
#Extraction of transformed normalized expression matrix.
#Caution: rld is an object containing TRANSFORMED normalized matrix rather than normalized expression matrix. (Normalized expression values = original ones/size factor of each sample)
#If you would like to get the untransformed normalized matrix:
#expr_norm <- counts(dds_DE, normalized = TRUE)
rld <- vst(dds_DE, blind = FALSE)
expr_vst <- assay(rld)
#visualization of data normalization.
par(cex = 0.7)
#the meaning of the parameters: the lower/left/upper/right margins of the graph.
par(mar=c(1,3,1,1))
n.sample=ncol(expr_vst)
cols <- rainbow(n.sample*1.2)
par(mfrow = c(2, 1))
boxplot(expr_vst, col = cols,main="expression value",las=2)
而下图是未vst转换的标准化counts数据的可视化（随机选取部分样本）：
图片说明DESeq2的标准化+转换功力是非常强大的。
- 绘制主成分分析图（PCA plot）
主成分分析是一种线性降维的方法，将数据集中的基因按照表达特征归入到少数互不相关的主成分中，用主成分解释两组间差异。
输入标准化矩阵，使用DESeq2包中的plotPCA函数绘图。
plotPCA(rld, intgroup = 'sample_type')
图片说明两个最主要的主成分分别可解释组间32%和10%的差异，且normal和cancer样本可较好地根据主成分分开。
- 绘制火山图（volcano plot）。
火山图横轴为log2FC, 纵轴为校正后p值，可以直观反映各基因的数据分布状况。火山图有两种做法。
（1）基础plot作图：校正p值<0.01的基因表现为蓝色点，校正p值 < 0.01 & abs(log2FC) > 2表现为红色点。
par(mfrow=c(1,1,1,1))
# Make a basic volcano plot
with(res_DE, plot(log2FoldChange, -log10(padj), pch=20, main="Volcano plot", xlim=c(-3,3)))
# Add colored points: blue if padj<0.05, red if log2FC>1 and padj<0.05)
with(subset(res_DE, padj<0.05), points(log2FoldChange, -log10(padj), pch=20, col="blue"))
with(subset(res_DE, padj<0.05 & abs(log2FoldChange)>2), points(log2FoldChange, -log10(padj), pch=20, col="red"))
（2）ggplot2作图：校正p值< 0.05 & log2FC > 2为红色点，校正p值< 0.05 & log2FC < -2为蓝色点。
for_volcano <- data.frame('log2FoldChange'=res_DE$log2FoldChange,
'padj'=res_DE$padj,
'trend' = rep('no', length(res_DE$log2FoldChange)))
up_sig_indices <- intersect(which(for_volcano$log2FoldChange > 2), which(for_volcano$padj < 0.05))
down_sig_indices <- intersect(which(for_volcano$log2FoldChange < -2), which(for_volcano$padj < 0.05))
for_volcano[up_sig_indices,'trend'] <- 'up'
for_volcano[down_sig_indices, 'trend'] <- 'down'
for_volcano$trend <- as.factor(for_volcano$trend)
for_volcano$padj <- -log10(for_volcano$padj)
p <- ggplot(for_volcano, aes(x=log2FoldChange, y=padj, colour=trend))+
geom_point(size=I(0.7))+
scale_color_manual(values = c('no'='black', 'up'='red', 'down'='blue'))+
geom_vline(xintercept = c(2, -2), lty=2, size=I(0.4), colour = 'grey11')+
geom_hline(yintercept = c(-log(x=0.05, base=10)), lty=2, size=I(0.4), colour = 'grey11')+
theme_bw()+
theme(panel.background = element_rect(colour = 'black', size=1, fill = 'white'),
panel.grid = element_blank())+
labs(x='log2FoldChange', y='-log10Pvalue')
p
对比有点惨烈，看来ggplot2大行其道不是没有道理的。
- 绘制差异表达基因热图。
热图可以用每个色块的色阶表现样本之间某基因的差异表达情况。色阶代表了每个基因的标准化表达值在样本之间的z-score =(某基因中的表达值-所有样本基因表达均值)/所有样本基因表达值的标准差，反映了基因表达值在样本间的相对离散状态。
本例中，选取所有差异基因作热图，列为样本，行为差异基因，并根据基因表达特征对样本和基因进行聚类。按照第二级聚类（将基因或样本可分为两类）划分热图。由于大量样本导致基因表达离散性较大，z-score在极少部分样本基因中表现出极高或极低值（如绝对值>15），这会导致整体可视化颜色差异不显著。所以设置|z-score|对应色阶极值为2，超过该值均用2对应色阶进行可视化。
DEG_expr_matr <- expr_vst[resSigAll@rownames,]
#scale() will calculate the z-scores of each [column, that is gene] (that's why we use t function twice!).
z_score_matrix <- t(scale(t(DEG_expr_matr)))
annotation_col <- data.frame(Sample_type = condition_table$sample_type)
rownames(annotation_col) <- condition_table$TCGA_IDs
pheatmap(mat = z_score_matrix, cluster_rows = T, show_colnames = F, cluster_cols = T,
annotation_col = annotation_col,
annotation_colors = list(Sample_type=c(cancer = 'orange', normal = 'blue')),
show_rownames=F, color = colorRampPalette(c("navy", "white", "firebrick3"))(100),
legend = T, legend_breaks = c(-2, 0, 2), breaks = unique(seq(-2, 2, length=100)),
cutree_cols = 2, cutree_rows = 2,
lagend_labels = c('≤2','0','≥2'))
图中可看出差异基因筛选较为有效。此外对样本的聚类分析也较为准确，仅2个肿瘤样本落入癌旁组织区，说明筛选出的差异基因能较好地反映组间的差异。
- 可视化组间某单基因的read counts分布情况。
对某一特定基因感兴趣的朋友可以输入自己想查询的intended_gene获得组间的normalized read counts小提琴图（比箱线图更高大上一点的玩意儿，还可以反映样本分布密度）。这里查询一下一个EMT transcription factor TWIST1的组间表达情况。
TIPS：注意gene symbol需要和表达矩阵中的保持一致。如果关注的基因有多个别名，自己得额外关注一下~
（1）用vioplot包作图。
#Violin plot can display the density distribution better.
intended_gene <- 'TWIST1'
normal_gene <- expr_vst[intended_gene, condition_table$sample_type == 'normal']
cancer_gene <- expr_vst[intended_gene, condition_table$sample_type == 'cancer']
vioplot(normal_gene, cancer_gene, names = c('normal', 'cancer'), col = c('blue','orange'))
title(xlab = 'groups', ylab = 'normalized read counts')
（2）用ggplot2包作图。
intended_gene <- 'TWIST1'
normal_gene <- data.frame('values'=expr_vst[intended_gene, condition_table$sample_type == 'normal'], 'group'= rep('normal', sum(condition_table$sample_type == 'normal')))
cancer_gene <- data.frame('values'=expr_vst[intended_gene, condition_table$sample_type == 'cancer'], 'group'= rep('cancer', sum(condition_table$sample_type == 'cancer')))
for_violin <- rbind(normal_gene, cancer_gene)
#To reorder the rank of groups. Otherwise cancer group will be visualized first by alphabetical order.
for_violin$group <- as.factor(for_violin$group) %>% relevel(ref = 'normal')
p <- ggplot(for_violin, aes(x=group, y=values, fill=group))+
geom_violin(position = position_dodge(width = 1), scale = 'width')+
scale_fill_manual(values = c('normal'='blue', 'cancer'='orange'))+
geom_boxplot(position = position_dodge(width = 1), outlier.size = 0.7, width = 0.2, show.legend = FALSE) +
labs(x = 'groups', y = 'normalized read counts')+theme_bw()+
theme(panel.background = element_rect(colour = 'black', size=1, fill = 'white'),
panel.grid = element_blank())
p
这个两者倒没有太大差异，按喜好选用即可。
All right, 先写到这里吧，刚不住了！信息量应该足够慢慢看了~
