一、在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)