转录组分析——差异表达分析

本文介绍转录组差异表达分析,数据名为“GSE124133_gene_count_matrix.txt”,整理分类数据集condition type.txt,利用R语言中DESeq包进行分析,结果文件可根据实验情况筛选处理,数据来源于相关文献。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

1.数据的类型,数据名称为“GSE124133_gene_count_matrix.txt”
在这里插入图片描述
2.整理分类数据集condition type.txt,数据如下:
在这里插入图片描述
3.利用R语言中DESeq包进行差异表达分析

#查看R语言当前所在的文件夹
getwd()
#进入到数据所在的文件夹所在的路径(这一步很重要,如果不正确,可能导致数据文件无法正常的导入)
setwd("E:/shengwuxinxi/基因组重测序数据操作步骤/作业")
#如果以前没有使用过DESeq,则需要到Bioconductor上查找你的R版本所对应的安装包方式[Bioconductor](http://www.bioconductor.org/)
#加载安装包
library(DESeq)
library(cluster)
library(DESeq2)
#载入"GSE124133_transcript_count_matrix.txt“数据
database_all <- read.table(file = "GSE124133_transcript_count_matrix.txt", sep = "\t", header = T, row.names = 1)
head(database_all)
#载入分组情况文件condition type.txt
colDataMOTSHFD<-read.table("condition type.txt",header=T,sep="\t")
head(colDataMOTSHFD)
##DESeq利用进行差异表达分析分析
ddsMOTSHFD  <-  DESeqDataSetFromMatrix( countData = database_all, colData = colDataMOTSHFD, design =~sex+pdj+status)
dds <- DESeq(ddsMOTSHFD)
res <- results(dds)
resdata<merge(as.data.frame(res),as.data.frame(counts(dds,normalized=TRUE)),by="row.names",sort=FALSE)
head(resdata)
#将计算结果文件,导出表格形式
write.csv(resdata,file = "3-tumorous_v_survive.csv")

4.结果文件为,可以根据实验的具体情况对数据进行筛选,和处理
在这里插入图片描述
说明:本文数据来源于文献《Integrated analysis of lncRNA and mRNA repertoires in Marek’s disease infected spleens identifies genes relevant to resistance》

### 有参转录组分析流程中的表达矩阵与TPM矩阵 #### 表达矩阵的构建 在有参转录组分析中,表达矩阵是通过将测序读段映射到已知参考基因组来生成的。每行代表一个基因,而每列则对应不同的样本。此矩阵记录了各基因在不同样本中的原始计数(counts),即每个基因被检测到的读段数量。 对于这一过程,通常会采用如下方法: - **比对工具的选择**:可以利用`Hisat2`或`STAR`这样的短读段比对器将RNA-seq数据映射至参考基因组上[^4]。 - **量化步骤**:一旦完成了读段的定位工作之后,则可以通过诸如`featureCounts`之类的程序统计落在各个特征区域内的唯一匹配上的读取数目,从而形成初步的count matrix——也就是所谓的“counts矩阵”。 ```bash # 使用 featureCounts 进行定量计算 featureCounts -a annotation.gtf \ -o output_counts.txt \ aligned_reads.bam ``` #### TPM矩阵的转换 为了更好地比较跨样品间的基因表达水平差异,在得到初始的counts矩阵基础上还需要进一步标准化处理,其中一种常见的方式就是转化为Transcripts Per Million (TPM) 值的形式。这样做不仅考虑到了总测序量的影响因素,同时也校正了因目标片段长度不一致所带来的偏差问题。 具体来说,TPM 的计算公式为: \[ \text{TPM}_i = (\frac{\sum_j C_{ij}}{L_i})\times10^6 \] 这里 \(C\) 是指某个特定基因\(j\)下的read counts;\(L\)表示相应转录本的有效长度(通常是减去接头序列后的实际碱基数)。因此,当我们将所有的reads分配给它们所属的transcript后就可以按照上述方程求得最终的结果并构成新的TPM矩阵。 ```r library(edgeR) # 加载 count data 文件 d <- read.table("output_counts.txt", header=TRUE, row.names="Geneid") # 计算有效长度向量 L (假设我们有一个名为 'lengths' 的文件) gene_lengths <- scan('lengths', what='numeric') # 创建 DGEList 对象用于后续操作 y <- DGEList(counts=d, genes=gene_lengths) # 应用 TMM 归一化因子调整 y_tmm <- calcNormFactors(y, method="TMM") # 输出 TPM 矩阵 tpm_matrix <- cpm(y_tmm, log=FALSE)*1e6 / colSums(cpm(y_tmm)) write.csv(tpm_matrix,"tpm_output.csv") ```
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值