从Fastq文件开始的转录组上下游分析 -下

Mirai233

已于 2024-09-01 19:56:34 修改

阅读量447

点赞数 4

文章标签：云计算 linux

于 2024-09-01 19:54:08 首次发布

本文链接：https://blog.csdn.net/Mirai233/article/details/141601534

版权

两篇重要文章，借鉴思路：

你确定你的差异基因找对了吗？

差异分析的时候到底是p值重要还是变化倍数重要？

RNA-seq入门实战（四）：差异分析前的准备——数据检查 - 简书

以下展示了样本hclust 图、距离热图、PCA图、前500差异性大的基因热图、相关性热图（选取了500高表达基因，防止低表达基因造成的干扰），确定我们不同样本间确实是有差异的。这些图并不是全都是必须的，它们全都是为了说明一个问题：我们的不同分组间确实存在差异。

如果这三张图出现的结果是各个处理组之间，能聚集起来，和其他组分开。

再进行PCA分组等之前，先使用R语言合并上一步中features输出的txt文件。以下使用R语言。

1.合并所有文件

#更改一下当前的文件路径设置
setwd("D:/R/6h") 
getwd()   #查看文件路径
#加载数据，CC900011*2, control*3, GSKLSD1*3, LSD1IN20*2, OGD*3, SP2509*2, SP2577*2
options(stringsAsFactors = FALSE)

control1<-read.table("count/6h_control_1.txt",sep = "\t",header =T)
#这个操作是把下面的数据集第七列的名称改成control1/实验组等,是所有文件哦
colnames(control1)[7]<-("control1")

control2<-read.table("count/6h_control_2.txt",sep = "\t",header =T)
colnames(control2)[7]<-("control2")

control3<-read.table("count/6h_control_3.txt",sep = "\t",header =T)
colnames(control3)[7]<-("control3")

OGD1<-read.table("count/6h_OGD_1.txt",sep = "\t",header =T)
colnames(OGD1)[7]<-("OGD1")

OGD2<-read.table("count/6h_OGD_2.txt",sep = "\t",header =T)
colnames(OGD2)[7]<-("OGD2")

OGD3<-read.table("count/6h_OGD_3.txt",sep = "\t",header =T)
colnames(OGD3)[7]<-("OGD3")

………………




#删除一些不太重要的信息比如基因所在的位置和长度等信息，即数据集第2-6列，第一列是GNEE ID,第7列是样品名，比如control2，或者treat3
control1 <- control1[,-(2:6)]
control2 <- control2[,-(2:6)]
control3 <- control3[,-(2:6)]
OGD1 <- OGD1[,-(2:6)]
OGD2 <- OGD2[,-(2:6)]
OGD3 <- OGD3[,-(2:6)]
……………………

#合并上述6个数据集，但是merge这个函数只可以合并两个数据集，所以采用了下述代码
raw_count1<- merge(control1, control2, by="Geneid")
raw_count1<-merge(raw_count1, control3, by="Geneid")

raw_count2<-merge(OGD1, OGD2, by="Geneid")
raw_count2 <- merge(raw_count2, OGD3, by="Geneid")

raw_count  <- merge(raw_count1, raw_count2, by="Geneid")

……………………
#由于数据库里的基因ID没有小数点，所以需要进一步替换为整数的形式。
#将匹配到的.以及后面的数字连续匹配并替换为空，并赋值给ENSEMBL
ENSEMBL <- gsub("\\.\\d*", "", raw_count$Geneid) 
#参考链接里是把 ENSEMBL 作为列名加入数据集，这儿做了改动，使ENSEMBL作为一列的内容插入到数据集第2列中，以便后续需要时进行更改，再把数据集第一列，也就是带有小数点之后数字的一列去掉。
raw_count <- cbind(raw_count[,1],ENSEMBL,raw_count[,2:ncol(raw_count)])
raw_count<-raw_count[,-1]

