Monthly Archives: 1月 2016

一 31

画基因的外显子覆盖度图

Posted on by
一般情况下，我们得到了测序reads在基因组的比对情况文件bam格式的，里面的信息非常多，如果我想特定的查看某个基因的情况，那么我们可以选择IGV等可视化工具，但它并不是万能的，因为即使是一个基因，它也会有多个转录本，多个外显子。
所以，我们可以画它的外显子覆盖图，如下：横坐标是外显子的长度，纵坐标是测序深度，每一个小图都是一个外显子
 DMD.NM_000109
根据这个图，我们就可以很明显的看出，DMD基因NM_000109转录本的1，10-17号外显子缺失，用IGV一个个的看这些外显子区域，是同样的结果！可能是芯片捕获不到，也可能是样本本身变异，造成的大片段缺失。但是这个图的信息就非常有用！
那么，我们该如何画这样的图呢？
首先，我们需要找到需要探究的基因的全部转录信息，及外显子信息！
在hg19_refGene.txt里面会有，在UCSC里面可以下载，新手可能会比较麻烦，实在不行你去annovar的目录也可以找到！
1
那么，我们根据这个信息，就可以判断该基因的起始终止位点啦
然后用samtools的depth命令去找这个基因的全部片段的测序深度信息
最后再格式化成下面的三列数据
第一列是该外显子的坐标，从1到该外显子的成都
第二列是该外显子在该坐标的测序深度，通过samtools的depth命令得到
最后一列是该外显子的标记，从exon:79一直倒推到exon:1，因为该基因在染色体的负链，所以外显子顺序是反着的！
1 84 exon:79
2 84 exon:79
3 84 exon:79
4 85 exon:79
5 85 exon:79
6 86 exon:79
7 85 exon:79
8 87 exon:79
9 89 exon:79
10 91 exon:79
11 92 exon:79
12 95 exon:79
13 96 exon:79
14 96 exon:79
15 99 exon:79
16 99 exon:79
17 97 exon:79
最后根据这个txt文档，用R语言，很容易就画出上面那样的图片了！
这里面的信息量还是蛮大的！
APOB.NM_000384

 

一 28

使用可视化工具MutationMapper来看看基因上面突变的分布

Posted on by
how to generate lollipop diagrams ?
这个工具是网页版的，不需要下载，只需上传两列数据即可，就是基因上面的第几个氨基酸是如何突变的这样的信息
网页就可以把这些信息画在基因上面，而且注释到domain和cosmic数据库，挺好用的！
Hugo_Symbol HUGO symbol for the gene TP53
Protein_Change Amino acid change V600E
结果也很容易理解！
我输入了4个基因，网页就会对每个基因画一个图，其实ABCA1这个基因有6个突变，就都画在该基因上面了，基因突变分布在基因的哪个domain也看得一清二楚！
1
我的输入文件是这样的！
2
非常好用：
  • Support mutation data with annotated protein effects
  • Mutation diagram/lollipop view
  • Mutation table view
  • 3D structure view if available
而且pfam数据库本身也提供了这个功能：

Pfam provides an online tool to not only generate the domain information in JSON format, but to draw the lollipop diagram using javascript as well. They have more information here: http://pfam.xfam.org/help#tabview=tab9

IMHO, not as pretty as cBioPortal's but it gets you close to a solution.

EDIT / SHAMELESS PLUG: After seeing the data available and how easy it'd be, I made my own quick tool to fetch the data and draw the diagram for me in a style similar to cBioPortal - feel free to fork it and add features: https://github.com/pbnjay/lollipops

Example output (w/ labels per the comments)

3

甚至还有高手用D3.JS写了一个能实现同样需求的模块

We found ourselves in the same need, we wanted such a plot (JavaScript). Thus, I add our solution, Mutations Needle Plot. The library creates an SVG image (with D3), which then may be downloaded.

You will npm in order to be able to install & run the library.

Examples may be found in the snippets folder or also the index.html - The one displayed here below

类似的高手很多：
My colleague @SolenaLS recently asked me to write something like this: (uniprot+SVG+javascript : )http://lindenb.github.io/pages/uniprot/paintsvg.html
而且ensembl数据库也提供类似的功能
Ensemble also gives similar plot.

 

一 27

gene的symbol与entrez ID并不是绝对的一一对应的

Posted on by

很多时候，我们都无法确定到底是基于symbol来进行分析，还是基于entrez ID，当我们要进行ID转换的时候也想当然的以为它们的一一对应的， 但是最近我写了一个脚本来分析CCLE的数据的时候，发现其实有一些特例:

suppressMessages(library(org.Hs.eg.db))

all_symbol=mappedkeys(org.Hs.egSYMBOL2EG)

all_EGID =mappedkeys(org.Hs.egSYMBOL)
tmp=as.list(org.Hs.egSYMBOL2EG[all_symbol])
#tmp=as.list(org.Hs.egSYMBOL[all_EGID ])
tmp=unlist(lapply(tmp,length))
tmp=tmp[tmp>1]
as.list(org.Hs.egSYMBOL2EG[names(tmp)])
有多个entrez ID对应一个symbol的现象出现,但是没有一个symbol对应多个entrez ID的现象。而且entrez ID也会过期！
$CSNK1E
[1] "1454"      "102800317"
$HBD
[1] "3045"      "100187828"
$RNR1
[1] "4549" "6052"
$RNR2
[1] "4550" "6053"
$SFPQ
[1] "6421"   "654780"
$TEC
[1] "7006"      "100124696"
$MEMO1
[1] "7795"  "51072"
$KIR3DL3
[1] "115653"    "100133046"
$MMD2
[1] "221938"    "100505381"
$`LSAMP-AS1`
[1] "100506708" "101926903"
通过下面的链接可以看到具体情况
一 22

