本系列文章是简单的 Rna-seq 分析全流程，主要是代码和结果展示。
1. 简介
分析数据可以从 NCBI，GEO 等网站上下载，进行数据挖掘，本次使用数据从 GEO 下载。
GEO（Gene Expression Omnibus）高通量基因表达数据库
网址：GEO官网
其下两个子数据库：
* GEO DataSets：收录整个实验的数据集，一般可以从中挖掘。
* GEO Profiles：用以储存基因在数据集中的表达谱。
2.数据信息
下载的数据为：
训练集：* GSE104954-GPL22945_series_matrix.txt.gz（储存表达量文件）
* GPL22945.soft.gz （储存样本信息文件）
验证集：* GSE108113_series_matrix.txt.gz（储存表达量文件）
* GPL19983.soft.gz （储存样本信息文件）
以上文件都需要用，为之后数据筛选分析做准备。
3. 思路
1. 通过 R 包 “GEOquery” 将上述文件的基因表达量信息，样本信息提取到数据集中
2. 根据样本信息将基因表达量分组，例如本次研究目标是 抗中性粒细胞胞浆抗体（Antineutrophil cytoplasmic antibody，ANCA）相关性血管炎（ANCA - associated vasculitis，AAV）之后简写为 AAV。就可以通过样本信息将AAV相关的样本和正常人的样本提取出来分组，并添加对应分组列如：AAV和Control。
3. 确保之前的表达量数据即的样本顺序与样本信息的样本顺序一致。
4. 通过 R 包 “GEOquery” 提取该实验集获取探针号（ID）与基因 Symbol 和 ENTREZ_GENE_ID 的对应关系。
5. 将样本信息和对应的分组单独提到一个 dataframe 里，后续操作会经常用到分组信息。
6. 保存之前处理过的数据文件，方便后续分析。
4. 代码
##---- 1. 已下载数据读取 ----
suppressMessages(library(GEOquery))
file_path <- "GSE104954-GPL22945_series_matrix.txt.gz"
eSet <- getGEO(filename = file_path)
exp <- exprs(eSet)
pd <- pData(eSet)
gpl <- eSet@annotation
eSet_v <- getGEO(filename = "GSE108113_series_matrix.txt.gz")
exp_v <- exprs(eSet_v)
pd_v <- pData(eSet_v)
fd_v <- fData(eSet_v)
gpl_v <- eSet_v@annotation
###---- 2. 提取临床信息 ----
pd_AAV = subset(pd, pd['diagnosis:ch1'] == 'ANCA-associated vasculitis')
pd_C = pd[is.na(pd$`diagnosis:ch1`), ]
pd_train = rbind(pd_AAV, pd_C)
pd_v_AAV0 = pd_v[grep("ANCA Associated Vasculitis", pd_v$title),]
pd_v_AAV = pd_v_AAV0[grep("tissue: Tubulointerstitium from kidney biopsy", pd_v_AAV0$characteristics_ch1),]
pd_v_C = pd_v[grep("Healthy Living Donor", pd_v$title),]
pd_verify = rbind(pd_v_AAV, pd_v_C)
###---- 3. 调整 pd 的行名顺序与 exp 列名完全一致 ----
exp_train = exp[, c(rownames(pd_train))]
exp_verify = exp_v[, c(rownames(pd_verify))]
p = identical(rownames(pd_train),colnames(exp_train));p
if(!p) exp_train = exp_train[,match(rownames(pd_train),colnames(exp_train))]
p1 = identical(rownames(pd_verify),colnames(exp_verify));p1
##---- 4.ids 平台注释 ------------
GPL22945 <-getGEO(gpl,destdir =".")
GPL22945_anno <- Table(GPL22945)
probe_symbols_ENid = GPL22945_anno[c('ID','Symbol','ENTREZ_GENE_ID')] # 获取探针号与基因Symbol的对应关系
exp_train = merge(rownames_to_column(as.data.frame(exp_train),'ID'), probe_symbols_ENid, by = "ID", all.x = T)
exp_verify = merge(rownames_to_column(as.data.frame(exp_verify),'ID'), fd_v, by = "ID", all.x = T)
exp_verify = merge(exp_verify, probe_symbols_ENid, by = "ENTREZ_GENE_ID", all.x = T)
# lncRNA
options(stringsAsFactors = F)
lncRNA_genecode <- read.table("anno.txt")
tmp<-sapply(strsplit(lncRNA_genecode$V1,'.',fixed = T),function(x)x[1])
lncRNA_genecode$GENEID<-tmp
colnames(lncRNA_genecode)[2]<-'GENEBIOTYPE'
colnames(lncRNA_genecode)[3]<-'gene_name'
tmp1<-lncRNA_genecode[lncRNA_genecode$gene_name %in% exp_train$Symbol,]
lncRNA<-exp_train[match(tmp1$gene_name,exp_train$Symbol),]
filtered_exp_train <- exp_train[!exp_train$Symbol %in% tmp1$gene_name, ]
filtered_exp_train <- filtered_exp_train[, !colnames(filtered_exp_train) %in% c("ID","ENTREZ_GENE_ID")]
rownames(filtered_exp_train) <- NULL
filtered_exp_train <- column_to_rownames(filtered_exp_train, "Symbol")
exp_train_old <-exp_train
exp_train<-filtered_exp_train
##---- 5. 分组 -----
train_group_AVV = data.frame(Sample = rownames(pd_AAV), label = "AAV")
train_group_C = data.frame(Sample = rownames(pd_C), label = "C")
train_group = rbind(train_group_AVV,train_group_C)
valid_group_AVV = data.frame(Sample = rownames(pd_v_AAV), label = "AAV")
valid_group_C = data.frame(Sample = rownames(pd_v_C), label = "C")
valid_group <- rbind(valid_group_AVV,valid_group_C)
###---- 6. save ----
write.csv(pd, file = "pd_train.csv")
write.csv(exp_train, file = "exp_train.csv")
write.csv(probe_symbols_ENid, file = "probe.csv")
write.csv(train_group, "train_group.csv")
write.csv(pd_v, file = "pd_verify.csv")
write.csv(exp_verify, file = "exp_verify.csv")
write.csv(valid_group, "valid_group.csv")
