04-转录组下游分析-标准化、聚类、差异分析

准备工作

1.数据标准化

标准化前需要进行数据预处理

过滤低表达的基因，并检查是否有异常样本
以下是常见的几种过滤方式（过滤的标准都可以自己调整）
1：在至少在75%的样本中都表达的基因（表达是指在某个样本中count值＞0）
2：过滤平均值count<10的基因
3：过滤平均cpm <10 的基因

为什么做标准化？

计数结果的差异的影响因素：落在参考区域上下限的read是否需要被统 计，按照什么样的标准进行统计。
标准化的主要目的是去除测序数据的测序深度和基因长度。
• 测序深度：每个样本的测序深度（产生的数据量）不完全一样，同一条件下，测序深度越深，基因表达的read读数越多。同一条件下，测序深度越深，基因表达的read读数越多。
• 基因长度：同一条件下，不同的基因长度产生不对等的read读数，基因越长，该基因的read读数越高。视基因类型而定，比如microRNA稳定在20+长度，则不用针对长度进行标准化。

FPKM值

FPKM，Fragments per kilobase per Million mapped reads 片段数量per千碱基长度per百万比对上的reads
FPKM = cDNA Fragments/Mapped Fragments(Millions)×Transcript Length(kb)
• cDNA Fragments：比对到某一转录本上的片段数目，即双端Reads数目，针对双端数据；
• Mapped Fragments (Millions)：比对到转录本上的片段总数，以百万为单位，即10^6；
• Transcript Length(kb)：转录本长度，kb为单位，即10^3。

RPM值

RPM：Reads of exon model per Million mapped reads
RPM=total exon reads / mapped reads (Millions)
RPM方法：标准化了测序深度的影响，但没有考虑转录本的长度的影响。
RPM适合于产生的read读数不受基因长度影响的测序方法，比如miRNAseq测序，miRNA的长度一般在18~32个碱基之间。

TPM值

TPM，Transcripts Per Kilobase Million
TPMi=(Ni/Li)*1000000/sum(Ni/Li+………+ Nm/Lm)
和FPRM值一样都矫正了深度和长度，但TPM先对每个基因的read数用基因的长度进行校正，之后再用校正后的这个基因read数(Ni/Li)与校正后的这个样本的所有read数（sum(Ni/Li+………+ Nm/Lm)）求商，最后再乘以10^6。最后每个样本的TPM加和是相等的
Ni：mapping到基因i上的read数；
Li：基因i的外显子长度的总和

各种值的使用场合

1.差异表达分析：原始count值，算法输入要求（针对二代测序差异分析算法，算法内部一般有标化方法。）
2.标化后的值：基因表达值在样本与样本之间具有可比性。比如PCA分析，样本表达总体分布，生存分析，热图绘制，相关性分析等

实际代码

第一步：Step01-airwayCount.R
rm(list = ls())
library(stringr)
 
## ====================1.读取数据
# 读取raw count表达矩阵
rawcount <- read.table("data/raw_counts.txt",row.names = 1, 
                       sep = "\t", header = T)
colnames(rawcount)

# 查看表达谱
rawcount[1:4,1:4]
# 去除前的基因表达矩阵情况
dim(rawcount)
# 获取分组信息
# 差异分析方案为：Dex vs untreated
# 自己做一个分组信息的表格，第一列是样本，第二列是分组信息
group <- read.table("data/group.txt",
                    header = T,sep = "\t", quote = "\"")
group
## =================== 2.表达矩阵预处理
# 过滤低表达基因
# rawcount＞0判断后生成一个和rawcount表达矩阵对应的逻辑向量
# rawSums统计每个基因在每个样本中生成的逻辑向量的和，F为0，T为1，则为每个基因在多少个样本中有表达
# ncol指样本的数量，0.75*ncol(rawcount)指至少在75%的样本中有表达，再floor向下取整
keep <- rowSums(rawcount>0) >= floor(0.75*ncol(rawcount))
table(keep)
# 生成一个判断满足在75%的样本中有表达的逻辑向量矩阵
# 如果按平均值count<10，则需要改成rowMeans(rawcount)>10，更为严格
# 按keep逻辑向量矩阵取rawcount的子集，生成一个过滤后的矩阵filter_count
filter_count <- rawcount[keep,]
# 简单看一下确认是对的
filter_count[1:4,1:4]
dim(filter_count)
结果如下：
> keep <- rowSums(rawcount>0) >= floor(0.75*ncol(rawcount))
> table(keep)
keep
FALSE  TRUE 
53504 24794 
> 
> filter_count <- rawcount[keep,]
> filter_count[1:4,1:4]
                   AJV_34 AJV_35 AJV_49 JV_42