用RankComp的思想来做差异基因分析

Posted on by
是福建医科大学的学者开发的，
文章里面详细讲解了他们的这个差异分析的统计学原理
大意就是找到同一组织的normal样本的表达量数据，几百个，这样就可以分析2万基因之间的互相配对，检测表达量是否在几百个样本里面稳定的不一样！

我现在还不是很确定这个方法，只是试一试，欢迎与我交流对该方法的讨论！

文章是：

Wang H, Sun Q, Zhao W, et al. Individual-level analysis of differential expression of genes and pathways for personalized medicine[J]. Bioinformatics, 2014: btu522.

他们把它写成了一个R包，可以下载使用，但是必须用R2.15.2版本，我用了一下，不好用！

We can download the R code for in http://bioinformatics.oxfordjournals.org/content/31/1/62/suppl/DC1

他们这个程序真心不好用，但是很容易看懂算法，可以自己用R语言写一个来实现同样的过程！

比如A基因在几百个样本里面表达量都是3左右，而B基因都是5左右，而且满足99%的A表达量高于B，那么这就是一个稳定的基因对！
一般2万基因之间可以配成2亿个基因对，其中稳定的大概有10%~40%
然后我们对每个疾病样本都可以进行检验，看看这样稳定的基因对是否被改变！
比如，我们拿到一个疾病样本的2万个基因的表达量，我们挑取一个，如果它有100个稳定的up的基因对，100个稳定的down的基因对
那么，我看看这些基因对是否被改变了，如果这样还有70基因对在该疾病样本里面仍然是up，60个是down，那么我用Fisher精确检验的结果是
4
这个基因在该疾病样本，相对于normal pool并没有差异表达！当然检验得到的P值最后可以做FDR校验。
依次这样，把所有的gene都分析完，就知道这个样本有哪些差异的gene了。
介绍完原理，我们拿一个具体的例子来看看吧：
首先我们下载一个2008年的一个人的肝脏表达数据（Gene Expression in Human Liver），都是正常组织，共427个样本。
不过这个芯片比较小众，是默克医药公司定制化的， 需要下载探针对应gene的文件！
我们读取GSE9588这个数据，得到表达矩阵，然后计算rank矩阵，然后计算得到comp矩阵
5
> table(rank_comp)
rank_comp
     down        no        up
    58479 465752098     58479
>
不知道为什么这个数据，stable的那么少，不知道是不是出了什么问题！
其实我的程序都是对的了，只是因为这个数据集已经不是纯粹的表达量数据了，而是这427个样本的数据都减去了某个样本的表达量。
这样每个个体的基因之间的表达量排序就会被干扰，导致得到的稳定基因对非常少！！！
但是，我后来下载了GTEx的表达数据，拿那里面的normal组织样本表达量来做，可以得到非常多的稳定基因对。
实际代码大概是：
得到正常组织的表达矩阵：
然后计算表达矩阵的rank，得到各个样本自己的基因排序情况，得到排序矩阵！
处理排序矩阵，每个基因对之间都算一下是否稳定，得到稳定性描述矩阵！
然后根据每个疾病个体的基因表达情况，来循环每个基因， 看看该基因是否差异！

 

一 22

用TCGA数据做cox生存分析的风险因子（比例风险模型）

Posted on by
再次强调一下，R里面实现生存分析非常简单！
用my.surv <- surv(OS_MONTHS,OS_STATUS=='DECEASED')构建生存曲线。
用kmfit2 <- survfit(my.surv~TUMOR_STAGE_2009)来做某一个因子的KM生存曲线。用 survdiff(my.surv~type, data=dat)来看看这个因子的不同水平是否有显著差异，其中默认用是的logrank test 方法。
用coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung) 来检测自己感兴趣的因子是否受其它因子(age,gender等等)的影响。

Continue reading →

一 14

R语言画网络图三部曲之sna

Posted on by

如果只是画网络图，那么只需要把所有的点，按照算好的坐标画出来，然后把所有的连线也画出即可！

其中算法就是，点的坐标该如何确定？Two of the most prominente algorithms (Fruchterman & Reingold’s force-directed placement algorithm and Kamada-Kawai’s)
有一个networkD3的包可以直接画图，但是跳过了确定点的坐标这个步骤，我重新找了一包，可以做到！
作者只是用sna包来得到数据，其实用的是ggplot来画网络图！
需要熟悉network包里面的network对象的具体东西，如何自己构造一个，然后数学sna包如何计算layout即可
解读这个包，也可以自己画网络图，代码如下：
plot(plotcord)
text(x=plotcord$X1+0.2,y=plotcord$X2,labels = LETTERS[1:10])
for (i in 1:10){
  for (j in 1:10){
    if(tmp[i,j]) lines(plotcord[c language="(i,j),1"][/c],plotcord[c language="(i,j),2"][/c])
  }
}
当然，我们还没有涉及到算法，就是如何生成plotcord这个坐标矩阵的！
大家看下面这个示意图就知道网络图是怎么样画出来的了，首先我们有一些点，它们之间有联系，都存储在networData这个数据里面，是10个点，共9个连接，然后我用reshape包把它转换成连接矩阵，理论上10个点的两两相互作用应该有100条线，但是我们的数据清楚的说明只有9条，所以只有9个1，其余的0代表点之间没有关系。接下来我们用sna这个包对这个连接矩阵生成了这10个点的坐标（这个是重点），最后很简单了，把点和线画出来即可！

