測序上游分析系列：
mRNA-seq轉(zhuǎn)錄組二代測序從raw reads到表達(dá)矩陣：上中游分析pipeline
miRNA-seq小RNA高通量測序pipeline：從raw reads州弟，鑒定已知miRNA-預(yù)測新miRNA，到表達(dá)矩陣【一】
miRNA-seq小RNA高通量測序pipeline：從raw reads讯榕，鑒定已知miRNA-預(yù)測新miRNA农猬，到表達(dá)矩陣【二】
其他文章系列：ggplot2作圖篇：
（1）基于ggplot2的RNA-seq轉(zhuǎn)錄組可視化：總述和分文目錄
（2）測序結(jié)果概覽：基因表達(dá)量rank瀑布圖拟枚，高密度表達(dá)相關(guān)性散點(diǎn)圖，注釋特定基因及errorbar的表達(dá)相關(guān)性散點(diǎn)圖繪制
（3）單/多個(gè)基因在組間同時(shí)展示的多種選擇：分組小提琴圖、分組/分面柱狀圖查近、單基因蜂群點(diǎn)圖拼圖的繪制
（4）配對(duì)樣本基因表達(dá)趨勢：優(yōu)化前后的散點(diǎn)連線圖+拼圖繪制
上次說到了TCGA高通量轉(zhuǎn)錄組測序數(shù)據(jù)的下載锨咙，整合與臨床數(shù)據(jù)的清洗语卤。接下來的必要流程必定和肯定是癌組織與癌旁組織的差異表達(dá)分析了。
R語言中有兩個(gè)封裝R包廣泛用于差異基因的篩選：DESeq2和EdgeR酪刀。由于測序數(shù)據(jù)屬于負(fù)二項(xiàng)分布粹舵，而且不同測序批次間總體深度會(huì)有區(qū)別，因此他們不能直接使用t檢驗(yàn)進(jìn)行差異分析骂倘！封裝包幫我們解決了這個(gè)問題：本教程使用DESeq2進(jìn)行篩選眼滤，不需要知道原理，操作過程非常簡潔历涝，就可以得到結(jié)果诅需。
本例仍使用之前教程選擇的的TCGA-LUSC白色人種肺鱗癌表達(dá)譜數(shù)據(jù)。下載整合操作見前文：從TCGA數(shù)據(jù)庫下載并整合清洗高通量腫瘤表達(dá)譜-臨床性狀數(shù)據(jù)
絕對(duì)的TIPS：DESeq2以及EdgeR包要求的輸入文件為【未經(jīng)標(biāo)準(zhǔn)化的原始read counts表達(dá)矩陣】荧库！DESeq2會(huì)對(duì)原始counts進(jìn)行標(biāo)準(zhǔn)化（該方法優(yōu)于CPM/TPM/FPKM/RPKM等方法）堰塌，去除不同測序批次導(dǎo)致的批次效應(yīng)（batch effect），因此不能使用FPKM等已經(jīng)標(biāo)準(zhǔn)化過的矩陣分衫。
所以在TCGA-portal中需要下載的文件是HTSeq-counts.txt场刑。另外，上篇教程提到過的丐箩，counts文件在鏡像操作的同時(shí)摇邦，需要在獲得最初的整合count matrix時(shí)去除末尾5行總結(jié)性（mei）注（yong）釋數(shù)據(jù)。
在經(jīng)過整合篩選后屎勘，我們從之前的教程中獲得（或者說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.
注：從上個(gè)教程到這里施籍，別忘了先把mRNA_exprSet的第一列'gene_name'并入rownames！否則后續(xù)分析columns數(shù)對(duì)不上概漱。
rownames(mRNA_exprSet) <- mRNA_exprSet$gene_name
mRNA_exprSet <- mRNA_exprSet[,-1]
原本392個(gè)樣本保留386個(gè)丑慎，其中包含344個(gè)癌組織樣本，以及42個(gè)癌旁組織樣本。
分析部分
- 下載并加載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的分組信息進(jìn)行因子化进副，并設(shè)置reference。
DESeq2要求condition_table的格式為：至少含一列樣本的分組信息悔常，且rownames需要是可以與表達(dá)矩陣colnames匹配的樣本名稱影斑。此外grouping information需要是factor格式。
我們有必要讓R知道我們設(shè)置組別的意義机打。如果不手動(dòng)設(shè)置分組的reference（差異分析中的分母淮逻，也就是cancer與normal比較中的normal）毕箍，DESeq2默認(rèn)按照分組的首字母順序設(shè)置reference抽高，顯然是不準(zhǔn)確的萨赁。
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較低的基因。
如果某基因在樣本間的整體表達(dá)水平很低芥挣，則小幅度的read counts擾動(dòng)即可造成較大的fold change差異驱闷，但這種差異屬于noise，并不具有生物學(xué)意義九秀。所以需要設(shè)置一定的threshold將這些基因去除遗嗽，減少假陽性。
過濾的原則：原理類似鼓蜒，具體操作各自不同痹换。有分析使用如下原則：若某基因在所有樣本中的read counts之和小于樣本總數(shù)（樣本均read counts<1）則去除。
本例使用相對(duì)更苛刻的原則：若某基因在大于80%的樣本中read counts數(shù)<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,]
- 將已經(jīng)調(diào)整好的輸入文件輸入DESeq2，一步自動(dòng)分析畅厢。
封裝函數(shù)幫我們一步完成了以下步驟：
（1）estimation of size factor of each sample for normalization.
（2）estimation of dispersion.
（3）negative binomial GLM fitting and Wald tests.
由于TCGA樣本量比較大冯痢，所以這一步會(huì)需要較長時(shí)間。作為參考框杜，macbook Rstudio跑392個(gè)樣本用了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)
- 以表格的形式，按照自己的標(biāo)準(zhǔn)分別輸出所有咪辱、上調(diào)和下調(diào)的差異基因振劳。
dds_DE對(duì)象中包含了差異分析的所有結(jié)果。我們需要自己設(shè)定cut-off值進(jìn)行輸出：
（1）首先在result函數(shù)中設(shè)定alpha第一類錯(cuò)誤率（即p值threshold）, 以及再次申明比較類型（contrast參數(shù)設(shè)置如下：c(condition_table中分組變量名油狂，分子組名历恐，分母組名)）并按校正后p值（padj寸癌，p.adjusted, q-value, 也是False Discovery Rate, FDR, 通常可以根據(jù)差異基因數(shù)量設(shè)置0.1弱贼，0.05或0.01）大小排序蒸苇。
（2）subset函數(shù)設(shè)定校正后p值，以及l(fā)og2FoldChange閾值（相對(duì)表達(dá)差異的絕對(duì)值超過該閾值視為差異表達(dá)）吮旅。
本例中樣本較多溪烤，為控制差異基因數(shù)量，使用padj < 0.05和|log2FoldChange| > 2（FoldChange大于4）作為cut-off值庇勃。
TIPS：測序數(shù)據(jù)差異分析需要用校正p值FDR而不是p值作為cut-off氛什！用于控制假陽性率。
注：最后本地保存差異基因表我選擇以txt格式保存匪凉！因?yàn)閏sv等格式常默認(rèn)用excel打開，隨后某些基因名會(huì)被Excel自動(dòng)轉(zhuǎn)換為年月...血淚教訓(xùn)稗嗝省再层！
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')
所以我們就獲得了代表所有/上調(diào)/下調(diào)差異基因的三個(gè)csv表格！我們只需要關(guān)注log2FoldChange, padj和表達(dá)趨勢即可堡纬。
可視化部分
- 繪制各樣本標(biāo)準(zhǔn)化后read counts的箱線圖聂受。
用于可視化表達(dá)數(shù)據(jù)標(biāo)準(zhǔn)化的程度。若各樣本normalized read counts的分布水平相對(duì)穩(wěn)定烤镐，則可以認(rèn)為標(biāo)準(zhǔn)化的操作以及后續(xù)的分析是可信有效的蛋济。
使用vst函數(shù)對(duì)數(shù)據(jù)進(jìn)行variance stablizing transformation，并獲得轉(zhuǎn)換后的表達(dá)矩陣炮叶。vst轉(zhuǎn)換可以使表達(dá)矩陣適用于可視化碗旅、聚類以及機(jī)器學(xué)習(xí)等操作分析。
注：vst比rlog函數(shù)適用于大樣本（如本例n=392）镜悉。rlog轉(zhuǎn)換速度非常慢祟辟！大樣本（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轉(zhuǎn)換的標(biāo)準(zhǔn)化counts數(shù)據(jù)的可視化（隨機(jī)選取部分樣本）：
圖片說明DESeq2的標(biāo)準(zhǔn)化+轉(zhuǎn)換功力是非常強(qiáng)大的侣肄。
- 繪制主成分分析圖（PCA plot）
主成分分析是一種線性降維的方法旧困，將數(shù)據(jù)集中的基因按照表達(dá)特征歸入到少數(shù)互不相關(guān)的主成分中，用主成分解釋兩組間差異稼锅。
輸入標(biāo)準(zhǔn)化矩陣吼具，使用DESeq2包中的plotPCA函數(shù)繪圖。
plotPCA(rld, intgroup = 'sample_type')
圖片說明兩個(gè)最主要的主成分分別可解釋組間32%和10%的差異矩距，且normal和cancer樣本可較好地根據(jù)主成分分開拗盒。
- 繪制火山圖（volcano plot）。
火山圖橫軸為log2FC, 縱軸為校正后p值剩晴，可以直觀反映各基因的數(shù)據(jù)分布狀況锣咒∏肿矗火山圖有兩種做法。
（1）基礎(chǔ)plot作圖：校正p值<0.01的基因表現(xiàn)為藍(lán)色點(diǎn)毅整，校正p值 < 0.01 & abs(log2FC) > 2表現(xiàn)為紅色點(diǎn)趣兄。
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為紅色點(diǎn)，校正p值< 0.05 & log2FC < -2為藍(lán)色點(diǎn)悼嫉。
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
對(duì)比有點(diǎn)慘烈艇潭，看來ggplot2大行其道不是沒有道理的。
- 繪制差異表達(dá)基因熱圖戏蔑。
熱圖可以用每個(gè)色塊的色階表現(xiàn)樣本之間某基因的差異表達(dá)情況蹋凝。色階代表了每個(gè)基因的標(biāo)準(zhǔn)化表達(dá)值在樣本之間的z-score =(某基因中的表達(dá)值-所有樣本基因表達(dá)均值)/所有樣本基因表達(dá)值的標(biāo)準(zhǔn)差，反映了基因表達(dá)值在樣本間的相對(duì)離散狀態(tài)总棵。
本例中鳍寂，選取所有差異基因作熱圖，列為樣本情龄，行為差異基因迄汛，并根據(jù)基因表達(dá)特征對(duì)樣本和基因進(jìn)行聚類。按照第二級(jí)聚類（將基因或樣本可分為兩類）劃分熱圖骤视。由于大量樣本導(dǎo)致基因表達(dá)離散性較大鞍爱，z-score在極少部分樣本基因中表現(xiàn)出極高或極低值（如絕對(duì)值>15），這會(huì)導(dǎo)致整體可視化顏色差異不顯著专酗。所以設(shè)置|z-score|對(duì)應(yīng)色階極值為2睹逃，超過該值均用2對(duì)應(yīng)色階進(jìn)行可視化。
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'))
圖中可看出差異基因篩選較為有效祷肯。此外對(duì)樣本的聚類分析也較為準(zhǔn)確沉填，僅2個(gè)腫瘤樣本落入癌旁組織區(qū)，說明篩選出的差異基因能較好地反映組間的差異佑笋。
- 可視化組間某單基因的read counts分布情況拜轨。
對(duì)某一特定基因感興趣的朋友可以輸入自己想查詢的intended_gene獲得組間的normalized read counts小提琴圖（比箱線圖更高大上一點(diǎn)的玩意兒，還可以反映樣本分布密度）允青。這里查詢一下一個(gè)EMT transcription factor TWIST1的組間表達(dá)情況橄碾。
TIPS：注意gene symbol需要和表達(dá)矩陣中的保持一致。如果關(guān)注的基因有多個(gè)別名颠锉，自己得額外關(guān)注一下~
（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
這個(gè)兩者倒沒有太大差異琼掠，按喜好選用即可拒垃。
All right, 先寫到這里吧，剛不住了瓷蛙！信息量應(yīng)該足夠慢慢看了~