ENSMUSG00000102135     12      3      1     2
ENSMUSG00000100764      3      1      2     5
ENSMUSG00000100635      2      0      1     3
ENSMUSG00000100480      3      1      4    10
> dim(filter_count)
[1] 24794     6
## 数据预处理,如果采用2.过滤平均值count<10的基因
> keep <- rowMeans(rawcount)>=10
> table(keep)
keep
FALSE  TRUE 
62689 15609

# 此处选择使用cpm值进行标准化
# 加载edgeR包计算counts per millio(cpm) 表达矩阵
library(edgeR)
express_cpm <- cpm(filter_count)
express_cpm[1:6,1:6]

# 保存表达矩阵和分组结果
save(filter_count, express_cpm, group, file = "data/Step01-airwayData.Rdata")

异常样本与重复性检测

why需要这一步？

1.检查组内样本的重复性是否好
2.检查组间样本是否可以较好的分开

步骤

1.样本表达水平总体分布图
2.样本相关性检验
3.PCA分析
4.层次聚类分析

样本表达总体分布图

使用箱线图/小提琴图/概率密度分布图，观察：
1.单个样本的数据分布特点
2.多个样本之间的差异

1：画箱线图

箱式图：主要是从四分位数的角度来描述数据的分布
从箱线图中不仅可以查看单个样品表达水平分布的离散程度，还可以直观地比较不同样品整体表 达水平。

第二步：Step02-sampleDistribution.R
rm(list = ls())
options(stringsAsFactors = F)

# 加载包，设置绘图参数
library(ggplot2)
library(ggsci)
library(tidyverse)

# 设置一个喜欢的主题
mythe <- theme_bw() + theme(panel.grid.major=element_blank(),
                            panel.grid.minor=element_blank())

# 加载原始表达的数据
lname <- load(file = "data/Step01-airwayData.Rdata")
lname

## [1] "filter_count" "express_cpm"  "group"

# 使用标准化后的数据进行画箱式图，cpm值的矩阵取1og10后+1，因为log的底数不能为0，生成用于比较画图的表达矩阵
exprSet <- log10(as.matrix(express_cpm)+1)
exprSet[1:6,1:6]

## 1.样本表达总体分布-箱式图
# 构造绘图数据：宽数据转为长数据
# 把所有数据都放到一列中，列名的值的新列名为“sample”，值的新列名为“expression”，基因编号丢弃
#cols = everything()，把所有的列转换为行，names_to = "sample"，转换为行的时候的列名，values_to = "expression")，转换为行后值得列名

data <- exprSet %>% 
  as.data.frame() %>% 
  pivot_longer(cols = everything(), names_to = "sample",values_to = "expression")

head(data)

p <- ggplot(data = data, aes(x=sample,y=expression,fill=sample))
p1 <- p + geom_boxplot() + 
  mythe+ theme(axis.text.x = element_text(angle = 90)) + 
  xlab(NULL) + ylab("log10(CPM+1)") + scale_fill_lancet()

p1
# 保存图片
png(file = "result/1.sample_boxplot.png",width = 800, height = 900,res=150)
print(p1)
dev.off()

2.样本表达总体分布-小提琴图

简要展示分布“密度”，突出数据分布的密集 区域
从小提琴图中可以查看单个样品表达水平分 布的密集程度，也可以直观地比较不同样品整体表达水平。
宽—>密度高

## 2.样本表达总体分布-小提琴图
p2 <- p + geom_violin() +  mythe +
  theme(axis.text = element_text(size = 12),
        axis.text.x = element_text(angle = 90)) + 
  xlab(NULL) + ylab("log10(CPM+1)")+scale_fill_lancet()
p2