1

另外一个例子：
net=network(150, directed=FALSE, density=0.03)
m <- as.matrix.network.adjacency(net) # get sociomatrix  
# get coordinates from Fruchterman and Reingold's force-directed placement algorithm.
plotcord <- data.frame(gplot.layout.fruchtermanreingold(m, NULL)) 
# or get it them from Kamada-Kawai's algorithm: 
# plotcord <- data.frame(gplot.layout.kamadakawai(m, NULL)) 
colnames(plotcord) = c("X1","X2")  ###所有点的坐标，共150个点
edglist <- as.matrix.network.edgelist(net) ##所有点之间的关系-edge ##共335条线
edges <- data.frame(plotcord[edglist[,1],], plotcord[edglist[,2],]) ##两点之间的连线的具体坐标，335条线的起始终止点点坐标
原始代码如下：
library(network)
library(ggplot2)
library(sna)
library(ergm)
clipboard
大家可以试用这个代码，因为它用的ggplot，肯定比我那个简单R作图要好看多了
参考算法文献：
 

 

一 14

R语言画网络图三部曲之networkD3

Posted on by
首先，我们需要了解一下基础知识：
网络图是为了展示数据与数据之间的联系，在生物信息学领域，一般是基因直接的相互作用关系，或者通路之间的联系！
通俗点说，就是我有一些点，它们之间并不是两两相互联系，而是有部分是有连接的，那么我应该如何把这些点画在图片上面呢？因为这些都并没有X,Y坐标，只有连接关系，所以我们要根据一个理论来给它们分配坐标，这样就可以把它们画出来了，然后也可以对这些点进行连线，连完线后，就是网络图啦！！！
而给它们分配坐标的理论有点复杂，大概类似于物理里面的万有引力和洛仑磁力相结合来给它们分配一个位置，使得总体的能量最小，就是最稳定状态！而通常这个状态是逼近，而不是精确，所以我们其实对同样的数据可以画出无数个网络图，只需使得网络图合理即可！
看这个ppt：http://statmath.wu.ac.at/research/friday/resources_WS0708_SS08/igraph.pdf

基本上就都理解了
画网络图真的好复杂呀！！！
大部分人都是用D3JS来画图：http://bl.ocks.org/mbostock/2706022
一步步讲受力分析图是如何画的：https://github.com/mbostock/d3/wiki/Force-Layout

接下来， 我们直接看看R里面是如何画网络图的，我们首推一个包：networkD3/

这个包非常好用，只需要做好data，然后用它提供的几个函数即可！
重要的是熟悉输入数据是什么，还可以结合shinny包来展示数据，非常好用！！！
具体怎么安装这个R包，我就不讲了，它的github主页上面其实也有说明书，很容易看懂的
最简单的用法，就是构造一个联结矩阵，把所有的连接关系都存储起来，然后直接用一个函数就可以画图啦！
# Plot
library(networkD3)
simpleNetwork(networkData,fontSize=25)
可以看出网络图，就是把所有的点，按照算好的坐标画出来，然后把所有的连线也画出即可！
其中算法就是，点的坐标该如何确定？Two of the most prominente algorithms (Fruchterman & Reingold’s force-directed placement algorithm and Kamada-Kawai’s)

clipboard

但是这个软件有一个弊端就是，生成图片之后，并没有给出这些点的坐标，当然，这种坐标基本上也很少有人需要，可视化就够了！
如果这些点还分组了，而且连接还有权重，那么网络图就复杂一点，比如下面这个：
就不仅仅是需要提供所有的link的信息，link还多了一列，是value，而且还需提供所有的node信息，node多了一列是分组。
clipboard
当然，还有好几个图，大家可以自己慢慢用，自己体会！
sankeyNetwork
sankeyNetwork
radialNetwork
diagonalNetwork
dendroNetwork
也可以把生成的网络图直接保存成网页，动态显示！
一 14

用broad出品的软件来处理bam文件几次遇到文件头错误

Posted on by

报错如下：ERROR MESSAGE: SAM/BAM file input.marked.bam is malformed: SAM file doesn't have any read groups defined in the header.  The GATK no longer supports SAM files without read groups ！

有些人遇到的是bam的染色体顺序不一样，还有可能是染色体的名字不一样，比如>1和>chr1的区别，虽然很傻，但是遇到这样问题的还不少！
还有一些人是遇到基因组没有dict文件，也是用picard处理一下就好。

大部分人是在GATK遇到的，我是在RNA-SeQC遇到的，不过原理都是一样的。
都是因为做alignment的时候并未添加头信息，比如：
bwa samse ref.fa my.sai my.fastq > my.sam
samtools view -bS my.sam > my.bam
samtools sort my.bam my_sorted
java -jar ReordereSam.jar I=/path/my_sorted.bam O=/path/my_reordered.bam R=/path/ref.fa
通过这个代码可以得到排序好的bam，但是接下来用GATK就会报错
java -jar GenomeAnalysisTK.jar -T DepthOfCoverage -R /paht/ref.fa -I /path/aln_reordered.bam
就是因为没有头信息，group相关信息，解决方法有两种：
bwa samse -r @RG\tID:IDa\tSM:SM\tPL:Illumina ref.fa my.sai my.fastq > my.sam
java -jar AddOrReplaceReadGroups I=my.bam O=myGr.bam LB=whatever PL=illumina PU=whatever SM=whatever
一种是比对的时候就加入头信息，这个需要比对工具的支持。
第二种是用picard工具来修改bam，推荐用这个！虽然我其实并不懂这些头文件信息是干嘛的， 但是broad开发的软件就是需要！希望将来去读PHD能系统性的学习一些基础知识！

 

