生信上游分析案例 Ｒ代码

由于生信分析案例种类繁多，无法一一列举所有的R代码。以下是一个简单的RNA-seq数据分析案例，供参考：

1. 导入数据并进行质控

library(DESeq2)

# 导入表达矩阵
counts <- read.table("counts.txt", header=TRUE, row.names=1)

# 导入样本信息
colData <- read.table("sample_info.txt", header=TRUE, row.names=1)

# 创建DESeq2对象
dds <- DESeqDataSetFromMatrix(countData=counts, colData=colData, design=~condition)

# 进行质控
plotPCA(rlog(dds), intgroup="condition")
plotHeatmap(counts(dds), cluster_rows=TRUE)

2. 进行差异表达分析

# 进行差异表达分析
dds <- DESeq(dds)

# 获取差异表达基因
res <- results(dds)

# 过滤p值和fold change
res_sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)

# 导出差异表达基因列表
write.table(res_sig, "DE_genes.txt", sep="\t", quote=FALSE)

3. 进行富集分析

# 导入基因列表和物种信息
gene_list <- read.table("DE_genes.txt", header=TRUE, row.names=1)
species <- "Homo sapiens"

# 进行富集分析
library(clusterProfiler)
enrich_res <- enrichGO(gene_list=row.names(gene_list), 
                       OrgDb=org.Hs.eg.db, 
                       keyType="SYMBOL", 
                       universe=row.names(counts), 
                       ont="BP",
                       pvalueCutoff=0.05,
                       qvalueCutoff=0.2,
                       readable=TRUE)

# 导出富集分析结果
write.table(as.data.frame(enrich_res), "enrichment_results.txt", sep="\t", quote=FALSE, row.names=FALSE)

4. 进行基因表达量可视化

# 绘制基因表达量热图
library(pheatmap)
pheatmap(counts, cluster_rows=TRUE, show_rownames=FALSE)

# 绘制差异表达基因的表达量热图
library(RColorBrewer)
cols <- colorRampPalette(brewer.pal(9, "Set1"))(length(levels(factor(colData$condition))))
pheatmap(counts[row.names(res_sig),], cluster_rows=TRUE, show_rownames=FALSE, col=cols)