# 保存图片
png(file = "result/1.sample_violin.png",width = 800, height = 900,res=150)
print(p2)
dev.off()

密度曲线图

从概率密度的角度描述基因表达量总体分 布图，能反映样品中基因的整体表达模式
图中不同颜色的曲线代表不同的样品，横 坐标表示对应样品 log2(cpm+1)的对数值， 纵坐标表示概率密度。

## 3.样本表达总体分布-概率密度分布图
m <- ggplot(data=data, aes(x=expression)) 
p3 <- m +  geom_density(aes(fill=sample, colour=sample),alpha = 0.1) + 
  xlab("log10(CPM+1)") + mythe +scale_fill_npg()
p3

# 保存图片
png(file = "result/1.sample_density.png",width = 800, height = 700, res=150)
print(p3)
dev.off()

层次聚类

层次聚类树

通过计算点与点之间的空间距离对样本进行类别划分

# 魔幻操作，一键清空
rm(list = ls()) 
options(stringsAsFactors = F)

library(FactoMineR)
library(factoextra)
library(corrplot)
library(pheatmap)
library(tidyverse)
 
# 加载数据并检查
lname <- load(file = 'data/Step01-airwayData.Rdata')
lname

## 1.样本之间的相关性-层次聚类树----
dat <- log10(express_cpm+1)
dat[1:4,1:4]
dim(dat)
sampleTree <- hclust(dist(t(dat)), method = "average")
plot(sampleTree)

PCA图

通过提取样本的综合特征，即主成分（第一主成分，第二主成分…）来对样本进行分类

## 2.样本之间的相关性-PCA----
# 第一步，数据预处理
dat <- log10(express_cpm+1)
dat[1:4,1:4]

dat <- as.data.frame(t(dat))
dat_pca <- PCA(dat, graph = FALSE)

group_list <- group[match(group$run_accession,rownames(dat)), 2]
group_list

# geom.ind: point显示点，text显示文字
# palette: 用不同颜色表示分组
# addEllipses: 是否圈起来
mythe <- theme_bw() + 
  theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank()) +
  theme(plot.title = element_text(hjust = 0.5))

p <- fviz_pca_ind(dat_pca,
                  geom.ind = "point",  #text
                  col.ind = group_list, 
                  palette = c("#00AFBB", "#E7B800"), 
                  addEllipses = T,  
                  legend.title = "Groups") + mythe
p

# 保存结果
pdf(file = "result/2.sample_PCA.pdf",width = 6.5,height = 6)
plot(p)
dev.off()

相关性分析

通过计算样本与样本之间的相关性系数来对样本进行分类，相关性系数可以是pearson, spearman, kendall，不同的方法要求的数据预处理有所不同，如pearson要取log

## 3.样本之间的相关性-cor----
# 选择差异变化大的基因算样本相关性
exprSet <- express_cpm
exprSet = exprSet[names(sort(apply(exprSet, 1, mad),decreasing = T)[1:800]),]
dim(exprSet)

# 计算相关性
M <- cor(exprSet,method = "spearman")
M

# 构造注释条
anno <- data.frame(group=group$sample_title,row.names = group$run_accession )

# 保存结果
pheatmap(M,display_numbers = T,
         annotation_col = anno,
         fontsize = 10,cellheight = 30,
         cellwidth = 30,cluster_rows = T,
         cluster_cols = T,
         filename = "result/2.sample_Cor.pdf",width = 7.5,height = 7)

差异表达分析

步骤：

1,创建设计矩阵和对比
2.构建edgeR的DGEList对象，并归一化，拟合模型
3.提取分析结果并筛选显著差异基因

rm(list = ls())
options(stringsAsFactors = F)

# 加载包
library(edgeR)
library(ggplot2)
 
# 读取基因表达矩阵信息并查看分组信息和表达矩阵数据
lname <- load(file = "data/Step01-airwayData.Rdata")
lname

# 表达谱
filter_count[1:4,1:4]

# 分组信息
group_list <- group[match(colnames(filter_count),group$run_accession),2]
group_list

# treat vs control
comp <- unlist(strsplit("A1KO_vs_WT",split = "_vs_"))
group_list <- factor(group_list,levels = comp)
group_list
table(group_list)


