1. 加载包
2. 导入数据
3. 注释基因
4. 去重
5. 过滤低表达基因
6. 无监督聚类MDS图
7. 保存数据

参考链接：
https://www.bioconductor.org/packages/devel/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow_CHN.html
https://cloud.tencent.com/developer/article/1422739

1. 加载包

if (!requireNamespace("BiocManager", quietly = TRUE))install.packages("BiocManager")
if (!require("airway", quietly = TRUE)) #使用数据集BiocManager::install("airway")
if (!require("org.Hs.eg.db", quietly = TRUE)) #注释基因BiocManager::install("org.Hs.eg.db")
if (!require("clusterProfiler", quietly = TRUE)) #ID转换BiocManager::install("clusterProfiler")
if (!require("edgeR", quietly = TRUE)) #过滤低表达基因BiocManager::install("edgeR")

2. 导入数据

options(stringsAsFactors = F)
library(airway)
data(airway)
# 获取基因counts矩阵
expr <- assay(airway)
expr[1:5,1:5]

#获取分组信息
group_list <- colData(airway)$dex
group_list

3. 注释基因

参考链接：http://www.360doc.com/content/19/0506/00/30846661_833639624.shtml

#查看支持的ID转换类型
keytypes(org.Hs.eg.db)

#转换ID
geneid<-bitr(rownames(expr), fromType='ENSEMBL', toType='SYMBOL', OrgDb='org.Hs.eg.db', drop = TRUE) #去除空值
geneid

#注释ID
expr = expr[ rownames(expr) %in% geneid[ , 1 ], ]
geneid = geneid[ match(rownames(expr), geneid[ , 1 ] ), ]
rownames(expr)<-geneid$SYMBOL
expr[1:5,1:5]

常用ID
1）Ensemble id： 欧洲生物信息数据库提供，一般以ENSG开头，后边跟11位数字。如TP53基因：ENSG00000141510
2）Entrez id： 美国NCBI提供，通常为纯数字。如TP53基因：7157
3）Symbol id： 文献中报道的基因名称。如TP53基因的symbol id为TP53
4）Refseq id： NCBI提供的参考序列数据库：可以是NG、NM、NP开头，分别代表基因，转录本和蛋白质。如TP53基因的某个转录本信息可为NM_000546
基因组注释包： http://bioconductor.org/packages/release/BiocViews.html

4. 去重

table(duplicated(rownames(expr))) #统计重复基因数目
#对重复基因名取平均表达量
if(sum(duplicated(rownames(expr)))>0) #判断是否有重复基因expr<-avereps(expr,ID=rownames(expr))

去重前：
去重后：

5. 过滤低表达基因

查看样本表达基因数目

apply(expr,2,fivenum)
apply(expr,2, function(x) sum(x>1))

从五个分位数统计量可以发现，每个样本含有部分表达值为0的基因，需要过滤

从上图可以发现平均每个样本被检测到的基因数量是15600左右

过滤方法1
过滤所选用的阈值具有主观性

#过滤在所有样本都是零表达的基因
expr=expr[apply(expr,1, function(x) sum(x>1) > 1),]
dim(expr)
apply(expr,2,fivenum)

过滤了31%的低表达基因

过滤方法2

x <-edgeR::DGEList(counts=expr)          #构建DEGList对象
keep.exprs <- edgeR::filterByExpr(x, group=group_list)#判断基因是否低表达
x <- x[keep.exprs,, keep.lib.sizes=FALSE]#过滤
head(keep.exprs)
dim(x)

ps：
1）默认选取最小的组内的样本数量为最小样本数，保留至少在这个数量的样本中有10个或更多序列片段计数的基因
2）对基因表达量进行过滤时使用CPM值而不是表达计数，以避免对总序列数大的样本的偏向性

可视化
使用的是过滤方法2得到的数据

#得到数据
library(RColorBrewer)
x <-edgeR::DGEList(counts=expr,group=factor(group_list))#构建DGE对象
L <- mean(x$samples$lib.size)/10^6
M <- median(x$samples$lib.size)/10^6
lcpm <- cpm(x, log=TRUE, prior.count=2)#log2(CPM + 2/L)
lcpm.cutoff <- log2(10/M + 2/L)        #过滤阈值
nsamples <- ncol(x)                 #提取样本数
col <- brewer.pal(nsamples, "Paired")  #颜色

lib.size为各个样本的总read counts

