生物医学数据库：
核酸，变异，表达，蛋白结构，功能，通路，表型，疾病，诊疗，药物
实验原始数据(公共数据)：
1KG：目标是发现人群中频率大于1%的变异位点。下载
GEO：GEO数据库全称GENE EXPRESSION OMNIBUS，是由美国国立生物技术信息中心NCBI创建并维护的基因表达数据库。它创建于2000年，收录了世界各国研究机构提交的高通量基因表达数据。GEO数据库是一个储存芯片、二代测序以及其他高通量测序数据的一个数据库。
GEO芯片数据分析
获得了芯片的GSE编号，三种不同的方法下载。
- 下载芯片的原始数据，不推荐
- 下载表达矩阵（series matrix）
rt <- read.table("GSE39582_series_matrix.txt.gz", sep = 't', header = T, comment.char = '!')例如在芯片中，包含了54675个基因探针，586个患者。 - 使用R的GEOquery包里面的getGEO函数直接读取进来（推荐）
library(GEOquery)
gset <- getGEO("GSE42872", destdir = ".")
考虑到网速问题，我们可以对参数进行设置，选择不下载平台的注释文件
gset <- getGEO("GSE42872", destdir = ".", AnnotGPL = F, getGPL = F)
3.1 通过pData函数获取分组信息，临床信息
pdata <- pData(gset[[1]])table(pdata$source_name_ch1)
3.2 通过exprs函数获取表达矩阵并校正
整理完临床信息提取对应的表达数据。对于表达数据除了下载Series Matrix后直接使用read.table函数读取，也可直接从GEOquery下载得到的变量gset中提取。
exp <- exprs(gset[[1]])
boxplot(exp,outline=FALSE, notch=T,col=group_list, las=2)
使用exprs函数可以从gset[[1]]提取表达信息；同时，我们可以使用boxplot函数先简单看一下整体样本表达情况。
避免误差，正式分析前需要对其进行人工校正一下。这里我们用limma包内置的一个函数，normalizeBetweenArrays函数。
library(limma)
exp=normalizeBetweenArrays(exp)
boxplot(exp,outline=FALSE, notch=T,col=group_list, las=2)
range(exp)
使用range函数查看表达数据exp的取值范围；一般而言范围在20以内的表达值基本已经经过了log转换。
将探针的id转换成为基因的Gene symbol
a. R包转换
index = gset[[1]]@annotation
通过提取列表gset[[1]]中的注释信息，可以看到，该芯片使用的是我们最常见的平台，GPL570。GPL570应的R包是hgu133plus2.db包
ids <- toTable(hgu133plus2SYMBOL)
head(ids)
length(unique(ids$symbol))
table(sort(table(ids$symbol)))
经过去重复一共存在20174个不同的Gene symbol，且部分基因存在多条探针的对应关系。
library(tidyverse) exp <- as.data.frame(exp) exp <- exp %>% mutate(probe_id=rownames(exp)) %>% inner_join(ids,by="probe_id") %>% select(probe_id, symbol, everything) exp <- exp[!duplicated(exp$symbol),] rownames(exp) <- exp$symbol exp <- exp[,-(1:2)]
经过id匹配，并去重复，最终得到了20174个基因的表达结果；
b. 使用soft文件注释
c. 手工注释
4. PCA分析
表达矩阵整理完成。在正式差异分析之前，首先可做一个PCA分析，整体水平查看正常和肿瘤样品是否存在显著的差异。
library(FactoMineR)
library(factoextra)
dat= as.data.frame(t(exp))
dat.pca <- PCA(dat, graph = FALSE)
pca_plot <- fviz_pca_ind(dat.pca, geom.ind = "point", col.ind = group_list, palette = c( "#00AFBB", "#E7B800"), addEllipses = TRUE, legend.title = "Groups" )
- 差异分析
芯片数据的差异分析，一般使用limma包来进行。差异分析的输入文件，主要是两个，一是整理好的表达矩阵，行名为基因名，列名为样本名；二是分组信息（group list）
library(limma)
design=model.matrix(~group_list)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
使用topTable函数提取所有基因的差异分析。deg包含了6块的内容，包括常见的logFC值，P.Value，adj.P.Val等等。
logFC=1.5
P.Value = 0.05
type1 = (deg $P.Value < P.Value)&(deg $logFC< -logFC) type2 = (deg $P.Value < P.Value)&(deg $logFC> logFC)
deg $change= ifelse(type1, "down",ifelse(type2, "up", "stable"))
table(deg $change)
根据设定的阈值，｜logFC｜>1.5和P值<0.05，将显著差异表达的基因进行分组。结果显示：显著下调的基因有427个，显著上调的基因有242个。
6. 可视化分析【差异分析结果，火山图和热图是最常见展示方式。】
a. 基于ggpubr包绘制火山图
#install.packages( "ggpubr")
#install.packages( "ggthemes")
library(ggpubr)
library(ggthemes)
deg$logP <- - log10(deg$P.Value)
ggscatter(deg, x = "logFC", y = "logP", color = "change", palette = c( "blue", "black", "red"), size = 1) + theme_base +geom_hline(yintercept = - log10(P.Value), linetype = "dashed") + geom_vline(xintercept = c(-logFC, logFC), linetype = "dashed")
进一步给火山图添加标签，把显著上调和显著下调基因中前5名的基因名进行标注；
添加标签deg $Label= ""#新加一列labeldeg <- deg[order(deg $P.Value), ]
up.genes <- head(rownames(deg)[ which(deg $change== "up")], 5)
down.genes <- head(rownames(deg)[ which(deg $change== "down")], 5)
#将up.genes和down.genes合并，并加入到Label中deg.top5.genes <- c(as.character(up.genes), as.character(down.genes))
deg $Label[match(deg.top5.genes, rownames(deg))] <- deg.top5.genes
ggscatter(deg, x = "logFC", y = "logP", color = "change", palette = c( "blue", "black", "red"), size = 1, label = deg$ Label, font.label = 8, repel = T, xlab = "log2FoldChange", ylab = "-log10(P-value)") + theme_base +geom_hline(yintercept = -log10( P. Value), linetype = "dashed") + geom_vline(xintercept = c(-logFC, logFC), linetype = "dashed")
b. 基于ggplot2包绘制火山图?再调整
library(ggplot2)
p < -ggplot( data= deg,aes( x= logFC,y= -log10(P.Value)))+ geom_point( alpha= 0.4,size= 3.5,aes( color= change))+ ylab(" -log10( Pvalue)")+ scale_color_manual( values= c(" blue", " grey"," red"))+ geom_vline( xintercept= c(-logFC,logFC),lty= 4,col= "black", lwd= 0.8)+ geom_hline( yintercept= -log10(P.Value),lty= 4,col= "black", lwd= 0.8)+ theme_bwp)
添加标签的操作；
x1 = deg %>% filter(change == "up") %>% head(5)x2 = deg %>% filter(change == "down") %>% head(5)label = rbind(x1,x2)volcano_plot <- p +geom_point(size = 3, shape = 1, data = label) +ggrepel::geom_label_repel(data = label, aes(label = rownames(label)), color= "black")
c. 热图
cg= rownames(deg)[deg $change!= “stable”] diff=exp[cg,]
提取差异表达基因的表达情况；
annotation_col= data.frame(group=group_list)
rownames(annotation_col)= colnames(diff) pheatmap(diff,annotation_col= annotation_col,scale= "row",show_rownames= F,show_colnames= F,color= colorRampPalette(c("navy", "white", "red"))(50),fontsize= 10,fontsize_row= 3,fontsize_col= 3)
GEO数据库芯片数据的下载，probe id的转换，差异分析已经基本完成.针对这些差异基因的常规分析，包括GO分析，KEGG分析，GSEA分析，蛋白-蛋白互作网络（Protein-protein interaction，PPI）
- TCGA数据库
TCGA收录的了全面的癌症基因组数据，包括突变，拷贝数变异，mRNA表达，miRNA表达，甲基化数据等
TCGA上储存数据分为三个级别，level-1为原始测序数据（fasta,fastq等），level-2为比对的bam格式文件，（这两级别为controlled-access，需要向TCGA申请使用权限），level-3为经过处理及标准化的数据(三级也分为controlled-access和open-access)
level-1/2 数据的下载需要向TCGA申请使用权限，并且由于文件较大，推荐使用官方提供的小小软件：gdc.cancer.gov/access-data/gdc-data-transfer-tool
最常用的是level-3数据，一般文件较小直接在网页上下载就可以。
目前主要有两个网站可以下载TCGA level-3的数据:
TCGA官网的data-portal: portal.gdc.cancer.gov
优点：数据最全，更新最快
缺点：每个样本的数据都单独储存在一个文件中，如果要下载RNA表达量数据的话，可能同一种癌症需要下载好几百个文件，并且需要排队下载，有时候很慢很慢很慢
Firehose服务器：gdac.broadinstitute.org
（这里的数据也来源于portal.gdc.cancer.gov，经过了简单的处理）
优点：经过了简单的合并，将每种癌症相同类型的数据合并到了一个文件中（例如443个胃癌样本的RNA表达量数据都合并到了一个文件中，非常适合用R进行后续的分析）
level-3的数据是仍需要一定的分析能力来提取感兴趣的信息
如果你仅需要看感兴趣的基因在某种癌症中的突变谱，表达量，或者甲基化情况，那么以下三个在线可视化网站：
1.>c-Bioportal: www.cbioportal.org
整合和简化了包括TCGA，ICGC以及GEO等多个癌症基因组的内容，提供友好可视化的界面，可供下载。
主要展示基因的somatic 突变谱，拷贝数变化，mRNA&miRNA表达量变化，DNA甲基化以及蛋白质表达的情况，并结合患者的临床资料，展示了KM生存曲线。
2.>OncoLnc: www.oncolnc.org
这是一个整合了TCGA的各种RNA数据和患者临床数据，提供生存分析的网站，灰常简单好用。
3.>MEXPRESS:mexpress.be/about
整合了TCGA中的DNA甲基化，表达量及临床数据，主要用来探索甲基化，基因表达和临床表型之间的关联，看界面也很友好，但我没怎么用过%>_<%。
BIGDDBJEMBL-EBI
扫一扫
495
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
