Loading [MathJax]/jax/output/CommonHTML/config.js
首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >专栏 >如何去掉数据中的离群样本?

如何去掉数据中的离群样本?

作者头像
生信菜鸟团
发布于 2024-05-11 10:00:59
发布于 2024-05-11 10:00:59
83400
代码可运行
举报
文章被收录于专栏:生信菜鸟团生信菜鸟团
运行总次数:0
代码可运行

引言

当我们拿到一组数据想要开始分析时,做的第一件事情就是质控,看一下数据怎么样,是否适用于我们的分析流程,以及某些低表达或极端表达的基因和样本是否应该删除更利于分析结果。今天分享一下如何删除离群样本,并探索一下是否有生物学意义。

自己的表达量矩阵数据绘制主成分分析图

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
#加载Rlibrary("FactoMineR")
library("factoextra")  
#载入数据
load(file = 'symbol_matrix.Rdata') 
symbol_matrix[1:4,1:4] ## 基因名字的样品,癌症转录组count矩阵     
## TCGA-97-7938-01A TCGA-55-7574-01A TCGA-05-4250-01A
## A1BG                   15               17               18
## A1BG-AS1              191               90              104
## A1CF                    4                5                0
## A2M                 87438            61893            55347
## TCGA-95-A4VK-01A
## A1BG                   17
## A1BG-AS1               53
## A1CF                    2
## A2M                 42664
#转换成cpm数据
dat = log2(edgeR::cpm(symbol_matrix)+1)
dat[1:4,1:4]
## TCGA-97-7938-01A TCGA-55-7574-01A
## A1BG            0.4834302        0.6586411
## A1BG-AS1        2.6013824        2.0225987
## A1CF            0.1455475        0.2267243
## A2M            11.1807753       11.0413363
## TCGA-05-4250-01A TCGA-95-A4VK-01A
## A1BG            0.3420593       0.45300345
## A1BG-AS1        1.3481920       1.10437674
## A1CF            0.0000000       0.06129024
## A2M             9.6860041       9.85607747
table(group_list)
## group_list
## Squamous_cell_carcinoma          Adenocarcinoma 
## 501                     526 
pro = 'test'
exp=t(dat)#画PCA图时要求是行名时样本名,列名时探针名,因此此时需要转换
exp=as.data.frame(exp)#将matrix转换为data.frame 
#主成分分析
dat.pca <- PCA(exp , graph = FALSE)
this_title <- paste0(pro,'_PCA')
p<-fviz_pca_ind(dat.pca,
                   geom.ind = "point",  #c( "point", "text" ), # show points only (nbut not "text")
                   col.ind = group_list, # color by groups
                   palette = "Dark2",
                   addEllipses = TRUE, # Concentration ellipses
                   legend.title = "Groups")+
  ggtitle(this_title)+ 
  theme(plot.title = element_text(size=12,hjust = 0.5))
p

我们可以看到有几个样本很明显散在椭圆之外,我们现在通过第一次pca分析的结果将其删除,看是否会对后续的分析有影响。

02

PCA删除离群样本

删除距离太远的样本,上面的pca绘图的时候其实也返回来了横纵坐标信息:

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
#筛选离群样本名称
name<-as.character(p2$data$name[p$data$x>600|p$data$x>600])#PCA图中x或y轴大于600的视为离群样本
name
## [1] "TCGA-44-5645-01B" "TCGA-44-3918-01B" "TCGA-44-2656-01B"
## [4] "TCGA-44-6775-01C" "TCGA-44-6146-01B" "TCGA-44-2662-01B"
## [7] "TCGA-44-3917-01B" "TCGA-44-2668-01B" "TCGA-44-4112-01B"
## [10] "TCGA-44-2666-01B" "TCGA-44-6147-01B" "TCGA-21-5782-01A"
name_index <- which(rownames(exp) %in% name)
#在基因矩阵及分组中删除离群样本
exp1 <- exp[-name_index, ]
group_list1<-group_list[-name_index]
#重新绘制PCA图
dat.pca1 <- PCA(exp1 , graph = FALSE)
this_title1 <- paste0(pro,'1_PCA')
p1<- fviz_pca_ind(dat.pca1,
                   geom.ind = "point",  #c( "point", "text" ), # show points only (nbut not "text")
                   col.ind = group_list1, # color by groups
                   palette = "Set1",
                   addEllipses = TRUE, # Concentration ellipses
                   legend.title = "Groups")+
  ggtitle(this_title)+ 
  theme(plot.title = element_text(size=12,hjust = 0.5))
p1

删除离群样本重新绘图之后已经没有距离很远的样本了,也更好看一些。

其实除了PCA图,还有WGCNA的层次聚类也可以实现这一过程。

03

层次聚类可视化

