DESeq2分析预处理

一、在R中去除低表达基因的一个常用方法是设定一个阈值,比如,一个基因在一定数量的样本中至少需要达到某个计数值。以下是一个示例R代码,展示如何去除在少于一半的样本中计数少于10的基因:

# Load necessary libraries
library(DESeq2)

# Read the count data from the CSV file without setting row names
counts_data <- read.csv("merged_count.csv", header = TRUE)

# Assume the first column contains the gene identifiers, so we skip it when applying the threshold
# Here we use a threshold of 10 counts, and the gene must be expressed at this level in at least half the samples
threshold <- 10
minimum_samples <- ceiling(ncol(counts_data) / 2)
filtered_counts <- counts_data[rowSums(counts_data[, -1] >= threshold) >= minimum_samples, ]

# Save the filtered data to a new CSV file
write.csv(filtered_counts, "filtered_merged_count.csv", row.names = FALSE)

在上面的代码中,rowSums(counts_data[, -1] >= threshold) 计算每一行(基因)有多少个样本满足计数值大于或等于阈值,然后检查这个数量是否至少为总样本数的一半。根据您的具体情况,您可能需要调整阈值和样本数的比例。

请注意,这个方法是基于一个简单的阈值,而实际中可能需要一个更复杂的方法来确定低表达基因,这可能涉及到基因的生物学知识和数据的具体情况。如果您的数据处理需要更专业的生物信息学方法,请咨询相关专家。

二、将csv文件中第一列的ENSG ID中的小数点及小数点后数字去掉

# 读取数据
counts_data <- read.csv("filtered_merged_count.csv")

# 使用sub函数去除第一列中的小数点和随后的数字
# 这里假设第一列的名称是"Geneid"
counts_data$Geneid <- sub("\\..*", "", counts_data$Geneid)

# 保存修改后的数据到新的CSV文件
write.csv(counts_data, "modified_filtered_merged_count.csv", row.names = FALSE)

三、将这个文件第一列的ENSG ID转化为GeneSymbol和GeneID,分别将每行转化后的结果放在本行的第二列和第三列

#######记得在packages中把biocManager和biomaRt加载(这个代码有时候不行,重启下Rstudio又可以了)

# 首先,安装并加载必要的包。如果您还没有安装biomaRt,请取消注释并运行下面的两行代码。
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("biomaRt")

# 加载biomaRt包
library(biomaRt)

# 设置ensembl mart
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")

# 读取CSV文件
counts_data <- read.csv("modified_filtered_merged_count.csv", stringsAsFactors = FALSE)

# 提取ENSG IDs
ensg_ids <- counts_data[, 1]

# 查询Ensembl数据库获取Gene Symbols 和 Entrez Gene IDs
annotations <- getBM(attributes = c('ensembl_gene_id', 'hgnc_symbol', 'entrezgene_id'),
                     filters = 'ensembl_gene_id', values = ensg_ids, mart = ensembl)

# 因为有时候annotations的顺序可能和原始的ensg_ids不一致,所以我们需要以ensg_ids为基础进行左连接
counts_data_annotated <- merge(counts_data, annotations, by.x = 'Geneid', by.y = 'ensembl_gene_id', all.x = TRUE, sort = FALSE)

# 将Gene Symbol和Gene ID插入到第一列后面
counts_data_annotated <- counts_data_annotated[, c(1, ncol(counts_data_annotated)-1, ncol(counts_data_annotated), 2:(ncol(counts_data_annotated)-2))]

# 为新插入的列重命名
names(counts_data_annotated)[2] <- "GeneSymbol"
names(counts_data_annotated)[3] <- "GeneID"

# 将新的数据框保存为CSV文件
write.csv(counts_data_annotated, "updated_modified_filtered_merged_count.csv", row.names = FALSE)

四、检测这个文件第一列各行有没有重复内容,并将重复行的内容只保留一行,把剩余重复行删除

# 读取CSV文件
data <- read.csv("/mnt/data/updated.csv", header = TRUE)

# 查找第一列中的重复项,并排除第一次出现的重复项
data_unique <- data[!duplicated(data[[1]]), ]

# 将处理后的数据写回到一个新的CSV文件
write.csv(data_unique, "/mnt/data/updated_unique.csv", row.names = FALSE)

把前四项代码合并优化了一下