#对基因进行注释，其实这一步可以在这个地方省略，这一步在这儿做的目的是提前对比一下对照组和实验组的reads读数差别，提前进行一个验证，做差异分析之后还需要进行一次注释）
library('biomaRt')
#数据集的第一列的ID名为参照去搜索对应的基因名称
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl")) #注意使用人的基因库
my_ensembl_gene_id<-raw_count[,1]
mms_symbols<- getBM(attributes=c('ensembl_gene_id','external_gene_name'),filters = 'ensembl_gene_id', values = my_ensembl_gene_id, mart = mart)
#为了把基因的名称和数据集进行合并，要把两个数据集其中一列的名称改成一致的，比如都改成  ensembl_gene_id
colnames(raw_count)[1]<-("ensembl_gene_id")
readout=merge(mms_symbols,raw_count,by="ensembl_gene_id")

write.csv(raw_count,"6h_readout.csv",row.names=FALSE)

这样我们就获得了所有处理组的基因表达量文件。第一列是基因名字，第一行是处理组的名称。

2.开始分组对比

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("DESeq2")
library(tidyverse)
library(DESeq2)
setwd("D:/R/6h/tirmmed") 
getwd()   #查看文件路径
mycounts<-read.csv("6h_trimmed_readout.csv")
#把列里面的内容换成dataframe的列，然后删掉ensemble-id 和 gene-name两列内容
rownames(mycounts)<-mycounts[,1]
mycounts <- mycounts[, -1]
head(mycounts)


#condition这里是因子，不是样本名称；CC900011*2, control*3, GSKLSD1*3, LSD1IN20*2, OGD*3, SP2509*2, SP2577*2
condition <- factor(c(rep("OGD",3),rep("LSD1IN20",3),rep("GSKLSD1",3),rep("control",3),rep("CC900011",3),rep("SP2509",3),rep("SP2577",3)))
condition
colData <- data.frame(row.names=colnames(mycounts), condition)
colData

#注释：dds=DESeqDataSet Object
dds <- DESeqDataSetFromMatrix(mycounts, colData, design= ~ condition)
dds <- DESeq(dds)
dds
# 检查 colData 内容  
colData(dds)

#下面的根据不同组别修改，换两个组对比就修改这里！！！
res = results(dds, contrast=c("condition", "组1", "组2")) #最前面的condition不变，后面是想比较的水平，寻找两个组之间的差异基因
#根据p-value进行重新排序
res = res[order(res$pvalue),]
head(res)
##统计一下上调和下调的基因数
summary(res)
#观察矫正p值
table(res$padj<0.05)

#获取padj（p值经过多重校验校正后的值）小于0.05，表达倍数取以2为对数后大于1或者小于-1的差异表达基因
diff_gene_deseq2 <-subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
#看看筛出来了几个基因
dim(diff_gene_deseq2) 
#看看头几个结果
head(diff_gene_deseq2)

#用bioMart对差异表达基因进行注释
#两个文件没有共同的列名，所以要先给'diff_gene_deseq2'添加一个‘ensembl_gene_id’的列名
# 检查行名  
rownames(diff_gene_deseq2)  # 查看行名
# 假设第一列是基因 ID  
ensembl_gene_id<-rownames(diff_gene_deseq2)
diff_gene_deseq2<-cbind(ensembl_gene_id,diff_gene_deseq2)
colnames(diff_gene_deseq2)[1]<-c("ensembl_gene_id")
diff_name<-merge(diff_gene_deseq2,mms_symbols,by="ensembl_gene_id")
head(diff_name)

write.csv(diff_name,"6h_OGD_CC900011_readout.csv",row.names=FALSE)

这里通过res定义寻找哪两个不同处理组之间的差异基因，还有阈值的设定，如果前面的三件套没有把不同处理分开的话，这里的差异基因会非常难寻找的，数量很少。这时候输出的列表第一列就是差异基因的编号。

3.绘制热图与PCA

plotMA(res,ylim=c(-2,2))
topGene <- rownames(res)[which.min(res$padj)]
with(res[topGene, ], {
  points(baseMean, log2FoldChange, col="dodgerblue", cex=6, lwd=2)
  text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})

BiocManager::install("apeglm")  
library(apeglm)
res_order<-res[order(row.names(res)),]
res = res_order
res.shrink <- lfcShrink(dds,contrast = c("condition","OGD", "LSD1IN20"), res=res)
#但结果提示Error in lfcShrink(dds, contrast = c("condition", "control", "SP2577"),  :   type='apeglm' shrinkage only for use with 'coef'，目前不确定原因，故采取了下边两种方式
resLFC <- lfcShrink(dds, coef=2)
plotMA(resApe, ylim = c(-5,5))
topGene <- rownames(res)[which.min(res$padj)]
with(res[topGene, ], {
  points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
  text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})