# 构建线性模型。0代表x线性模型的截距为0
design <- model.matrix(~0+group_list)
rownames(design) <- colnames(filter_count)
colnames(design) <- levels(factor(group_list))
design

# 构建edgeR的DGEList对象
DEG <- DGEList(counts=filter_count, 
               group=factor(group_list))

# 归一化基因表达分布
DEG <- calcNormFactors(DEG)

# 计算线性模型的参数
DEG <- estimateGLMCommonDisp(DEG,design)
DEG <- estimateGLMTrendedDisp(DEG, design)
DEG <- estimateGLMTagwiseDisp(DEG, design)

# 拟合线性模型
fit <- glmFit(DEG, design)

# 进行差异分析
lrt <- glmLRT(fit, contrast=c(1,-1)) 

# 提取过滤差异分析结果
DEG_edgeR <- as.data.frame(topTags(lrt, n=nrow(DEG),adjust.method = "BH"))
head(DEG_edgeR)

# 筛选上下调，设定阈值
fc_cutoff <- 1.5
pvalue <- 0.05

DEG_edgeR$regulated <- "normal"

loc_up <- intersect(which( DEG_edgeR$logFC > log2(fc_cutoff) ),
                    which( DEG_edgeR$PValue < pvalue) )

loc_down <- intersect(which(DEG_edgeR$logFC < (-log2(fc_cutoff))),
                      which(DEG_edgeR$PValue<pvalue))

DEG_edgeR$regulated[loc_up] <- "up"
DEG_edgeR$regulated[loc_down] <- "down"

table(DEG_edgeR$regulated)

## 添加一列gene symbol
# 方法1：使用包
library(org.Mm.eg.db)
#library(org.Hs.eg.db)
keytypes(org.Mm.eg.db)
#keytypes(org.Hs.eg.db)

library(clusterProfiler)
id2symbol <- bitr(rownames(DEG_edgeR), 
                  fromType = "ENSEMBL", 
                  toType = "SYMBOL", 
                  OrgDb = org.Mm.eg.db)
head(id2symbol)

DEG_edgeR <- cbind(GeneID=rownames(DEG_edgeR),DEG_edgeR)
DEG_edgeR_symbol <- merge(id2symbol,DEG_edgeR,
                          by.x="ENSEMBL",by.y="GeneID",all.y=T)
head(DEG_edgeR_symbol)


# 方法2：gtf文件中得到的id与name关系
# Assembly: GRCh37(hg19) Release: ？
# 使用上课测试得到的count做

# 选择显著差异表达的结果
library(tidyverse)
DEG_edgeR_symbol_Sig <- filter(DEG_edgeR_symbol,regulated!="normal")

# 保存
write.csv(DEG_edgeR_symbol,"result/4.DEG_edgeR_all.csv", row.names = F)
write.csv(DEG_edgeR_symbol_Sig,"result/4.DEG_edgeR_Sig.csv", row.names = F)
save(DEG_edgeR_symbol,file = "data/Step03-edgeR_nrDEG.Rdata")

##====== 检查是否上下调设置错了
# 挑选一个差异表达基因
head(DEG_edgeR_symbol_Sig)

exp <- c(t(express_cpm[match("ENSMUSG00000000001",rownames(express_cpm)),]))
test <- data.frame(value=exp, group=group_list)
dev.off()
ggplot(data=test,aes(x=group,y=value,fill=group)) + geom_boxplot()

如果想在服务器上使用RStudio怎么操作？

1.如果服务器上安装了RStudio，可以用在线的RStudio Server，和本地一样使用
2.传参脚本
Rscript edgeR.R --count filter_count.txt --group group.txt --comp Dex_vs_untreated --fc 1.5 –
pvalue 0.05 --od ./

nohup sh DEG.sh >DEG.log &

火山图(Volcano Plot)

通过火山图可查看基因 在两个（组）样品中表达水平的差异 即FC值，以及差异的统计学显著性即FDR。
差异表达火山图中的每一个点表示一个基因, 横坐标表示某一 个基因在两样品中表达量差异倍数的对数值;纵坐标表示FDR 的负对数值。横坐标绝对值越大，说明表达量在两样品间的 表达量倍数差异越大;纵坐标值越大，表明差异表达越显著， 筛选得到的差异表达基因越可靠。

