依旧是曾老师的学徒作业,分享作业过程中发现的现象。
作业链接如下:
《circRNA芯片也是同样的差异分析》circRNA芯片也是同样的差异分析 (qq.com)
总结:做一个差异分析和原文比较一下
1- 文章机构,是兰州第一医院
2- 数据集 GSE97332
3- 看一下差异分析图
ok,开始处理数据吧。
看文献知道是用Agilent的设备扫描的,不存在双面的情况。
#先简单加载一堆包
rm(list=ls())
options(stringsAsFactors=F)
#获取路径,方便复制粘贴
require("limma")
getwd()
raw.files <- dir(path="./GSE97332_RAW", pattern=".txt.gz$")
raw.files
# [1] "GSM2561829_L72763.txt.gz" "GSM2561830_N69749.txt.gz" "GSM2561831_N73936.txt.gz" "GSM2561832_N73970.txt.gz"
# [5] "GSM2561833_N74000.txt.gz" "GSM2561834_N74009.txt.gz" "GSM2561835_N74013.txt.gz" "GSM2561836_T72763.txt.gz"
# [9] "GSM2561837_T69749.txt.gz" "GSM2561838_T73936.txt.gz" "GSM2561839_T73970.txt.gz" "GSM2561840_T74000.txt.gz"
# [13] "GSM2561841_T74009.txt.gz" "GSM2561842_T74013.txt.gz"
Eraw <- read.maimages(files=raw.files, source="agilent", path="GSE97332_RAW", green.only=T)#单色,所以T
# Read GSE97332_RAW/GSM2561829_L72763.txt.gz
# Read GSE97332_RAW/GSM2561830_N69749.txt.gz
# Read GSE97332_RAW/GSM2561831_N73936.txt.gz
# Read GSE97332_RAW/GSM2561832_N73970.txt.gz
# Read GSE97332_RAW/GSM2561833_N74000.txt.gz
# Read GSE97332_RAW/GSM2561834_N74009.txt.gz
# Read GSE97332_RAW/GSM2561835_N74013.txt.gz
# Read GSE97332_RAW/GSM2561836_T72763.txt.gz
# Read GSE97332_RAW/GSM2561837_T69749.txt.gz
# Read GSE97332_RAW/GSM2561838_T73936.txt.gz
# Read GSE97332_RAW/GSM2561839_T73970.txt.gz
# Read GSE97332_RAW/GSM2561840_T74000.txt.gz
# Read GSE97332_RAW/GSM2561841_T74009.txt.gz
# Read GSE97332_RAW/GSM2561842_T74013.txt.gz
E.bg <- backgroundCorrect(Eraw, method="normexp",offset=50) #背景清除
# Array 1 corrected
# Array 2 corrected
# Array 3 corrected
# Array 4 corrected
# Array 5 corrected
# Array 6 corrected
# Array 7 corrected
# Array 8 corrected
# Array 9 corrected
# Array 10 corrected
# Array 11 corrected
# Array 12 corrected
# Array 13 corrected
# Array 14 corrected
E.norm <- normalizeBetweenArrays(log2(E.bg$E), method="quantile") #表达矩阵取log
range(E.norm)
# [1] 5.655481 16.505287
E.norm[1:4,1:4]
# GSM2561829_L72763.txt GSM2561830_N69749.txt GSM2561831_N73936.txt GSM2561832_N73970.txt
# [1,] 16.442645 16.442645 16.468603 16.356482
# [2,] 5.955997 6.072796 5.920535 5.978874
# [3,] 5.907458 6.093137 5.929130 6.071116
# [4,] 7.309498 7.614550 7.352912 7.113876
library("stringr")
sample_name <- str_extract(raw.files,"GSM\\d*")
sample_name
# [1] "GSM2561829" "GSM2561830" "GSM2561831" "GSM2561832" "GSM2561833" "GSM2561834" "GSM2561835" "GSM2561836" "GSM2561837"
# [10] "GSM2561838" "GSM2561839" "GSM2561840" "GSM2561841" "GSM2561842"
colnames(E.norm)=sample_name
dim(E.norm)
# [1] 15743 14
E.norm[1:4,1:4]
# GSM2561829 GSM2561830 GSM2561831 GSM2561832
# [1,] 16.442645 16.442645 16.468603 16.356482
# [2,] 5.955997 6.072796 5.920535 5.978874
# [3,] 5.907458 6.093137 5.929130 6.071116
# [4,] 7.309498 7.614550 7.352912 7.113876
rownames(E.norm) <- E.bg[["genes"]][["SystematicName"]]## 改成探针名
dim(E.norm)
# [1] 15743 14
E.norm[1:4,1:4]
# GSM2561829 GSM2561830 GSM2561831 GSM2561832
# GE_BrightCorner 16.442645 16.442645 16.468603 16.356482
# DarkCorner 5.955997 6.072796 5.920535 5.978874
# DarkCorner 5.907458 6.093137 5.929130 6.071116
# ASCRP000245 7.309498 7.614550 7.352912 7.113876
E.nc <- E.norm[Eraw$genes$ControlType==0,] #选择ControlType为0的行
dim(E.nc)
# [1] 15207 14
E.nc[1:4,1:4]
# GSM2561829 GSM2561830 GSM2561831 GSM2561832
# ASCRP000245 7.309498 7.614550 7.352912 7.113876
# ASCRP003158 6.774156 6.840387 6.688771 6.368704
# ASCRP002115 5.869470 5.995311 5.869927 5.923229
# ASCRP001695 11.890401 12.229166 11.843689 12.888761
exp=E.nc
dim(exp);range(exp);exp[1:4,1:4]
# [1] 15207 14
# [1] 5.655481 16.505287
# GSM2561829 GSM2561830 GSM2561831 GSM2561832
# ASCRP000245 7.309498 7.614550 7.352912 7.113876
# ASCRP003158 6.774156 6.840387 6.688771 6.368704
# ASCRP002115 5.869470 5.995311 5.869927 5.923229
# ASCRP001695 11.890401 12.229166 11.843689 12.888761
###### 用soft.gz获取ids注释文件
require(GEOquery)
GPL_data <- getGEO(filename = 'GPL19978_family.soft.gz', AnnotGPL = F)
ids_data <- Table(GPL_data)
ids_data[1:4,1:4];dim(ids_data)
ids=ids_data
colnames(ids)
# [1] "ID" "circRNA" "SPOT_ID" "PROBE_TYPE" "BUILD" "SEQUENCE"
ids=ids[,c("ID","circRNA")]
colnames(ids)=c("probe_id","symbol")
head(ids)
### 获取分组
group=factor(c(rep("ctrl",7),rep("HCC",7)),levels = c("ctrl","HCC"));group
# [1] ctrl ctrl ctrl ctrl ctrl ctrl ctrl HCC HCC HCC HCC HCC HCC HCC
# Levels: ctrl HCC
#### limma差异分析
design=model.matrix(~group)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
#去重复
deg$logFC_abs=abs(deg$logFC)
deg=deg[order(deg$ID,deg$logFC_abs,decreasing = T),] ##最终目的,保留最大的logFC值对应的探针
deg<- deg[!duplicated(deg$ID),] ## 去除重复ID
### 加change列
library(dplyr)
logFC_t=1
p_t = 0.05
k1 = (deg$adj.P.Val < p_t)&(deg$logFC < -logFC_t)
k2 = (deg$adj.P.Val < p_t)&(deg$logFC > logFC_t)
deg <- mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)
### 获取目标基因
rna_up=deg[deg$change == "up",]; rna_down=deg[deg$change == "down",]
top10_up=head(rna_up[order(rna_up$logFC_abs,decreasing = T),],10)
top10_down=head(rna_down[order(rna_down$logFC_abs,decreasing = T),],10)
top10_all=rbind(top10_up,top10_down)
#### 完成注释
heatmap_exp=exp[top10_all$ID,]
heatmap_exp=as.data.frame(heatmap_exp)
probe_id=rownames(heatmap_exp)
heatmap_exp$probe_id=probe_id
heatmap_exp=merge(ids,heatmap_exp,by="probe_id")
rownames(heatmap_exp)=heatmap_exp$symbol
heatmap_exp=heatmap_exp[,c(-1,-2)]
#### 绘制差异基因热图
library(pheatmap)
annotation_col=data.frame(group=group)
rownames(annotation_col)=colnames(heatmap_exp)
heatmap_plot <- pheatmap(heatmap_exp,show_colnames =F,
scale = "row",
#cluster_cols = F,
annotation_col=annotation_col,
breaks = seq(-3,3,length.out = 100))
.....没有可能是一样的吧.. 连样本数量都不一样
emm,作者是不是真的上传了数据?
那正确的GSE号是多少呢?
按照GSE97332的GPL19978号搜索使用同个芯片获取的数据集。
搜索breast关键词,看到了如下结果
往下直接有3个数据集与乳腺癌相关
回到文章查看样本信息
应该有10个样本做了这个芯片实验,发现上面所显示的样本都对不上。
顺带看一眼芯片的参数,
再看一眼GPL19978的注释文件
文章说有5639个circular RNAs的探针位置,但GPL19978的注释文件一共5490行。
芯片号应该也不对。
但在搜索的过程中,看见了V1这个版本序号
觉得应该有V2版本 于是顺着找到了GPL21825
去arraystar官网找这个芯片的数据,查了一下文章发表目录,确实存在这篇文章。
看来确实使用了arraystar的系列芯片。
顺着GPL21825这条线索继续,搜索breast
搜到了4个
逐一查看,最后发现了有一个GSE182471
顺着样本数量看下来
在GSE182471的时候
样本和文章的描述能对上诶~
网页往上翻,看作者的名字和机构
再贴一下作者名单吧
嗯,一模一样,应该就是这个数据集了。
就是这个数据上传时间略晚。
文章2018年发表,数据2021年上传。
拿新发现的数据跑了一下流程。
......也不是很对得上....
作者分析可能用的不是limma包,所以前10的差异基因才不一样吧。
懒得再验证了。
看到这篇文章发表和数据上传的时间差,找数据集真的需要学会放弃。