在GSE46224文件夹里,定量以后,获得(15+8)个fc.txt文件,这个是RNA-seq测序数据的表达矩阵,RNA-seq的count数据符合负二项分布,因此用SVA包时,要用ComBat-Seq函数(针对RNA-seq等离散型数据),不能用combat函数(针对芯片等连续型变量数据)。在RUVSeq包里,用RUVg函数。
分析流程如下:
1.R
读取二十三个样本、合并,获取未进行id转换的原始表达矩阵:

然后用biomaRt包整理基因注释信息,将geneid转换为gene名,获得genematrix表达矩阵

并储存。
save(genematrix,file = "genematrix.data.Rdata")
2.R
在genematrix里,行名为阿拉伯数字,第一列为gene名,要想进行后续分析,要把行名设置为基因名,但是获得的表达矩阵里,基因名有很多重复,如果设置为行名会报错,因此就要去重复。
去重复后获得genematrix2表达矩阵,行数从67103到38570,储存。

save(genematrix2,file = "genematrix2.data.Rdata")
3.R
这一步主要是基于RUVseq这个R包去批次以及差异分析
原始数据为genematrix2,这个矩阵是原始矩阵经过去重复及ID转换后的。
过滤数据:#去除表达量为0的,行数从38570到18561,得到genematrix3表达矩阵。
save(genematrix3,file = "genematrix3.data.Rdata")
标准化:使用上四分位数法,标准化前后分别画PCA图
用edgeR的方法找到DEG,挑选排名靠前的DEG,然后将排名靠后的DEG作为negative genes进行后续RUVg的normalization。
这时再画PCA图,HF和normal的样本已经被区分的非常开了。
最后再用edgeR进行GLM建模分析DEG,最后获得差异基因的矩阵DEG_result46224。

将行名变为第一列,最后变为下列形式:

4.R
这一步的主要目的:画出差异最大的几个基因;对任意基因画在两组里表达量的箱式图
所需输入文件:
1.分组文件:

2.差异基因

3.表达矩阵并进行归一化

利用这三者去画图

网友评论