par(mfrow=c(1,2)) #一式两图
#过滤前图形
plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="")
abline(v=lcpm.cutoff, lty=3) #阈值线
for (i in 2:nsamples){
den <- density(lcpm[,i])
lines(den$x, den$y, col=col[i], lwd=2)
}
title(main="A. Raw data", xlab="Log-cpm") #标题和x轴名称
legend("topright",colnames(x), text.col=col, bty="n",cex=0.5) #图例
#过滤后图形
keep.exprs <- edgeR::filterByExpr(x, group=group_list)#判断基因是否低表达
x <- x[keep.exprs,, keep.lib.sizes=FALSE]#过滤
lcpm <- cpm(x, log=TRUE, prior.count=2)
plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="")
title(main="B. Filtered data", xlab="Log-cpm")
abline(v=lcpm.cutoff, lty=3)
for (i in 2:nsamples){
den <- density(lcpm[,i])
lines(den$x, den$y, col=col[i], lwd=2)
}
legend("topright",colnames(x), text.col=col, bty="n",cex=0.5)

注：卡阈值方法

6. 无监督聚类MDS图

参考链接：https://blog.csdn.net/yanta0/article/details/84798062

基本思想： 将高维坐标中的点投影到低维空间中，保持点彼此之间的相似性尽可能不变，而点与点之间的距离和对象之间的相似度高度相关
目的：
1）展示出样品间的相似性和不相似性
2）定性检测差异表达基因数目，越相似则差异基因越少
3）检查有无分组错误情况出现
理想情况： 样本会在不同的实验组内很好的聚类，且可以鉴别出远离所属组的样本

par(mfrow=c(1,1))
lcpm <- cpm(x, log=TRUE)
col.group <- group_list
levels(col.group) <- brewer.pal(nlevels(col.group), "Set1")
col.group<-droplevels(col.group) #去除不使用的水平因子
col.group <- as.character(col.group)
plotMDS(lcpm, labels=group, col=col.group)
title(main="Sample groups")

7. 保存数据

保留的是使用过滤方法1得到的数据

save(group_list,expr,file="expr.RData")

差异分析流程（一）数据预处理相关推荐

sklearn实战-----3.数据预处理和特征工程
1 概述 1.1 数据预处理与特征工程想象一下未来美好的一天,你学完了菜菜的课程,成为一个精通各种算法和调参调库的数据挖掘工程师了.某一天你从你的同事,一位药物研究人员那里,得到了一份病人临床表现 ...
python部分引入total值的问题_Python数据分析基础与过程综述，关键数据预处理异常点的发现与处理,python,及,流程,回顾,重点,之,值...
一. python数据分析基础库的导入基本是固定搭配 import numpy as np #科学计算基础库,多维数组对象ndarray import pandas as pd #数据处理库,Dat ...
珞珈1号-数据预处理流程
珞珈1号-数据预处理流程 1.重投影Albers 2.重采样 3.辐射校正–将INT32转化为浮点型真实数据 4.统一量纲(eg:和NPP同一量纲) 5.去噪 1.重投影Albers 参考这篇文章 2 ...
【mmdetection3d】——03自定义数据预处理流程
教程 3: 自定义数据预处理流程数据预处理流程的设计遵循一般惯例,我们使用 Dataset 和 DataLoader 来调用多个进程进行数据的加载.Dataset 将会返回与模型前向传播的参数所对 ...
基于FSL的DTI数据预处理流程
最近在学习处理DTI数据,总结了一份应用FSL做DTI数据预处理的流程与大家交流交流.如果有错误的地方欢迎大家指正! 我用的数据是Philips的数据,如果是GE或者西门子的数据可能会有所不同. 原始 ...
数据分析问题（异常值识别）中数据预处理部分流程（含2022年全国服务外包大赛实例）
博主个人理解的数据预处理主要包括个方面:读取文件 => 数据概览 => 缺失值填补 => 数据分布预览 => 衍生特征设计.这套流程在完成异常值识别时作为数据预处理时没有 ...
WGS完整流程介绍(原始数据质控、数据预处理、变异检测、数据注释)
一.原始数据质控 1.原始测序数据(也是reads) 从测序仪中直接取下来的数据,它包括了所有的碱基,无论是测序质量低的,还有可能包含测错的,可能还会包含实验误差. 2.数据质控 ...
【理论】数据预处理流程
文章目录 1.找数据集 2.理解数据 3.数据处理 1.找数据集已经有数据集的跳过这一步. 找到合适的数据集.如何找数据集请查看一些其他教程. 2.理解数据这一步主要是对自己找到的数据集要有一个总 ...
matlab预处理光谱数据,一种近红外光谱数据预处理方法与流程
本发明公开了属于近红外光谱分析技术领域,尤其涉及一种近红外光谱数据预处理的新方法,主要用于建立近红外定量和定性模型时对近红外光谱数据的预处理. 背景技术: 近红外光谱技术具有分析速度快.样本制作简单的 ...
差异分析流程（二）limma
1. 制作三个矩阵 2. 建模前归一化 3. 建模及结果 4.检查参考链接: https://www.bioconductor.org/packages/devel/workflows/vignet ...

差异分析流程（一）数据预处理