绘制层次聚类图

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
#hclust
library(WGCNA)
EXP<-exp
rownames(EXP)<-str_sub(rownames(EXP), 6, 16)#简化一下样本名称
sampleTree = hclust(dist(EXP), method = "average")
#绘图
png(file = "sampleClustering.png", width = 26000, height =1500,res = 300)
par(cex = 0.6)
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)
dev.off()

数据样本量较大,所以截取一部分,只有这几个样本是单独一个分支,我们可以把这些异常样本的分支切除。

增加切除线

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
png(file = "sampleClustering1.png", width = 26000, height =1500,res = 300)
par(cex = 0.6)
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)
abline(h = 270, col = "red")
dev.off()

切除分支

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
clust = cutreeStatic(sampleTree, cutHeight = 270,minSize = 10)#切除高度设置为270
  table(clust) # 2代表切除的,1代表保留的
## clust
## 1    2 
## 1016   11 
  keepSamples = (clust!=2)
  EXP1 = EXP[keepSamples, ]
#看一下被切除的样本分支
rownames(EXP)[!keepSamples]
## [1] "44-5645-01B" "44-3918-01B" "44-2656-01B" "44-6775-01C"
## [5] "44-6146-01B" "44-2662-01B" "44-3917-01B" "44-2668-01B"
## [9] "44-4112-01B" "44-2666-01B" "44-6147-01B" 
#重新绘制层次聚类图
sampleTree1 = hclust(dist(EXP1), method = "average")
png(file = "sampleClustering2.png", width = 26000, height =1500,res = 300)
par(cex = 0.6)
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)
dev.off()

现在就没有单独分支的样本啦~而且和PCA图删除的样本几乎是一样的。

那么这个步骤到底有没有生物学意义呢?我们接下来继续探究。

04

差异分析结果比较

两组数据分别用的DESeq2包进行差异分析(这个代码省略,因为太简单了),有了差异结果矩阵,就可以比较一下删除离群样本之后是否会对差异分析的结果产生影响。

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
#导入差异分析结果
load(file = 'DEG_deseq2.Rdata')#原始数据
summary(DEG_deseq2)
deg_DESeq2 = na.omit(as.data.frame(DEG_deseq2[order(DEG_deseq2$padj),]))  
load(file = 'DEG1_deseq2.Rdata')#去除样本之后
summary(DEG1_deseq2)
deg1_DESeq2 = na.omit(as.data.frame(DEG1_deseq2[order(DEG1_deseq2$padj),]))  
#绘制散点图比较一下两者的log2FoldChange
colnames(deg1_DESeq2)
ids=intersect(rownames(deg_DESeq2),
              rownames(deg1_DESeq2))
df= data.frame(
  deg_DESeq2 = deg_DESeq2[ids,'log2FoldChange'],
  deg1_DESeq2 = deg1_DESeq2[ids,'log2FoldChange']
)
library(ggpubr)
ggscatter(df, x = "deg_DESeq2", y = "deg1_DESeq2",
          color = "black", shape = 21, size = 3, # Points color, shape and size
          add = "reg.line",  # Add regressin line
          add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
          conf.int = TRUE, # Add confidence interval
          cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
          cor.coeff.args = list(method = "pearson",  label.sep = "\n"))

使用的数据有1027个样本,只是删除了PCA中的12个样本,所以看起来影响不大,那么我们再考虑他的统计学意义,结合P值看一下对差异基因是否有影响。

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
#设置阈值
pvalue_t=0.05
logFC_t = 1
deg_DESeq2$g=ifelse(deg_DESeq2$padj > pvalue_t,'stable', #if 判断:如果这一基因的P.Value>0.01,则为stable基因
                    ifelse( deg_DESeq2$log2FoldChange > logFC_t,'up', #接上句else 否则:接下来开始判断那些P.Value<0.01的基因,再if 判断:如果logFC >1.5,则为up(上调)基因
                            ifelse( deg_DESeq2$log2FoldChange < -logFC_t,'down','stable') )#接上句else 否则:接下来开始判断那些logFC <1 的基因,再if 判断:如果logFC <1,则为down(下调)基因,否则为stable基因
) 
deg1_DESeq2$g=ifelse(deg1_DESeq2$padj > pvalue_t,'stable', #if 判断:如果这一基因的P.Value>0.01,则为stable基因
                    ifelse( deg1_DESeq2$log2FoldChange > logFC_t,'up', #接上句else 否则:接下来开始判断那些P.Value<0.01的基因,再if 判断:如果logFC >1.5,则为up(上调)基因
                            ifelse( deg1_DESeq2$log2FoldChange < -logFC_t,'down','stable') )#接上句else 否则:接下来开始判断那些logFC <1 的基因,再if 判断:如果logFC <1,则为down(下调)基因,否则为stable基因
) 
pdf('figures/balloonplot.pdf',width = 6,h=4)
gplots::balloonplot(
  table(  deg_DESeq2[ids,'g'],
          deg1_DESeq2[ids,'g'])
)
dev.off()