rm(list = ls())
options(stringsAsFactors = F)
library(ggplot2)
library(tidyverse)
 
# 读差异分析结果
lname <- load(file = "data//Step03-edgeR_nrDEG.Rdata")

# 根据需要修改DEG的值
data <- DEG_edgeR_symbol
colnames(data)


# 绘制火山图
colnames(data)
p <- ggplot(data=data, aes(x=logFC, y=-log10(PValue),color=regulated)) + 
  geom_point(alpha=0.5, size=1.2) + 
  theme_set(theme_set(theme_bw(base_size=20))) + theme_bw() +
  theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank()) +
  xlab("log2FC") + ylab("-log10(Pvalue)") +
  scale_colour_manual(values = c(down='blue',normal='grey',up='red')) +
  geom_vline(xintercept=c(-(log2(1.5)),log2(1.5)),lty=2,col="black",lwd=0.6) +
  geom_hline(yintercept = -log10(0.05),lty=2,col="black",lwd=0.6)
p


# 添加top基因
# 通过FC选取TOP10
label <- data[order(abs(data$logFC),decreasing = T)[1:10],]
# 通过pvalue选取TOP10
#label <- data[order(abs(data$PValue),decreasing = F)[1:10],]
label <- na.omit(label)

p1 <- p + geom_point(size = 3, shape = 1, data = label) +
  ggrepel::geom_text_repel( aes(label = SYMBOL), data = label, color="black" )

p1


# 保存结果
png(file = "result/5.Volcano_Plot.png",width = 900, height = 800, res=150)
plot(p1)
dev.off()

差异分析热图

不仅对差异表达基因进行了可视化， 还对筛选出的差异表达基因做层次聚 类分析，将具有相同或相似表达行为 的基因进行聚类

rm(list = ls())
options(stringsAsFactors = F)
# 加载包
library(pheatmap)
library(tidyverse)
 
# 加载原始表达矩阵
lname <- load(file = "data/Step01-airwayData.Rdata")
lname

express_cpm1 <- rownames_to_column(as.data.frame(express_cpm) ,var = "ID")

# 读取差异分析结果
lname <- load(file = "data/Step03-edgeR_nrDEG.Rdata")
lname

# 提取所有差异表达的基因名
edgeR_sigGene <- DEG_edgeR_symbol[DEG_edgeR_symbol$regulated!="normal",]
head(edgeR_sigGene)

data <- merge(edgeR_sigGene,express_cpm1,by.x = "ENSEMBL",by.y = "ID")
data <- na.omit(data)
data <- data[!duplicated(data$SYMBOL),]

# 绘制热图
dat <- data[, grep("^(AJV|JV)", names(data), value = TRUE)]
#dat <- select(data,starts_with("JV"))
rownames(dat) <- data$SYMBOL
dat[1:4,1:4]
anno <- data.frame(group=group$sample_title,row.names = group$run_accession)

pheatmap(dat,scale = "row",show_colnames =T,
         show_rownames = F, cluster_cols = T,
         annotation_col=anno,
         main = "edgeR's DEG")


# 显示指定symbol,这里随便展示10个基因symbol
labels <- rep(x = "",times=nrow(dat))
labels[1:10] <- rownames(dat)[1:10]
pheatmap(dat,scale = "row",show_colnames =T,
         show_rownames = T, 
         cluster_cols = T,
         annotation_col=anno,
         labels_row = labels,
         fontsize_row = 8,
         main = "edgeR's DEG")

# 按照指定顺序绘制热图
dex_exp <- express_cpm[,match(rownames(anno)[which(anno$group=="Dex")],
                              colnames(express_cpm))]

untreated_exp <- express_cpm[,match(rownames(anno)[which(anno$group=="WT")],
                              colnames(express_cpm))]

data_new <- cbind(dex_exp, untreated_exp)
dat1 <- data_new[match(edgeR_sigGene$ENSEMBL,rownames(data_new)),]

pheatmap(dat1, scale = "row",show_colnames =T,show_rownames = F, 
              cluster_cols = F, 
              annotation_col=anno,
              main = "edgeR's DEG")