有些文章的原始数据找不到就算了吧，真不一定有

#先简单加载一堆包
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))