一 14

用RNA-SeQC得到表达矩阵RPKM值

Posted on by
这个软件不仅仅能做QC，而且可以统计各个基因的RPKM值！尤其是TCGA计划里面的都是用它算的
一、程序安装
直接在官网下载java版本软件即可使用：http://www.broadinstitute.org/cancer/cga/tools/rnaseqc/RNA-SeQC_v1.1.8.jar
但是需要下载很多注释数据
clipboard

二、输入数据

clipboard

箭头所指的文件，一个都不少，只有那个rRNA.tar我没有用， 因为这个软件有两种使用方式，我用的是第一种
三、软件使用
软件的官网给力例子，很容易学习：
RNA-SeQC can be run with or without a BWA-based rRNA level estimation mode. To run without (less accurate, but faster) use the command:
java -jar RNASeQC.jar -n 1000 -s "TestId|ThousandReads.bam|TestDesc" -t gencode.v7.annotation_goodContig.gtf -r Homo_sapiens_assembly19.fasta -o ./testReport/ -strat gc -gc gencode.v7.gc.txt 
我用的就是这个例子，这个例子需要的所有文件里面，染色体都是没有chr的，这个非常重要！！！
代码如下：
 java -jar RNA-SeQC_v1.1.8.jar  \
-n 1000 \
-s "TestId|ThousandReads.bam|TestDesc" \
-t gencode.v7.annotation_goodContig.gtf \
-r ~/ref-database/human_g1k_v37/human_g1k_v37.fasta  \
-o ./testReport/ \
-strat gc \
-gc gencode.v7.gc.txt \
To run the more accurate but slower, BWA-based method :
java -jar RNASeQC.jar -n 1000 -s "TestId|ThousandReads.bam|TestDesc" -t gencode.v7.annotation_goodContig.gtf -r Homo_sapiens_assembly19.fasta -o ./testReport/ -strat gc -gc gencode.v7.gc.txt -BWArRNA human_all_rRNA.fasta
Note: this assumes BWA is in your PATH. If this is not the case, use the -bwa flag to specify the path to BWA
四、结果解读
运行要点时间，就那个一千条reads的测试数据都搞了10分钟！
出来一大堆突变，具体解释，官网上面很详细，不过，比较重要的当然是RPKM值咯，还有QC的信息
clipboard
TCGA数据里面都会提供由RNA-SeQC软件处理得到的表达矩阵！
Expression
  • RPKM data are used as produced by RNA-SeQC.
  • Filter on >=10 individuals with >0.1 RPKM and raw read counts greater than 6.
  • Quantile normalization was performed within each tissue to bring the expression profile of each sample onto the same scale.
  • To protect from outliers, inverse quantile normalization was performed for each gene, mapping each set of expression values to a standard normal.
软件的主页是：
 
 
 
 
一 14

华盛顿大学把所有的变异数据都用自己的方法注释了一遍，然后提供下载

Posted on by
华盛顿大学把所有的变异数据都用自己的方法注释了一遍，然后提供下载：
文献是：Kircher M, Witten DM, Jain P, O'Roak BJ, Cooper GM, Shendure J. 

A general framework for estimating the relative pathogenicity of human genetic variants.
Nat Genet. 2014 Feb 2. doi: 10.1038/ng.2892.
PubMed PMID: 24487276.

文中的观点是：现在大多的变异数据注释方法都非常单一，通常是看看该位点是否保守，对蛋白功能的改变，在什么domain上面等等。
但这样是远远不够的，所以他们提出了一个新的注释方法，用他们自己的CADD方法把现存的一些公共数据库的变异位点（约86亿的位点）都注释了一下，并对每个位点进行了打分。
C scores correlate with allelic diversity, annotations of functionality, pathogenicity, disease severity, experimentally measured regulatory effects and complex trait associations, and they highly rank known pathogenic variants within individual genomes.
总之，他们的方法是无与伦比的！
所有他们已经注释好的数据下载地址是：http://cadd.gs.washington.edu/download
这些数据在很多时候非常有用，尤其是想跟自己得到的突变数据做交叉验证，或者做一下统计分析的时候！
clipboard
人的基因组才300亿个位点，他们就注释了86亿！！！
所以有三百多G的压缩包数据，我想，一般的公司或者单位都不会去用这个数据了！
一 14

蛋白质相互作用（PPI）数据库大全

Posted on by
最近遇到一个项目需要探究一个gene list里面的基因直接的联系，所以就想到了基因的产物蛋白的相互作用关系数据库，发现这些数据库好多好多！
一个比较综合的链接是：A compendium of PPI databases can be found in http://www.pathguide.org/.

里面的数据库非常多，仅仅是对于人类就有

Your search returned 207 results in 9 categories with the following search parameters:

人类的六个主要PPI是：Analysis of human interactome PPI data showing the coverage of six major primary databases (BIND, BioGRID, DIP, HPRD, IntAct, and MINT), according to the integration provided by the meta-database APID.
BIND the biomolecular interaction network database died link
DIP the database of interacting proteins http://dip.doe-mbi.ucla.edu/ 
MINT the molecular interaction database http://mint.bio.uniroma2.it/mint/ 
STRING Search Tool for the Retrieval of Interacting Genes/Proteins http://string-db.org/  
HPRO Human protein reference database http://www.hprd.org/ 
BioGRID The Biological General Repository for Interaction Datasets http://thebiogrid.org/ 
这些数据库大部分都还有维护者，还在持续更新，每次更新都可以发一篇paper，而数据库收集的paper引用一般都上千，如果你做了一个数据库，才十几个人引用，那就说明你是自己在跟自己玩。
其中比较好用的是宾夕法尼亚州匹兹堡的大学的一个：http://severus.dbmi.pitt.edu/wiki-pi/
(a) PPI definition; a definition of a protein-to-protein interaction compared to other biomolecular relationships or associations.
(b)PPI determination by two alternative approaches: binary and co-complex; a description of the PPIs determined by the two main types of experimental technologies.
(c) The main databases and repositories that include PPIs; a description and comparison of the main databases and repositories that include PPIs, indicating the type of data that they collect with a special distinction between experimental and predicted data.
(d) Analysis of coverage and ways to improve PPI reliability; a comparative study of the current coverage on PPIs and presentation of some strategies to improve the reliability of PPI data.
(e) Networks derived from PPIs compared to canonical pathways; a practical example that compares the characteristics and information provided by a canonical pathway and the PPI network built for the same proteins. Last, a short summary and guidance for learning more is provided.
现在的蛋白质相互作用数据库的数据都很有限，但是在持续增长，一般有下面四种原因导致数据被收录到数据库
There are four common approaches for PPI data expansions:
1) manual curation from the biomedical literature by experts;
2) automated PPI data extraction from biomedical literature with text mining methods;
3) computational inference based on interacting protein domains or co-regulation relationships, often derived from data in model organisms; and
4) data integration from various experimental or computational sources.
Partly due to the difficulty of evaluating qualities for PPI data, a majority of widely-used PPI databases, including DIP, BIND, MINT, HPRD, and IntAct, take a "conservative approach" to PPI data expansion by adding only manually curated interactions. Therefore, the coverage of the protein interactome developed using this approach is poor.
In the second literature mining approach, computer software replaces database curators to extract protein interaction (or, association) data from large volumes of biomedical literature . Due to the complexity of natural language processing techniques involved, however, this approach often generates large amount of false positive protein "associations" that are not truly biologically significant "interactions".
The challenge for the integrative approach is how to balance quality with coverage.
In particular, different databases may contain many redundant PPI information derived from the same sources, while the overlaps between independently derived PPI data sets are quite low .
参考：
2009年发表的HIPPI数据库：http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-10-S1-S16#CR6_2544 （是对HPRD [11], BIND [20], MINT [21], STRING [26], and OPHID数据库的整合）
一 14

居然可以下载千人基因组计划的所有数据bam，vcf数据

Posted on by
它有两个ftp站点存储所有的数据！
直接看最新版的数据，共有NA编号开头的1182个人，HG开头的1768个人！
每个人的目录下面都有 四个数据文件夹
Oct 01 2014 00:00    Directory alignment
Oct 01 2014 00:00    Directory exome_alignment
Oct 01 2014 00:00    Directory high_coverage_alignment
Oct 01 2014 00:00    Directory sequence_read
这些数据实在是太丰富了！
也可以直接看最新版的vcf文件，记录了这两千多人的所有变异位点信息！
可以直接看到所有的位点，具体到每个人在该位点是否变异！
不过它的基因型信息是通过MVNcall+SHAPEIT这个程序call出来的，具体原理见：http://www.ncbi.nlm.nih.gov/pubmed/23093610
它有两个ftp站点存储所有的数据！
直接看最新版的数据，共有NA编号开头的1182个人，HG开头的1768个人！
每个人的目录下面都有 四个数据文件夹
Oct 01 2014 00:00    Directory alignment
Oct 01 2014 00:00    Directory exome_alignment
Oct 01 2014 00:00    Directory high_coverage_alignment
Oct 01 2014 00:00    Directory sequence_read
这些数据实在是太丰富了！
也可以直接看最新版的vcf文件，记录了这两千多人的所有变异位点信息！
可以直接看到所有的位点，具体到每个人在该位点是否变异！
不过它的基因型信息是通过MVNcall+SHAPEIT这个程序call出来的，具体原理见：http://www.ncbi.nlm.nih.gov/pubmed/23093610
clipboard

 

一 11

批量运行GSEA，命令行版本

Posted on by

之前用过有界面的那种，那样非常方便，只需要做好数据即可，但是如果有非常多的数据，每次都要点击文件，点击下一步，也很烦，不过，，它既然是java软件，就可以以命令行的形式来玩转它！

能够命令行运行了，就很容易批量啦

一、程序安装

直接在官网下载java版本软件即可：http://software.broadinstitute.org/gsea/downloads.jsp

二、输入数据

需要下载gmt文件，自己制作gct和cls文件，或者直接下载测试文件p53

见：http://www.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats

1
三、运行命令

hgu95av2的芯片数据，只有一万多探针，所以很快就可以出结果

 java -cp gsea2-2.2.1.jar  -Xmx1024m xtools.gsea.Gsea   -gmx c2.cp.kegg.v5.0.symbols.gmt  \
 -res P53_hgu95av2.gct  -cls P53.cls   -chip  chip/HG_U95Av2.chip  -out result -rpt_label p53_example
但是一般我们都用默认的即可

2
里面报错说有些探针找不到，不要管它

四、输出数据
3
自己看官网去理解这些结果咯！
需要下载的数据如下：

首先需要下载 Molecular Signatures Database (MSigDB)，一般选择C2的kegg，BioCarta 还有Reactome

