表达芯片分析[第一讲：GEO，表达芯片与R]

Andyhigh

先看看这一系列文章整体思路，文章见QQ群：

下面是一个中文的例子：
paper：http://d.wanfangdata.com.cn/Thesis/Y1504435
博士论文：http://www.doc88.com/p-116632873539.html

[研究目的] 分离在胃癌肝转移过程中起关键作用的关键基因(Hubs)，建立能预测胃癌发生及其转移的新的生物学标志，为未来开展胃癌分子流病学研究奠定基础。

[研究方法] 2006年1月到2009年2月收集67例胃癌病人的手术切除标本，所有标本均于术后2小时内取回、分装和冻存；

选择3例患者的胃癌肝转移及配对原发组织进行cDNA芯片检测，选择在3张芯片中基因表达变化倍数均≥1.6倍的基因作为差异表达基因；

利用MILANO、Chilibot、UniHI和Endeavour等多种生物信息学软件工具，依次对差异表达基因进行系统的文献挖掘、蛋白质网络和功能相似性分析，最后预测胃癌肝转移相关的新Hub基因(相互作用基因数≥30)；

根据预测结采选择兴趣基因，采用定量RT-PCR技术，检测这些兴趣基因在另外8对胃癌肝转移及配对原发组织、57对胃癌原发灶及配对癌旁组织间的相对表达变化，并分析这些基因的相对表达变化与疾病分期之间的关系；

采用免疫组化技术，检测兴趣基因所对应蛋白在11例胃癌肝转移、26例胃癌原位及26例胃癌癌旁组织中的表达模式，并将蛋白表达变化与转录水平表达变化进行比较。

[研究结果] 芯片分析共确定272个差异表达基因，通路分析显示差异表达基因在13条KEGG通路上显著富集，其中8条通路已证明与肿瘤转移相关；

文献注释结果显示272个差异表达基因中包含118个已知癌症相关基因和68个已知肿瘤转移相关基因；已知转移相关基因在118个癌症相关基因中的比例明显高于已知转移相关基因在Sanger中心的382个癌症普查基因中的比例(P＜0.001)，提示本研究所获得的272个差异表达基因显著富集了转移相关基因；

蛋白质网络分析发现所获得的差异表达基因中存在一个由63个蛋白形成的亚网络，该亚网络包含60个相互作用、29个已知转移相关基因和37个Hub基因。

统计分析显示，该网络不仅富集了已知转移相关基因(P＜0.05)，而且富集了Hub基因(P＜0.05)；利用Endeavour软件进行基因功能相似性分析，发现45个未报道与转移相关的基因与我们所设定的已知转移相关基因群的功能最为相似，其中8个基因为蛋白质亚网络中的Hub基因：NR3C1、NR4A2、HNRPA1、PSMB3，FBLN2、DARS、XAB2，CD8A。

根据预测结果和研究兴趣，9个基因(NR3C1、NR4A2、HNRPA1、XAB2、HSP90AA1、CCNEl、RPL17、FKBP1A和XAB2)被选择在临床样本中验证，定量RT-PCR分析结果显示：NR4A2在8例同时性胃癌肝转移样本中的表达水平显著低于其配对原发灶(P=0.001)，而HSP90AA1的表达水平显著高于其配对原发灶(P=0.029)；NR4A2、NR3C1、ARF3、XAB2以及NR4A2的两种替换剪切模式(SP8和SP-novel)在57例胃癌原发样本中的表达水平显著低于其配对癌旁组织(P＜0.001)，而CCNE1的表达水平显著高于其配对癌旁组织(P=0.001)。 HSP90AA1(P=0.043)和NR4A2(P=0.003)在胃癌肝转移组织中的蛋白表达水平分别高于和低于原发组织的表达水平，CCNE1(P=0.000)和NR3C1(P=0.005)在胃癌原发灶组织中的蛋白表达水平分别高于和低于癌旁正常组织的表达水平，上述结果与定量RT-PCR结果一致；HSP90AA1在胃癌原发灶中的蛋白表达水平高于癌旁正常组织(P=0.003)，而定量RT-PCR未发现差异。在定量RT-PCR分析中显示：癌旁正常组织中NR4A2基因的表达水平显著高于原位组织，但免疫组化显示NR4A2蛋白在癌旁正常组织间质细胞中的表达水平高于原位肿瘤间质细胞细胞中的表达水平，而在癌旁胃腺组织上皮细胞中的表达水平低于原位胃癌上皮组织中的表达水平。

[研究结论] NR4A2可能是胃癌肝转移的诊断和预后标志，NR4A2可能参与胃癌上皮细胞的上皮-间质转化(EMT)过程；NR3C1、ARF3和XAB2与胃癌发生相关，可能是胃癌发生中的新型Hub基因；HSP90AA1可能是胃癌肝转移的预后标志，而CCNE1是胃癌发生的诊断标志，但尚需进一步研究。

其它：