#DESeq2提供了一个plotCounts()函数来查看某一个感兴趣的gene在组间的差别。counts会根据groups分组
d1=plotCounts(dds, gene="ENSMUSG00000024045", intgroup="condition", returnData=TRUE) 
ggplot(d1,aes(condition, count)) + geom_boxplot(aes(fill=condition)) + scale_y_log10() + ggtitle("Akap8")



#variance stabilizing transformation(VST)可消除方差对mean均值的依赖，尤其是低均值时的高log counts的变异。
vsdata <- vst(dds, blind=FALSE)
plotPCA(vsdata, intgroup="condition")



#绘制热图
library("pheatmap")
select<-order(rowMeans(counts(dds, normalized = TRUE)),
              decreasing = TRUE)[1:200]
df <- as.data.frame(colData(dds)[,c("condition","sizeFactor")])
#log2(n + 1)
ntd <- normTransform(dds)
pheatmap(assay(ntd)[select,], cluster_rows=FALSE, show_rownames=FALSE,
         cluster_cols=FALSE, annotation_col=df)

#加载R包，并看看都有啥可以选择的颜色
library("RColorBrewer")
display.brewer.all()
#转换数据还可以做出样本聚类热图。用dist函数来获得sample-to-sample距离。
sampleDists <- dist(t(assay(vsdata)))
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vsdata$condition, vsdata$type, sep="-")
colnames(sampleDistMatrix) <- NULL

colors <- colorRampPalette(brewer.pal(7,"Set3"))(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows=sampleDists,
         clustering_distance_cols=sampleDists,
         col=colors)

这里根据第二步定义的分析因子，直接开始所有处理组之间的对比，可以输出PCA，另外第一步合并的文件也可以在一些其他的网站上作图。如果二维PCA俩个轴的解释度＜50%，可以考虑使用三维PCA。

BioLadder-生物信息在线分析可视化云平台

https://zhuanlan.zhihu.com/p/313365092

4.基因富集分析

setwd("D:/R/6h/tirmmed") 
getwd()   #查看文件路径
# 重新打开 R 或 RStudio，并尝试更新包  
install.packages("cli")  
install.packages("colorspace")  
install.packages("glue")  
install.packages("rlang")


if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("org.Mm.eg.db")

BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")  

library(clusterProfiler)
library(DOSE)
library(org.Mm.eg.db)
library(org.Hs.eg.db) #这里需要导入和ENSEMBL编号对应的库
library(ggplot2)
library(stringr)

#加载数据，并挑选表达差异的基因ID，也可以用biomaRt包转换成gene name,用gene name 进行分析，这样就和参考链接里的语法一样了
sig.gene=read.csv("6h_OGD_SP2509_readout.csv")
deg <-subset(sig.gene,  padj<0.05 & abs(log2FoldChange)>=1)
id<-deg[,1]
head(id)
#这里需要选择和测序物种相同的库，比如人类和小鼠
id_human<- bitr(id , fromType = "ENSEMBL", toType = c("SYMBOL", "ENTREZID"), OrgDb = org.Hs.eg.db)
head(id_human)
#这一步后可以直接获得带有NCBI ID的文件，合并新信息到原始数据框##### 
deg_with_ids <- deg %>%  
  left_join(id_human, by = c("ensembl_gene_id" = "ENSEMBL"))  
head(deg_with_ids)  
write.csv(deg_with_ids, "6h_OGD_SP2509_readout_id.csv", row.names = FALSE)  

