一、数据标准化
micro rna测序常用
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)
# 获取分组信息
group <- read.table("data/filereport_read_run_PRJNA229998_tsv.txt",
header = T,sep = "t", quote = """)
colnames(group)
# 提取表达矩阵对应的样本表型信息
group <- group[match(colnames(rawcount), group$run_accession),
c("run_accession","sample_title")]
group
# 差异分析方案为：Dex vs untreated
group$sample_title <- str_split_fixed(group$sample_title,pattern = "_", n=2)[,2]
group
write.table(group,file = "data/group.txt",row.names = F,sep = "t",quote = F)
## =================== 2.表达矩阵预处理
# 过滤低表达基因
keep <- rowSums(rawcount>0) >= floor(0.75*ncol(rawcount))
table(keep)
filter_count <- rawcount[keep,]
filter_count[1:4,1:4]
dim(filter_count)
# 加载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")
二、异常样本与重复性检测
1.样本总体分布
rm(list = ls())
options(stringsAsFactors = F)
# 加载包，设置绘图参数
library(ggplot2)
library(ggsci) #绘图的配色包
mythe <- theme_bw() theme(panel.grid.major=element_blank(),
panel.grid.minor=element_blank()) #主题
# 加载原始表达的数据
lname <- load(file = "data/Step01-airwayData.Rdata")
lname
exprSet <- log10(as.matrix(express_cpm) 1)
exprSet[1:6,1:6]
## 1.样本表达总体分布-箱式图
# 构造绘图数据
data <- data.frame(expression=c(exprSet),
sample=rep(colnames(exprSet),each=nrow(exprSet)))
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.样本表达总体分布-小提琴图
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()
## 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()
2.样本之间的相关性
1.层次聚类树
2.PCA主成分分析
3.相关性分析
pearson:对离异值非常敏感，如果有一个值与正常值差很远会导致数据相关性很低，所以通常进行log处理之后再进行pearson分析。 pearson计算出的相关性比spearman计算出的相关性要高。
spearman：对异常值不敏感，log转换不改变相关大小。
代码语言：javascript复制# 魔幻操作，一键清空
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) #log 1是因为log之后要求底数不为0，故而 1
dat[1:4,1:4]
dim(dat)
sampleTree <- hclust(dist(t(dat)), method = "average")
plot(sampleTree)
# 提取样本聚类信息
temp <- as.data.frame(cutree(sampleTree,k = 2)) %>%
rownames_to_column(var="sample")
temp1 <- merge(temp,group,by.x = "sample",by.y="run_accession")
table(temp1$`cutree(sampleTree, k = 2)`,temp1$sample_title)
# 保存结果
pdf(file = "result/2.sample_Treeplot.pdf",width = 7,height = 6)
plot(sampleTree)
dev.off()
## 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", #point
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()
## 3.样本之间的相关性-cor----
# 选择差异变化大的基因算样本相关性
exprSet <- express_cpm
exprSet = exprSet[names(sort(apply(exprSet, 1, mad),decreasing = T)[1:800]),]
dim(exprSet)
# 计算相关性
M <- cor(log10(exprSet 1),method = "pearson")
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,
width = 7.5,height = 7)
三、差异表达分析
1.edge 差异分析
p value 看显著程度 FDR：校正后的p值
logFC看差异程度 fold change，取log之后通过正负号来判断上调和下调
代码语言：javascript复制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("Dex_vs_untreated",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')) #这里adjust method是用于校正p值的方法
head(DEG_edgeR)
# 筛选上下调，设定阈值
fc_cutoff <- 1.5 #FC值的阈值
fdr <- 0.05
DEG_edgeR$regulated <- "normal" #加一列表达情况
loc_up <- intersect(which( DEG_edgeR$logFC > log2(fc_cutoff) ),
which( DEG_edgeR$PValue < fdr) ) #得到表达上调的位置坐标
loc_down <- intersect(which(DEG_edgeR$logFC < (-log2(fc_cutoff))),
which(DEG_edgeR$PValue< fdr)) #得到表达下调的位置坐标
DEG_edgeR$regulated[loc_up] <- "up"
DEG_edgeR$regulated[loc_down] <- "down"
table(DEG_edgeR$regulated)
## 添加一列gene symbol
# 方法1：使用包
library(org.Hs.eg.db)
keytypes(org.Hs.eg.db)
library(clusterProfiler)
id2symbol <- bitr(rownames(DEG_edgeR),
fromType = "ENSEMBL",
toType = "SYMBOL",
OrgDb = org.Hs.eg.db) #有23.08%转换失败，因为GRCH版本不一样
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("ENSG00000001626",rownames(express_cpm)),]))
test <- data.frame(value=exp, group=group_list)
ggplot(data=test,aes(x=group,y=value,fill=group)) geom_boxplot()
2.传参脚本
在diff里面 edgeR.R 中
需要通过Fillzilla传递传参脚本到服务器上
四、差异结果可视化
1.热图
代码语言：javascript复制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 <- select(data,starts_with("SRR"))
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") #label_rows 可以用来展示感兴趣的基因
# 按照指定顺序绘制热图
dex_exp <- express_cpm[,match(rownames(anno)[which(anno$group=="Dex")],
colnames(express_cpm))]
untreated_exp <- express_cpm[,match(rownames(anno)[which(anno$group=="untreated")],
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") #设置treeheight_row = 0，就可以隐去聚类树2.火山图
代码语言：javascript复制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()