GEO介绍：http://www.bio-info-trainee.com/bioconductor_China/software/GEOquery.html
所有bioconductor支持的芯片平台对应关系：通过bioconductor包来获取所有的芯片探针与gene的对应关系
gpl           organism                  bioc_package    https://www.bioconductor.org/packages/devel/data/annotation/
芯片探针注释基因ID或者symbol，并对每个基因挑选最大表达量探针
凡是在GEO里面有GPL信息的，在ensembl的arrayexpress里面都会有对应的！ https://www.ebi.ac.uk/arrayexpress/files/A-AFFY-104/A-AFFY-104.adf.txt
常用的在R的bioconductor里面有对应，但芯片厂商是肯定有的


没有必要用R包GEOquery  ,拼接URL，拼接URL，拼接URL

• 使用GEOmetadb包来获取对应GEO数据的实验信息
• 从GEO数据库下载矩阵数据-可以直接进行下游分析
  

2011-pnas-gene-expression-signature-NSCLC：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3084093/
涉及到3种数据，我们只看affymetrix的表达数据
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE22874


原始数据：15个非小细胞肺癌患者中提取出肺癌成纤维细胞和正常成纤维细胞样本，并做了Affymetrix平台下的芯片数据。
差异分析：6种方法(limma/SAM/t-test/~~~)分析肺癌成纤维细胞和正常成纤维细胞的基因表达水平差异，筛选出肺癌成纤维细胞中高表达或者低表达基因（通常称之为差异表达基因）
富集分析：在得到高表达和低表达基因后，对这些基因作了Gene ontology和KEGG通路富集分析
GSEA分析：运行软件即可
生存分析：需要收集非小细胞肺癌病人的临床数据
PPI网络构建：在STRING官网，或者用stringDB包
BioNet分析：寻找maximal-scoring subgraph (MSS)

GSE63067
GE37031
GSE20950
GSE26339
GSE15773
GSE30122
GSE23343
GSE21785
GSE1009
GSE15653
GSE25724
GSE20966

表达芯片并非只有mRNA的表达量
lncRNA : http://pubmedcentralcanada.ca/pmcc/articles/PMC3691033/ GSE33335  Affymetrix Human Exon 1.0 ST Array [transcript (gene) version
Expression data from gastric tissues: Cancer Samples vs. Matched Adjacent Noncancerous Samples

R包终极解决方案！   

还有两个脚本，一个是从GEO数据库里面获取所有的GPL信息：
#~/bin/bash
# This script will update data from GEO

DIR=~
DIR=`pwd`
{
        now=`date +"%m_%d_%Y"`
        mv $DIR/platform.csv $DIR/backup-platform.csv &&
        echo "Backup created.."

} || {

        echo "No local platform data found.. one will be created"


}       
echo "Update GEO data..."
data=$(curl --silent "http://www.ncbi.nlm.nih.gov/geo/browse/?view=platforms" | grep "total_count")
count=$(echo "$data" | grep -Eo '[0-9]*')
count=`python -c "from math import ceil; print int(ceil($count/5000.0))"`

echo "Getting GEO Data..."
for ((i=1; i <= count; i++)); do

if [ "$i != 1" ]
then
        wget --content-disposition "http://www.ncbi.nlm.nih.gov/geo/browse/?view=platforms&zsort=date&mode=csv&page=$i&display=5000" -q -O ->> $DIR/platform.csv
else
        wget --content-disposition "http://www.ncbi.nlm.nih.gov/geo/browse/?view=platforms&zsort=date&mode=csv&page=$i&display=5000" -q -a  ->> $DIR/platform.csv       

fi

done

echo "Finished update...data saved to -> platform.csv"

一个是从GEO数据库获取所有的study ID信息;
#~/bin/bash
# This script will update data from GEO

#DIR=/home/user
DIR=`pwd`
{

        mv $DIR/series.csv $DIR/backup-series.csv &&
        echo "Backup created.."

} || {

        echo "No local series data found.. one will be created"
}       

echo "Update GEO data..."
data=$(curl --silent "http://www.ncbi.nlm.nih.gov/geo/browse/" | grep "total_count")
count=$(echo "$data" | grep -Eo '[0-9]*')
count=`python -c "from math import ceil; print int(ceil($count/5000.0))"`

echo "Getting GEO Data..."
echo "Using count = $count"

for ((i=1; i <= count; i++)); do

if [ "$i != 1" ]
then
        wget --content-disposition "http://www.ncbi.nlm.nih.gov/geo/browse/?view=series&zsort=date&mode=csv&page=$i&display=5000" -q -O ->> $DIR/series.csv
else
        wget --content-disposition "http://www.ncbi.nlm.nih.gov/geo/browse/?view=series&zsort=date&mode=csv&page=$i&display=5000" -q -a  ->> $DIR/series.csv       

fi

done

echo "Finished update...data saved to -> series.csv"

Andyhigh

后面的脚本看不懂！！！

果子

你最后标注一下，这是jimmy的东西，你要分析要注明，我们只能分享自己的东西。

Andyhigh

果子 发表于 2016-12-19 19:32
你最后标注一下，这是jimmy的东西，你要分析要注明，我们只能分享自己的东西。 ...

嗯，当时是一个下午Jimmy叫我帮他博客上的内容转过来的，说晚上讲课要用，比较急，就没有注意到这些。
下次注意！！！

y461650833y

这不是群主的讲义