#提取非空非NA的Ensembl ID
valid_ids <- id_human$ENSEMBL[!is.na(id_human$ENSEMBL) & id_human$ENSEMBL != ""]  
head(valid_ids)
#如果换成SYMBOL的话
valid_ids <- id_human$SYMBOL[!is.na(id_human$SYMBOL) & id_human$SYMBOL != ""]  
head(valid_ids)
#接下来GO分析做图，GO可以在GO:BP(生物过程），GO:MF(分子功能），GO:CC(细胞组分）三个方面分别进行注释。这儿举一个MF的例子，换换语法就可以得到对应的BP & CC
ego_mf<-enrichGO(gene       = valid_ids,
                 OrgDb      = org.Hs.eg.db,  #这里选择匹配的库
                 keyType    = 'ENSEMBL',
                 ont        = "BP", #这里可以改成BP、MF、CC
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.05, #可以根据严格程度修改
                 qvalueCutoff = 0.05)
print(ego_mf)  
summary(ego_mf)
######运行后面的可以把每条GO和对应的基因合并成表格的一行####
ego_mf_readable <- setReadable(ego_mf, OrgDb = org.Hs.eg.db, keyType = "ENSEMBL")  
# 将可读对象转换为数据框  
ego_mf_df <- as.data.frame(ego_mf_readable)  
# 查看前几行结果  
head(ego_mf_df)  
# 保存结果到 CSV 文件  
write.csv(ego_mf_df, file = "6h_OGD_LSD1IN20_GO_BP_results_with_genes.csv", row.names = FALSE)

#绘制GO图###########
p <- ggplot(ego_mf, showCategory = 20) +  
  geom_bar(stat = "identity", aes(x = Description, y = Count, fill = p.adjust)) +  
  scale_fill_gradient(low = "red", high = "green", name = "p-value") +  
  scale_y_continuous(name = "Count") +  
  scale_x_discrete(name = "", labels = function(x) str_wrap(x, width = 50)) +  
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),  
        plot.title = element_text(hjust = 0.5),  
        panel.grid.major.y = element_blank(),  
        panel.grid.minor.y = element_blank()) +  
  coord_flip() +  
  ggtitle("The GO_BP enrichment analysis of all DEGs")  
print(p)  

#KEGG
kk<-enrichKEGG(gene = valid_ids,
               organism = 'hsa',#指定分析的生物体
               pvalueCutoff = 0.05)
kk[1:5]  #查看前5个KEGG结果

##### 将结果转换为数据框并保存为 CSV 文件######
kk_df <- as.data.frame(kk)  
write.csv(kk_df, file = "24h_control_SP2577_KEGG_analysis_results.csv", row.names = FALSE)  

# 绘制KEGG富集分析结果########  
p <- ggplot(kk, showCategory = 30) +  
  geom_bar(stat = "identity", aes(x = Description, y = Count, fill = p.adjust)) +  
  scale_fill_gradient(low = "red", high = "blue", name = "p-value") +  
  scale_y_continuous(name = "Count") +  
  scale_x_discrete(name = "", labels = function(x) str_wrap(x, width = 60)) +  
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),  
        plot.title = element_text(hjust = 0.5),  
        panel.grid.major.y = element_blank(),  
        panel.grid.minor.y = element_blank()) +  
  coord_flip() +  
  ggtitle("The KEGG enrichment analysis of all DEGs")  
print(p)  
#还可以绘制气泡图
p <- dotplot(kk, showCategory = 20) +  
  scale_size(range = c(2, 12)) +  
  scale_x_discrete(labels = function(x) str_trunc(x, width = 30)) +  
  theme(axis.text.y = element_text(angle = 45, hjust = 1, size = 8),  
        plot.title = element_text(hjust = 0.5),  
        panel.grid.major.x = element_blank(),  
        panel.grid.minor.x = element_blank())  
print(p)  

#Gene Set Enrichment Analysis（GSEA）
#需把gene name 转换好
genelist <- sig.gene$log2FoldChange
names(genelist) <- sig.gene[,8]  #我的数据集第8列是gene name
genelist <- sort(genelist, decreasing = TRUE)
head(genelist)
# GSEA分析，KeyType 也可以是 "ENSEMBL"，注意物种数据库的选择
gsemf <- gseGO(genelist, OrgDb = org.Hs.eg.db, keyType = "SYMBOL", ont="BP")
# 查看信息
head(gsemf)
# 画出GSEA图,从head(gsemf)里面选一条，注意GO:前后没有空格
gseaplot(gsemf, geneSetID="GO:0061635")

这下可以一并搞定GO和KEGG富集出图了，还可以做出特定GO条目的GSEA。