都是gmt格式的文件！
CP: Canonical pathways
(browse 1330 gene sets)		 Gene sets from the pathway databases. Usually, these gene sets are canonical representations of a biological process compiled by domain experts. details Download GMT Files
original identifiers
gene symbols
entrez genes ids
CP:BIOCARTA: BioCarta gene sets
(browse 217 gene sets)		 Gene sets derived from the BioCarta pathway database (http://www.biocarta.com/genes/index.asp). Download GMT Files
original identifiers
gene symbols
entrez genes ids
CP:KEGG: KEGG gene sets
(browse 186 gene sets)		 Gene sets derived from the KEGG pathway database (http://www.genome.jp/kegg/pathway.html). Download GMT Files
original identifiers
gene symbols
entrez genes ids
CP:REACTOME: Reactome gene sets
(browse 674 gene sets)		 Gene sets derived from the Reactome pathway database (http://www.reactome.org/). Download GMT Files
original identifiers
gene symbols
entrez genes ids
然后做出表达数据gct文件和cls表型文件~
然后就可以直接运行了
如果是芯片数据，第一列是芯片探针，那么还需要下载chip数据：ftp://ftp.broadinstitute.org/pub/gsea/annotations
一 11

关于芯片平台GPL15308和GPL570

Posted on by

它们虽然被GEO数据标记了不同的ID号，但是其实都是一种芯片，都是昂飞公司的U133++2芯片，分析过芯片数据的人肯定不会陌生了

http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL15308

事实上，这个平台应该是GPL570，但是被CCLE数据库给稍微变通了一下，就给了一个GPL15308的标签，平台主页也写的很清楚，它的探针ID是伪ID，其实就是entrez gene ID

1

本来这个芯片设计的是五万多个探针，最后只剩下了18926个基因
This array is identical to GPL570 but the data were analyzed with a custom CDF Brainarray Version 15, hgu133plus2hsentrezg.
一 11

对CCLE数据库可以做的分析

Posted on by
收集了那么多的癌症细胞系的表达数据，拷贝数变异数据，突变数据，总不能放着让它发霉吧!
这些数据可以利用的地方非常多，但是在谷歌里面搜索引用了它的文章却不多，我挑了其中几个，解读了一下别人是如何利用这个数据的，当然，主要是用那个mRNA的表达数据咯！
这篇文献对CCLE的数据进行了八个步骤的处理，一个合格的生物信息学分析着完全可以重写这个过程
step1:Affymetrix U133 Plus2 DNA microarray gene expressions of 27 gastric cancer cell lines (Kato-III, IM95, SNU-620, SNU-16, OCUM-1, NUGC-4, 2313287, HUG1N, MKN45, NCIN87, KE39, AGS, SNU-5, SNU-216, NUGC-3, NUGC-2, MKN74, MKN7, RERFGC1B, GCIY, KE97, Fu97, SH10TC, MKN1, SNU-1, Hs746 T, HGC27) were downloaded from Cancer Cell Line Encyclopedia (CCLE) [16] in March 2013.
step2: Robust Multi-array Average (RMA) normalization was performed. Principal component analysis plot show no obvious batch effect.
step3: The normalized data is then collapsed by taking the probe sets with highest gene expression.
前三步是为了得到27个胃癌相关细胞系的mRNA表达矩阵，方法是下载cel文件，用RMA归一化，对多探针基因去最大表达量探针！

step4:Unsupervised hierarchical clustering (1-Spearman distance, average linkage) was performed on the cell lines using the aCGH data.

Putative driver genes of which copy number aberrations correlated to mRNA gene expression were identified to determine subtypes or clusters that are driven by different mechanisms. This was done using Mann Whitney U-test with p<0.05, and Spearman Correlation Coefficient test with Rho >0.6.

step5:We then performed consensus clustering[17] on the gene expression data of the 27 gastric cancer cell lines from CCLE using these putative driver genes. We selected k = 2 as it gives sufficiently stable similarity matrix.

step6: In order to assign new samples to this integrative cluster, significance analysis of microarray (SAM) [18]with threshold q<2.0 was used to generate subtype signature based on the mRNA expression data of the 1762 genes from the 27 gastric cancer cell lines in CCLE.

先用甲基化数据来聚类，得到putative driver genes，然后再用这些基因的表达数据来再次聚类，分成两类，然后对这两类进行SAM找差异基因

step7:ssGSEA (single sample GSEA)was used to estimate pathway activities of the gastric cancer cell line in the Molecular Signature Database v3.1 (Msigdb v3.1) [19], [20]. The pathway activities are represented in enrichment scores which were rank normalized to [0.0, 1.0]. 
step8:SAM analysis was performed with threshold q<0.2, and fold change >2.0 (for up-regulated pathways), or <0.5 (for down-regulated pathways) to obtain subtype-specific pathways from the 27 gastric cell lines in CCLE.
这里既用来gene set的富集分析，又用来超几何分布的富集分析，结果去看看这篇文章就知道了！
 
这篇文章只用了CCLE的一个地方，就是看看不同cancer type里面的某个基因表达boxplot
这个图的数据用GEOquery可以得到，样本的分类信息也用GEOquery可以得到，这样就可以做下面这个图了，非常简单
Further, the Cancer Cell Line Encyclopedia (CCLE) database demonstrated that of 1062 cell lines representing 37 distinct cancer types, glioma cell lines express the highest levels of STK17A
1

结论就是：STK17A is highly expressed in glioma cell lines compared to other cancer types. Data was obtained through the Cancer Cell Line Encyclopedia (CCLE).

第三篇文献：http://www.nature.com/ncomms/2013/130709/ncomms3126/fig_tab/ncomms3126_F4.html

这篇文献更简单了，直接对这个表达矩阵进行聚类：
 
The 5,000 most variable genes were used for unsupervised clustering of cell lines by mRNA expression data. Cell lines are colour-coded (vertical bars) according to the reported tissue of origin (a PDF version that can be enlarged at high resolution is in Supplementary Information, Supplementary Fig. S4); horizontal labels at bottom indicate the dominating tissue types within the respective branches of the dendrogram. Most ovarian cancer cell lines (magenta) cluster together, interspersed with endometrial cell lines. However, some ovarian cancer cell lines cluster with other tissue types (*). Top right panels: neighbourhoods (1) of the top cell lines in our analysis, (2) of cell line IGROV1, and (3) of cell line A2780. For the ovarian cancer cell lines in these enlarged areas, the histological subtype as assigned in the original publication is indicated by coloured letters.
就直接拿整个表达矩阵即可，然后挑选变异最大的5000个基因来进行聚类，就可以得到类似的图

 

一 11

CCLE数据库几个知识点

Posted on by
Here we describe the Cancer Cell Line Encyclopedia (CCLE): a compilation of gene expression, chromosomal copy number, and massively parallel sequencing data from 947 human cancer cell lines. 
收集了三种数据：
The mutational status of >1,600 genes was determined by targeted massively parallel sequencing, followed by removal of variants likely to be germline events . 
Moreover, 392 recurrent mutations affecting 33 known cancer genes were assessed by mass spectrometric genotyping13 . 
DNA copy number was measured using high-density single nucleotide polymorphism arrays (Affymetrix SNP 6.0; Supplementary Methods). 
Finally, mRNA expression levels were obtained for each of the lines using Affymetrix U133 plus 2.0 arrays. 
These data were also used to confirm cell line identities .
一般用得最多的就是表达数据，因为表达数据最简单，大多数生物信息学分析着只会用这个数据！
而它的突变数据又不是通常意义的高通量测序得到的，snp6芯片数据很多人听都没听过
文章的附件有对cell lines的具体描述。
different_kinds_of_cancer_in_CCLE
CCLE的数据在broad institute里面可以下载，也放在GEO数据库里面，我比较喜欢GEO里面的数据
This SuperSeries is composed of the following SubSeries:
GSE36133 Expression data from the Cancer Cell Line Encyclopedia (CCLE)
GSE36138 SNP array data from the Cancer Cell Line Encyclopedia (CCLE)
GSE36133这个study的metadata里面有对每个cellline来源的cancer进行描述！
有人喜欢把这个metadata叫做是clinical data。
library(GEOquery)
ccleFromGEO <- getGEO("GSE36133")
annotBlock1 <- pData(phenoData(ccleFromGEO[[1]]))
>dim(annotBlock1)
[1] 917  38
exprSet=exprs(ccleFromGEO[[1]])
> dim(exprSet)
[1] 18926   917
##它的表达数据矩阵，包含了18926个基因，列名是917个细胞系的名字，行是基因的entrez ID
keyColumns <- c("title", "source_name_ch1", "characteristics_ch1", "characteristics_ch1.1", 
    "characteristics_ch1.2")
options(stringsAsFactors = F)
allAnnot=annotBlock1[,keyColumns]
##这几列信息是比较重要的metadata，里面详细记录了细胞系的收集公司单位，tissue，癌症分类等信息
Cell line （1035个细胞系简介）Gene Sets
1035 sets of genes with high or low expression in each cell line relative to other cell lines from the CCLE Cell Line Gene Expression Profiles dataset.
一些关于CCLE数据库的文章：
http://onlinelibrary.wiley.com/doi/10.1002/cncy.21471/pdf 介绍了几个类似的数据库资源
Anticancer drug sensitivity analysis: An integrated approach applied to Erlotinib sensitivity prediction in the CCLE database
一 08

TCGA数据里面的生存分析例子

Posted on by
我们知道了生存分析，就是随着时间的流逝，死亡率是如何增加的，一般是用KM法来估计生存函数，然后画个图即可！而根据某些因子把样本分组，可以看到他们死亡率的变化趋势显著的不同，这就说明了我们的这个因子是非常有效的分类方式，这个因子可以是一个biomarker，也可以某些其它指标！
甚至，我们还可以用cox模型来分析这个因子是如何影响生存函数的，那个稍后再讲
这里，我们就简单讲一个例子，是TCGA里面卵巢癌的数据，根据甲基化数据分成了4个组，那么我们就下载这四个组样本的临床数据，
看看这样分组后，他们的死亡率变化趋势是不是有显著区别！

Continue reading →

一 08

生存分析简介

Posted on by
一般我们谈生存分析，就是说的KM方法估计生存函数，并且画出生存曲线，然后还可以根据分组检验一下它们的生存曲线是否有显著的差异！
在R里面，非常的方便，有个包survival很容易就可以做生存分析了！
只需要记住三个函数即可：
Surv：用于创建生存数据对象
survfit：创建KM生存曲线或是Cox调整生存曲线
survdiff：用于不同组的统计检验

Continue reading →

一 08

别人写的代码运行真快！！！

Posted on by
最近需要做十几万个差异基因分析，每个分析都是对大约5万个探针，200个样本的数据量进行批量T检验计算P值
然后发现自己无论怎么用R来写，每个分析都要耗时半分钟左右，因为我必须循环所有的探针，即使不用for，而用R推荐的apply系列函数，也快不到哪里去，但是我搜索时候发现有一个package里面自带了矩阵T检验，直接对5万个探针进行T检验，而不需要循环处理它们
看下面代码！
dat=matrix(rnorm(10000000),nrow = 50000)
dim(dat) #50000   200
system.time(
  apply(dat,1,function(x){
    t.test(x[1:100],x[101:200])$p.value
  })
)
#用户  系统  流逝
#29.29  0.04 30.64
library(pi0)
system.time(matrix.t.test(dat,1,100,100))
#用户 系统 流逝
#0.48 0.03 0.53
差距真的是非常的明显呀！！！
然后，我解析了它的代码，发现里面调用了C写的代码，我想这就是问题所在咯，可是他们到底怎么写，才能把速度搞这么快？？？
  tmp = .C("tstatistic", dat = x, n1 = n1, n2 = n2, ntests = ntests, 
        MARGIN = MARGIN, pool = pool, tstat = rep(0, ntests), 
        df = rep(0, ntests), PACKAGE = "pi0")

源码在这个package的github里面可以找到，有兴趣的童鞋可以研究一下

 
 

 

一 06

R语言软件的各种旧版本下载

Posted on by
这种编程语言，我之前以为只有包，或者模块什么的才烦人，没想到最近还碰到了版本的问题。
一般来说， 去R的官网，可以下载到基于各种操作系统的最新版R软件，但是在某些特别的情况下，我们可能需要下载以前的旧版本。这时候，想用百度这个破东西搜索出来下载地址简直是痴心妄想。
比如，我想搜索2.15.2这个版本，在Google里面很容易找到了
如果进入bin目录，可以看到各种操作系统的软件https://cran.r-project.org/bin/ ，里面都是可执行的文件，有下面几种操作系统
[DIR] linux/ 23-Jan-2008 19:47 -
[DIR] macos/ 19-Apr-2005 09:45 -
[DIR] macosx/ 12-Dec-2015 09:04 -
[DIR] windows/ 24-Feb-2012 18:41 -
大多数情况下面我们会用window这个平台的，毕竟，可视化嘛！
很多时候，我们也会需要linux平台的https://cran.r-project.org/bin/linux/ 但是linux平台本身就比较复杂
[DIR] debian/ 15-Dec-2015 02:06 -
[DIR] redhat/ 27-Jul-2014 21:12 -
[DIR] suse/ 16-Feb-2012 15:09 -
[DIR] ubuntu/ 06-Jan-2016 04:05 -
我用的Ubuntu比较多，即使选择了ubuntu，里面还有好几个选项，因为Ubuntu也有各种版本！
[DIR] precise/ 06-Jan-2016 04:03 -
[DIR] trusty/ 06-Jan-2016 04:04 -
[DIR] vivid/ 06-Jan-2016 04:04 -
[DIR] wily/ 06-Jan-2016 04:05 -

所以如果，大家是想在linux系统里面安装旧版本的R，建议大家直接下载c源码，直接三部曲就可以安装啦！

[DIR] R-0/ 04-Oct-2004 10:20 -
[DIR] R-1/ 04-Oct-2004 19:02 -
[DIR] R-2/ 01-Mar-2013 09:11 -
[DIR] R-3/ 10-Dec-2015 09:13 -
windows系统里面一般是是exe的可执行文件，直接安装，不需要源码变异，只需要不停的下一步即可
R 3.2.2 (August, 2015)
R 3.2.1 (June, 2015)
R 3.2.0 (April, 2015)
R 3.1.3 (March, 2015)
R 3.1.2 (October, 2014)
R 3.1.1 (July, 2014)
R 3.1.0 (April, 2014)
R 3.0.3 (March, 2014)
R 3.0.2 (September, 2013)
R 3.0.1 (May, 2013)
R 3.0.0 (April, 2013)
R 2.15.3 (March, 2013)
R 2.15.2 (October, 2012)
R 2.15.1 (June, 2012)
R 2.15.0 (March, 2012)
R 2.14.2 (February, 2012)
R 2.14.1 (December, 2011)
R 2.14.0 (November, 2011)
R 2.13.2 (September, 2011)
R 2.13.1 (July, 2011)
R 2.13.0 (April, 2011)
R 2.12.2 (February, 2011)
R 2.12.1 (December, 2010)
R 2.12.0 (October, 2010)
R 2.11.1 (May, 2010)
R 2.11.0 (April, 2010)
R 2.10.1 (December, 2009)
R 2.10.0 (October, 2009)
R 2.9.2 (August, 2009)
R 2.9.1 (June, 2009)
R 2.9.0 (April, 2009)
R 2.8.1 (December, 2008)
R 2.8.0 (October, 2008)
R 2.7.2 (August, 2008)
R 2.7.1 (June, 2008)
R 2.7.0 (April, 2008)
R 2.6.2 (February, 2008)
R 2.6.1 (November, 2007)
R 2.6.0 (October, 2007)
R 2.5.1 (July, 2007)
R 2.5.0 (April, 2007)
R 2.4.1 (December, 2006)
R 2.4.0 (October, 2006)
R 2.3.1 (June, 2006)
R 2.3.0 (April, 2006)
R 2.2.1 (December, 2005)
R 2.2.0 (October, 2005)
R 2.1.1 (June, 2005)
R 2.1.0 (April, 2005)
R 2.0.1 (November, 2004)
R 2.0.0 (October, 2004)
R 1.9.1 (June, 2004)
R 1.8.1 (November, 2003)
R 1.7.1 (June, 2003)
R 1.6.2 (January, 2003)
Installer for R 1.5.1 (June, 2002)
Installer for R 1.4.1 (January, 2002)
Installer for R 1.3.1 (September, 2001)
Binary files for R 1.2.2 (March, 2001)
Binary files for R 1.0.0 (February, 2000)