从比较的表格中可以看出删除样本之后上调的差异基因减少了将近一半,下调的差异基因相差不大,那么删除的样本影响了什么导致的这个结果呢?

以我们最常研究的编码mRNA为起点,看一下是否也是有同样的结果。

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
library(AnnoProbe) 
ids=annoGene( ids,'SYMBOL','human')
head(ids) 
tail(sort(table(ids$biotypes)))
#筛选编码基因
pbGenes = unique(ids[ids$biotypes=='protein_coding',1])
deg_DESeq2_pb<-deg_DESeq2[pbGenes,]
deg1_DESeq2_pb<-deg1_DESeq2[pbGenes,]
pdf('figures/balloonplot3.pdf',width = 6,h=4)
gplots::balloonplot(
  table(  deg_DESeq2_pb[,'g'],
          deg1_DESeq2_pb[,'g'])
   )
dev.off()

上下调基因列表重合度很好,可见,异常样本基本上只是影响了非编码mRNA,对编码mRNA并没有太大的影响。

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2024-05-09,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信菜鸟团 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
暂无评论
推荐阅读
编辑精选文章
换一批
WGCNA仅仅是划分基因模块,其它都是附加分析
曾老师给我分享了一篇数据挖掘的文章,里面的WGCNA非常奇怪,我之前没见过这样的模块与表型的相关性热图
生信技能树
2023/09/04
1.5K0
WGCNA仅仅是划分基因模块,其它都是附加分析
三种转录组差异分析方法及区别你会了吗?
在做项目时,曾有小伙伴对我用edgeR进行差异分析筛选出的具体显著差异基因表示质疑,因为发表的文章清楚的说明某个基因是差异基因,但是我edgeR的分析结果并没有表明。在小伙伴的质疑下,我认真看了下文章,发现文章用的是DEseq2进行差异分析。值得注意的是该小伙伴关注的差异基因是一个离散比较大的基因,此处的离散较大可以理解为假定对照组为5,6,7;实验组则为14,13,3的情况。那为什么这个基因在edgeR分析下不是显著差异基因,然而在DEseq2的分析下是差异基因呢?这应该很大程度源于算法判定显著差异基因的区别。接着,我看了关于DEseq2与edgeR区别的描述,发现「edgeR与Deseq2都是基于负二项分布模型做的,两者处理同一组数据时,相同阈值处理大部分基因是一样的,但是也会有一部分基因会因为离散度不同导致差异不同」,如刚刚示例的基因离散度被DEseq2识别为差异,但是不被edgeR识别,所以两种算法获取的差异基因与数目是存在细微区别的。
生信菜鸟团
2022/10/31
6.1K1
三种转录组差异分析方法及区别你会了吗?
GEO数据挖掘之转录组测序数据流程-以GSE150392为例
这个包里可以画pca, 热图,火山图,韦恩图,具体每个图的算法,可以看生信技能树GEO芯片分析
生信技能树
2022/06/08
2.6K1
GEO数据挖掘之转录组测序数据流程-以GSE150392为例
比较不同流程(limma/voom,edgeR,DESeq2 )差异分析的区别
距离第一次听说生信已经十几年了,现在是邋遢大叔重新开始学代码,精力确实已不像从前,各位入坑还是要乘早。后来约莫在5年前,课题组当时有个RNA-Seq数据,lab meeting时听瑞典小哥在汇报DEGs筛选,当时感觉好是神奇。其实陆陆续续也有过学习的念头,但在对自己的各种纵容下,想法又逐渐隐没。直到2月前,机缘巧合参加了生信技能树培训,才进一步强化了自己学习生信技术的信念。
生信技能树
2021/01/05
5.3K0
初探mRNA、lncRNA联合分析之下游
虽然这个项目是在转录本水平上开展的研究,但既然我们拿到了基因表达矩阵,也干脆看一看一些基本情况,这个部分代码此处省略,基本上和后面的转录本水平对应代码,包括使用的封装函数,是一致的
生信菜鸟团
2023/08/23
8150
初探mRNA、lncRNA联合分析之下游
转录组差异分析—基本流程
读取RawCounts.csv文件,其文件形式如下图行名为ensembleid,列名为样本名称。
sheldor没耳朵
2024/07/29
3020
转录组差异分析—基本流程
GEO数据库中芯片数据分析思路
AnnoProbe是曾建明老师2020年开发的一款用于下载GEO数据集并注释的R包,收录在tinyarray里。 idmap##根据所给的GPL号,返回探针的注释 geoChina##根据所给的GSE号,下载对应的表达矩阵 annoGene##根据gencode中的GTF文件注释基因ID
小张小张
2023/05/25
2K0
GEO多数据集联合分析-文献复现
文献题目:基于生物信息学的新型铁死亡基因生物标志物和免疫浸润谱在糖尿病肾病中的应用Huang, Y., & Yuan, X. (2024). Novel ferroptosis gene biomarkers and immune infiltration profiles in diabetic kidney disease via bioinformatics. FASEB journal : official publication of the Federation of American Societies for Experimental Biology, 38(2), e23421. https://doi.org/10.1096/fj.202301357RR. IF: 4.8 Q1
用户11008504
2024/05/11
4452
TCGA数据整理-3
生信菜鸟团
2024/07/11
1470
TCGA数据整理-3
一起画个圈圈看差异基因
最近朋友看论文,看到了个展示差异基因的好看图,说想给自己的差异基因也来画一个,我研究了下,实现挺简单,现成的R包circlize 就可以做,那我们就一起来画一个圈圈吧!
生信菜鸟团
2022/04/08
1.1K0
一起画个圈圈看差异基因
并不需要得到去批次后的表达量矩阵
但是很多人理解的"去批次效应"(batch effect removal)这个操作应该是会输入一个表达量矩阵,然后输出一个表达量矩阵。其实在单细胞转录组数据分析里面并不是这样的,比如我们常见的harmony操作,它针对的就并不是原始的单细胞转录组表达量矩阵(几万个基因几万个细胞),而是pca分析结果(还是几万个细胞但是只有少量的pc)。这样的话,harmony操作后并没有修改我们的原始的单细胞转录组表达量矩阵,这一点可能会确实是让大家困惑。
生信技能树
2024/05/06
3090
并不需要得到去批次后的表达量矩阵
两分组差异分析的上下调基因跟PCA分析的主成分基因的交集如何
学员提出来了一个很有意思的问题, 就是表达量矩阵质量控制环节里面的PCA分析,可以看出来不同样品在二维图里面的距离,其中PC1是可以区分两个分组,如下所示:
生信技能树
2023/09/04
5650
两分组差异分析的上下调基因跟PCA分析的主成分基因的交集如何
生信技能树GEO数据挖掘直播配套笔记
二代测序(RNA_seq):如果是counts 可选择limma的voom算法或者edgeR或者DESeq2。 如果是FPKM或TPM可选择limma,注意:edgeR和DESeq2只能处理count注意:count做差异分析计算上下调,FPKM或TPM进行下游可视化
生信技能树
2022/06/08
2.1K0
生信技能树GEO数据挖掘直播配套笔记
生信技能树——GEO转录组RNA_seq_GSE162550
和生信技能树GEO转录组“GSE150392“分析类似,唯一区别就是在数据处理和ID转换这一环节略微有区别
生信技能树
2022/06/08
1.8K0
生信技能树——GEO转录组RNA_seq_GSE162550
GEO数据库挖掘
输入数据是数值型矩阵/数据框,颜色的变化表示数值的大小。有相关性热图和差异基因热图。
叮当猫DDM
2023/07/16
8491
GEO表达芯片数据分析
---title: "GEO表达芯片数据分析"output: html_documentdate: "2023-03-20"---关于该流程代码的说明:(1)本流程仅适用于GEO芯片表达数据,以"GSE56649"为例(2)先在GEO数据库中确定是否为"Expression profiling by array",不是的话不能使用本流程!(3)注意需要自行修改或判断的代码一般放在了两个空行之间(4)代码的注释有一丢丢多,目的是为了更好地帮助大家理解1.下载数据,提取表达矩阵、临床信息和GPL编号rm(lis
小叮当aka
2023/03/23
3.3K1
TCGA分析-数据整理2
素素
2023/10/31
3780
指定通路绘制gsea图热图和火山图
但是我们直接是对gsea分析结果的最终es值在可视化,所以是行是通路,列是癌症的,数值是gsea的es打分的矩阵。对初学者来说, 跳过了大量细节,所以跟这个教程会比较吃力,有粉丝就提问了希望可以对这些通路在在具体的癌症里面细化展示,比如绘制gsea图,热图和火山图。
生信技能树
2022/07/26
2.5K0
指定通路绘制gsea图热图和火山图
基于count数据的基因差异表达分析万能代码
关于差异分析的文章中【一文就会TCGA数据库基因表达差异分析】其实有推送过,这篇文章目前为止,有近千人付费学习。
DoubleHelix
2022/06/13
3.9K0
基于count数据的基因差异表达分析万能代码
转录组测序结果分析
#没有任何提示就是成功了,如果有warningxx包不存在,用library检查一下。
叮当猫DDM
2024/03/19
3031
相关推荐
WGCNA仅仅是划分基因模块,其它都是附加分析
更多 >
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档