counts_data <- read.csv("apln_count1.csv", header = TRUE)
# Assume the first column contains the gene identifiers, so we skip it when applying the threshold
# Here we use a threshold of 10 counts, and the gene must be expressed at this level in at least half the samples
minimum_samples <- ceiling(ncol(counts_data) / 2)
filtered_counts <- counts_data[rowSums(counts_data[, -1] >= 10) >= minimum_samples, ]

# 使用sub函数去除第一列中的小数点和随后的数字
# 这里假设第一列的名称是"Geneid"
filtered_counts$Geneid <- sub("\\..*", "", filtered_counts$Geneid)

# 设置ensembl mart
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")

ensg_ids <- filtered_counts[, 1]# 提取ENSG IDs

# 查询Ensembl数据库获取Gene Symbols 和 Entrez Gene IDs
annotations <- getBM(attributes = c('ensembl_gene_id', 'hgnc_symbol', 'entrezgene_id'),
                     filters = 'ensembl_gene_id', values = ensg_ids, mart = ensembl)
# 因为有时候annotations的顺序可能和原始的ensg_ids不一致,所以我们需要以ensg_ids为基础进行左连接
counts_data_annotated <- merge(filtered_counts, annotations, by.x = 'Geneid', by.y = 'ensembl_gene_id', all.x = TRUE, sort = FALSE)
# 将Gene Symbol和Gene ID插入到第一列后面
data_annotated <- counts_data_annotated[, c(1, ncol(counts_data_annotated)-1, ncol(counts_data_annotated), 2:(ncol(counts_data_annotated)-2))]
# 为新插入的列重命名
names(data_annotated)[2] <- "GeneSymbol"
names(data_annotated)[3] <- "GeneID"

checklwl <- data_annotated 
data_unique <- checklwl[!duplicated(checklwl[[1]]), ]
matrix <-data_unique[,c(2,9:23)]
sample_info <- data_unique[,1:3]
chr_info <- data_unique[,c(1,4:8)]

write.csv(matrix, "matrix1.csv", row.names = FALSE)
write.csv(sample_info, "sample_info1.csv", row.names = FALSE)
write.csv(chr_info, "chr_info1.csv", row.names = FALSE)

五、在R中使用 DESeq2 包进行差异表达分析通常涉及以下步骤:

加载必要的库(如 DESeq2 和其他数据处理工具)。
读取数据和样本信息。
创建 DESeq2 数据集对象。
运行差异表达分析。
提取结果并将其写入文件。
下面是一个基本的R脚本,包含执行上述步骤的代码:

注意第一列为基因ID,第二至四列为第一组三个重复,第五至七列为第二组三个重复。

代码中Control和Treatment分别代表第一组和第二组名称,第二至七列的第一行的建议改成Control1,Control2,Control3,Treatment1,Treatment2,Treatment3

建议将矩阵按照ENSG ID升序排列,删掉NA行

在运行此代码之前,请确保 updated_unique.csv 文件的格式符合要求。第一列应该是基因或特征的标识符,后面的列应该是计数数据。

# 加载所需的库
library(DESeq2)
library(tibble)

# 假设有两组样品,每组3个生物学重复
# 读取数据文件
data <- read.csv("matrix1.csv", header = TRUE)

# 保存基因ID列
Ensemble_id <- data[, 1]

# 样品信息,根据你的样品布局修改colData的内容
# 这里假设文件的列格式是:第一列为基因ID,后面六列分别为三个重复的两组样品
# 样本名称为 Sample1_rep1, Sample1_rep2, Sample1_rep3, Sample2_rep1, Sample2_rep2, Sample2_rep3
colData <- data.frame(
  sampleName = colnames(data)[-1],
  condition = rep(c("Control", "Treatment"), each = 3) 
)

# 为DESeq2创建矩阵,行为基因,列为样品
dds <- DESeqDataSetFromMatrix(
  countData = data[, -1],
  colData = colData,
  design = ~ condition
)

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

# 获取结果
res <- results(dds)

# 将基因ID与结果合并
output_data <- cbind(Ensemble_id, as.data.frame(res))

# 将合并后的数据保存为CSV文件
write.csv(output_data, file = "4v1.csv", row.names = FALSE)
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,332评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,508评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,812评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,607评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,728评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,919评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,071评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,802评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,256评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,576评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,712评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,389评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,032评论 3 316
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,798评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,026评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,473评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,606评论 2 350

推荐阅读更多精彩内容