二 28

在R里面操作SQLite

Posted on 2016年2月28日 by ulwvfje

我前面写到过如何把数据库写到mysql，但是发现其实msyql并不方便，需要连接数据库什么的，如果发布一个离线小网页，这时候sqlite的优点就显示出来了！

基础代码很简单：

library(RSQLite)

sqlite <- dbDriver("SQLite")

con <- dbConnect(sqlite,"hg19_bioconductor.sqlite") # makes a new file

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

kegg2ID=toTable(org.Hs.egPATH)

#[1] "gene_id" "path_id"

dbWriteTable(con,'keggID2geneID',kegg2ID,row.name=F,overwrite=T)

具体代码，可以看我的github主页：https://github.com/jmzeng1314/my-R/blob/master/3-get-hg19-gene-mapping/get-hg19-gene-mapping-bioconductor.R

做出来的数据，如下，就是几个table存储在文件里面！

最后这些数据都保存在了当前工作目录下的hg19_bioconductor.sqlite文件里面！

在其它程序里面就可以直接调用这个文件，而不需要加载一大堆的包了！

library(KEGG.db)

library(GO.db)

library(org.Hs.eg.db )

一个R语言package：https://cran.r-project.org/web/packages/RSQLite/RSQLite.pdf

参考：https://scottishsnow.wordpress.com/2014/08/14/writing-to-a-database-r-and-sqlite/

很好用，跟mysql一样！

示例代码：http://faculty.washington.edu/kenrice/sisg-adv/exampleSQLite.R

二 28

把bioconductor的gene mapping信息上传到mysql数据库

Posted on 2016年2月28日 by ulwvfje

在R语言里面的bioconductor系列包里面有一个物种注释信息包，其中人类是org.Hs.eg.db

里面有基于人的hg19基因组版本的大部分基因信息之间的转换数据！

包括基因的entrez ID，symbol，name，locus，refseq，kegg pathway，GO ontology对应关系！

但是它们散布在该包的各个数据结构里面！

> ls("package:org.Hs.eg.db")

[1] "org.Hs.eg" "org.Hs.eg.db"

[3] "org.Hs.eg_dbconn" "org.Hs.eg_dbfile"

[5] "org.Hs.eg_dbInfo" "org.Hs.eg_dbschema"

[7] "org.Hs.egACCNUM" "org.Hs.egACCNUM2EG"

[9] "org.Hs.egALIAS2EG" "org.Hs.egCHR"

[11] "org.Hs.egCHRLENGTHS" "org.Hs.egCHRLOC"

[13] "org.Hs.egCHRLOCEND" "org.Hs.egENSEMBL"

[15] "org.Hs.egENSEMBL2EG" "org.Hs.egENSEMBLPROT"

[17] "org.Hs.egENSEMBLPROT2EG" "org.Hs.egENSEMBLTRANS"

[19] "org.Hs.egENSEMBLTRANS2EG" "org.Hs.egENZYME"

[21] "org.Hs.egENZYME2EG" "org.Hs.egGENENAME"

[23] "org.Hs.egGO" "org.Hs.egGO2ALLEGS"

[25] "org.Hs.egGO2EG" "org.Hs.egMAP"

[27] "org.Hs.egMAP2EG" "org.Hs.egMAPCOUNTS"

[29] "org.Hs.egOMIM" "org.Hs.egOMIM2EG"

[31] "org.Hs.egORGANISM" "org.Hs.egPATH"

[33] "org.Hs.egPATH2EG" "org.Hs.egPFAM"

[35] "org.Hs.egPMID" "org.Hs.egPMID2EG"

[37] "org.Hs.egPROSITE" "org.Hs.egREFSEQ"

[39] "org.Hs.egREFSEQ2EG" "org.Hs.egSYMBOL"

[41] "org.Hs.egSYMBOL2EG" "org.Hs.egUCSCKG"

[43] "org.Hs.egUNIGENE" "org.Hs.egUNIGENE2EG"

[45] "org.Hs.egUNIPROT"

其实比较重要的信息，我们需要把它们统一关联起来，做成一个表格，这样可以上传到数据库里面，做成网页，可视化展现，查找！

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

all_EG=mappedkeys(org.Hs.egSYMBOL)

###下面的表格最后是基于entrez ID而关联起来的！

tmp=unlist(as.list(org.Hs.egSYMBOL))

EG2Symbol=data.frame(EGID=names(tmp),symbol=as.character(tmp))

tmp=unlist(as.list(org.Hs.egENSEMBL))

EG2ENSEMBL=data.frame(EGID=names(tmp),ENSEMBL=as.character(tmp))

tmp=unlist(as.list(org.Hs.egGENENAME))

EG2name=data.frame(EGID=names(tmp),name=as.character(tmp))

tmp=unlist(as.list(org.Hs.egMAP))

EG2MAP=data.frame(EGID=names(tmp),MAP=as.character(tmp))

###EG2GO and using mySQL

EG2path=as.list(org.Hs.egPATH)

EG2path=lapply(EG2path, function(x) paste(x,collapse = ":"))

tmp=unlist(EG2path)

EG2path=data.frame(EGID=names(tmp),path=as.character(tmp))

#as.list(head(org.Hs.egGO2ALLEGS))

GO2allEG=as.list(org.Hs.egGO2ALLEGS)

tmp=unlist(GO2allEG)

names(tmp)=substring(names(tmp),1,10) ## change GO:0000002.IMP to GO:0000002

EG2allGO <- tapply(tmp,tmp,function(x){names(x)})

EG2allGO=lapply(EG2allGO, function(x) paste(x,collapse = ":"))

tmp=unlist(EG2allGO)

EG2GO=data.frame(EGID=names(tmp),GO=as.character(tmp))

tmp=merge(EG2Symbol,EG2MAP,by='EGID',all=TRUE)

tmp=merge(tmp,EG2ENSEMBL,by='EGID',all=TRUE)

tmp=merge(tmp,EG2path,by='EGID',all=TRUE)

tmp=merge(tmp,EG2name,by='EGID',all=TRUE)

my_gene_mapping=merge(tmp,EG2GO,by='EGID',all=TRUE)

##[1] "EGID" "symbol" "MAP" "ENSEMBL" "path" "name" "GO"

因为我用的merge的all=TRUE，所以最后记录会越来越多

最后的结果就是下面这样，各种ID转换都储存到了my_gene_mapping这个数据对象里面！

可以看到有些基因是没有ensemble数据库对应的ID的，非常多的基因没有被注释到KEGG通路数据库！

我这里如果一个基因对应多个通路，我用冒号连接起来了，成了一个字符串！在后面会有用！

现在需要把它上传到数据库里面，这样我就可以用R的shiny可视化网页来查找数据库，做ID转换啦！！！

suppressMessages(library(RMySQL))

con <- dbConnect(MySQL(), host="127.0.0.1", port=3306, user="root", password="11111111")

dbSendQuery(con, "USE test")

dbWriteTable(con,'my_gene_mapping',my_gene_mapping)

dbDisconnect(con)

然后就可以在自己的mysql里面看到它啦！！！

其实 "org.Hs.egOMIM" 也很重要，不过我这里只是举个例子，就不深究了！

还有一点，就是那个pathway ID，因为我要以entrez ID为主键，所以把多个pathway给合并了！

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

kegg2ID=toTable(org.Hs.egPATH)

#[1] "gene_id" "path_id"

dbWriteTable(con,'kegg2ID_bioconductor',kegg2ID,row.name=F,overwrite=T)

go2id=toTable(org.Hs.egGO2ALLEGS)

## gene_id go_id Evidence Ontology

dbWriteTable(con,'go2id_bioconductor',go2id,row.name=F,overwrite=T)

dbDisconnect(con)

用这个代码，可以在数据库里面多创建两个表，会有用的！

二 28

用crossmap代替liftover做基因组坐标转换

Posted on 2016年2月28日 by ulwvfje

其实国际三大主流生物信息学数据库运营单位都出了自己的基因组坐标转换，它们分别是 (UCSC liftOver, NCBI Remap, Ensembl API)

Ensembl's Assembly Converter.是基于crossmap的，我觉得挺好用的，就介绍给大家！！！

This online tool currently uses CrossMap, which supports a limited number of formats (see our online documentation for details of the individual data formats listed below). CrossMap also discards metadata in files, so track definitions, etc, will be lost on conversion.

Important note: CrossMap converts WIG files to BedGraph internally for efficiency, and also outputs them in BedGraph format.

但是不知道为什么UCSC的liftover最出名，我也写过它的教程，（http://www.bio-info-trainee.com/?p=990）

Continue reading →

二 21

画基因结构图！

Posted on 2016年2月21日 by ulwvfje

一个基因有一个甚至多个转录本！根据转录及翻译现象可以把一个基因人为的定义成多种结构！

首先自己查资料搞明白转录本，外显子和内含子，5端UTR区域和CDS区域，还有3端UTR区域的概念。

真核生物的基因含有外显子和内含子，是区别原核生物的特征之一，能转录的就是外显子，不能转录的就是内含子！

内含子（英语：Intron）是一个基因中非编码DNA片段，它分开相邻的外显子。更精确的定义是：内含子是阻断基因线性表达的序列。DNA上的内含子会被转录到前体RNA中，但RNA上的内含子会在RNA离开细胞核进行转译前被剪除。在成熟mRNA被保留下来的基因部分被称为外显子。

通常我们看到基因结构示意图会像下面这样！

但是对于这个图，我之前有非常多的疑问，直到我自己写程序去仔细统计，画图探究后才真正搞明白！

从图中可以看到5端UTR区域后面接着的是第一外显子，但是事实上不是！

外显子和内含子是转录本上面的概念，是基于转录这个行为的定义。

而5端UTR区域和CDS区域，还有3端UTR区域,是基于翻译这种行为的定义！

把它们画成这样，是严重的干扰信息。

实际上，我们需要更正很多概念，我检验了一下，得到下面几个正确的：

1，如果基因有多个转录本，基因的起始坐标，就是该基因所有转录本的第一个外显子的起始坐标的最小值，同理基因的终止坐标，就是该基因的所有转录本的最后一个外显子的终止坐标的最大值。

2，通过这个概念，可以把基因分成闭合基因和非闭合基因。闭合基因:有一个最长转录本使得基因起始终止坐标等于该最长转录本的起始终止坐标。（这个是我乱说的，并没有这个定义）

3，如果基因只有一个转录本，那么基因的起始终止坐标，就是转录本的起始终止坐标！

4，一个基因的一个转录本的5'utr区域可以包括多个外显子区域，前者是翻译行为，后者是转录行为

5，起始密码子和终止密码子是CDS的起止处，是基于翻译的概念

6，一个基因的多个转录本的外显子坐标不一定会排列整齐，每个转录本的剪切位点并不一定要比其它转录本一致！

比如对这个ANXA1基因来说，非常多的转录本，但是基因的起始终止坐标，是所有转录本起始终止坐标的极大值和极小值！同时，它是一个闭合基因，因为它存在一个转录本的起始终止坐标等于该基因的起始终止坐标。可以看到它的外显子并不是非常整齐的，虽然多个转录本会共享某些外显子，但是也存在某些外显子比同区域其它外显子长的现象！

再比如下面这个例子：对DDX11L11这个基因来说，前两个外显子都不会翻译，直到第三个外显子才开始翻译，构成CDS。

有些转录本是没有utr的，所以该转录本的起始坐标，就是CDS的起始坐标

首先把gtf文件格式化导入mysql数据库

我用的是broadinstitute提供的GENCODE的gtf文件，是hg19版本的

http://www.broadinstitute.org/cancer/cga/rnaseqc_download

Modified GENCODE GTF file for human with contig names of form ("1","2", etc)

##description: evidence-based annotation of the human genome (GRCh37), version 7 (Ensembl 62)

##provider: GENCODE

##contact: gencode@sanger.ac.uk

##format: gtf

##date: 2011-03-23

很简单，下载数据，parse好之后，用R导入mysql即可！

创建表：http://www.yiibai.com/mysql/mysql_create_tables.html

下面是mysql代码：

--mysql code

use test;

drop table if exists hg19_gtf;

create table hg19_gtf (

gene_name VARCHAR(30),

transcript_name VARCHAR(30) ,

record VARCHAR(15) NOT NULL ,

chr VARCHAR(2) NOT NULL ,

start INT NOT NULL ,

end INT NOT NULL ,

source VARCHAR(10) NOT NULL ,

strand VARCHAR(1) NOT NULL ,

gene_id VARCHAR(30) NOT NULL ,

transcript_id VARCHAR(30) NOT NULL ,

gene_status VARCHAR(30) ,

gene_type VARCHAR(30) ,

transcript_type VARCHAR(30) ,

transcript_status VARCHAR(30)

);

--我的网页好像不支持mysql代码高亮，大家凑合着看吧，反正就是一个简单的建表语句！

select * from hg19_gtf limit 100;

select * from hg19_gtf where gene_name='DMD';

select count(*) from hg19_gtf where gene_name='DMD' and record='start_codon'; --18 start condon

select count(distinct(transcript_name)) from hg19_gtf where gene_name='DMD' ; --34 transcript

select count(distinct(transcript_name)) c ,gene_name from hg19_gtf where record='transcript' group by gene_name order by c desc;

我是随意设计的一个表，主要是为了画图方便！

接下来是R code来具体的对基因进行画图！

[perl]
suppressMessages(library(ggplot2))
suppressMessages(library(RMySQL))
con <- dbConnect(MySQL(), host="127.0.0.1", port=3306, user="root", password="11111111")
dbSendQuery(con, "USE test")
gene='SOX10'
#gene='DDX11L11'
if (T){
query=paste("select * from hg19_gtf where gene_type='protein_coding' and gene_name=",shQuote(gene),sep="")
structure=dbGetQuery(con,query)
tmp_min=min(c(structure$start,structure$end))
structure$new_start=structure$start-tmp_min
structure$new_end=structure$end-tmp_min
tmp_max=max(c(structure$new_start,structure$new_end))
num_transcripts=nrow(structure[structure$record=='transcript',])
tmp_color=rainbow(num_transcripts)
x=1:tmp_max;y=rep(num_transcripts,length(x))
#x=10000:17000;y=rep(num_transcripts,length(x))
plot(x,y,type = 'n',xlab='',ylab = '',ylim = c(0,num_transcripts+1))
title(main = gene,sub = paste("chr",tmp$chr,":",tmp$start,"-",tmp$end,sep=""))
j=0;
tmp_legend=c()
for (i in 1:nrow(structure)){
tmp=structure[i,]
if(tmp$record == 'transcript'){
j=j+1
tmp_legend=c(tmp_legend,paste("chr",tmp$chr,":",tmp$start,"-",tmp$end,sep=""))
}
if(tmp$record == 'exon') lines(c(tmp$new_start,tmp$new_end),c(j,j),col=tmp_color[j],lwd=4)
}
# legend('topleft',legend=tmp_legend,lty=1,lwd = 4,col = tmp_color);

}
[/perl]

通过这个代码，，就能简单的得到我上面显示的基因结构图啦！

二 15

用R获取芯片探针与基因的对应关系三部曲-bioconductor

Posted on 2016年2月15日 by ulwvfje

现有的基因芯片种类不要太多了！

但是重要而且常用的芯片并不多！

一般分析芯片数据都需要把探针的ID切换成基因的ID，我一般喜欢用基因的entrez ID。

一般有三种方法可以得到芯片探针与gene的对应关系。

金标准当然是去基因芯片的厂商的官网直接去下载啦！！！

一种是直接用bioconductor的包

一种是从NCBI里面下载文件来解析好！

首先，我们说官网，肯定可以找到，不然这种芯片出来就没有意义了！

然后，我们看看NCBI下载的，会比较大

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

这两种方法都比较麻烦，需要一个个的来！

所以我接下来要讲的是用R的bioconductor包来批量得到芯片探针与gene的对应关系！

一般重要的芯片在R的bioconductor里面都是有包的，用一个R包可以批量获取有注释信息的芯片平台，我选取了常见的物种，如下：

gpl           organism                  bioc_package

1     GPL32       Mus musculus                        mgu74a

2     GPL33       Mus musculus                        mgu74b

3     GPL34       Mus musculus                        mgu74c

6     GPL74       Homo sapiens                        hcg110

7     GPL75       Mus musculus                     mu11ksuba

8     GPL76       Mus musculus                     mu11ksubb

9     GPL77       Mus musculus                     mu19ksuba

10    GPL78       Mus musculus                     mu19ksubb

11    GPL79       Mus musculus                     mu19ksubc

12    GPL80       Homo sapiens                        hu6800

13    GPL81       Mus musculus                      mgu74av2

14    GPL82       Mus musculus                      mgu74bv2

15    GPL83       Mus musculus                      mgu74cv2

16    GPL85  Rattus norvegicus                        rgu34a

17    GPL86  Rattus norvegicus                        rgu34b

18    GPL87  Rattus norvegicus                        rgu34c

19    GPL88  Rattus norvegicus                         rnu34

20    GPL89  Rattus norvegicus                         rtu34

22    GPL91       Homo sapiens                      hgu95av2

23    GPL92       Homo sapiens                        hgu95b

24    GPL93       Homo sapiens                        hgu95c

25    GPL94       Homo sapiens                        hgu95d

26    GPL95       Homo sapiens                        hgu95e

27    GPL96       Homo sapiens                       hgu133a

28    GPL97       Homo sapiens                       hgu133b

29    GPL98       Homo sapiens                     hu35ksuba

30    GPL99       Homo sapiens                     hu35ksubb

31   GPL100       Homo sapiens                     hu35ksubc

32   GPL101       Homo sapiens                     hu35ksubd

36   GPL201       Homo sapiens                       hgfocus

37   GPL339       Mus musculus                       moe430a

38   GPL340       Mus musculus                     mouse4302

39   GPL341  Rattus norvegicus                       rae230a

40   GPL342  Rattus norvegicus                       rae230b

41   GPL570       Homo sapiens                   hgu133plus2

42   GPL571       Homo sapiens                      hgu133a2

43   GPL886       Homo sapiens                     hgug4111a

44   GPL887       Homo sapiens                     hgug4110b

45  GPL1261       Mus musculus                    mouse430a2

49  GPL1352       Homo sapiens                       u133x3p

50  GPL1355  Rattus norvegicus                       rat2302

51  GPL1708       Homo sapiens                     hgug4112a

54  GPL2891       Homo sapiens                       h20kcod

55  GPL2898  Rattus norvegicus                     adme16cod

60  GPL3921       Homo sapiens                     hthgu133a

63  GPL4191       Homo sapiens                       h10kcod

64  GPL5689       Homo sapiens                     hgug4100a

65  GPL6097       Homo sapiens               illuminaHumanv1

66  GPL6102       Homo sapiens               illuminaHumanv2

67  GPL6244       Homo sapiens   hugene10sttranscriptcluster

68  GPL6947       Homo sapiens               illuminaHumanv3

69  GPL8300       Homo sapiens                      hgu95av2

70  GPL8490       Homo sapiens   IlluminaHumanMethylation27k

71 GPL10558       Homo sapiens               illuminaHumanv4

72 GPL11532       Homo sapiens   hugene11sttranscriptcluster

73 GPL13497       Homo sapiens         HsAgilentDesign026652

74 GPL13534       Homo sapiens  IlluminaHumanMethylation450k

75 GPL13667       Homo sapiens                        hgu219

76 GPL15380       Homo sapiens      GGHumanMethCancerPanelv1

77 GPL15396       Homo sapiens                     hthgu133b

78 GPL17897       Homo sapiens                     hthgu133a

这些包首先需要都下载

gpl_info=read.csv("GPL_info.csv",stringsAsFactors = F)

### first download all of the annotation packages from bioconductor

for (i in 1:nrow(gpl_info)){

  print(i)

  platform=gpl_info[i,4]

  platform=gsub('^ ',"",platform) ##主要是因为我处理包的字符串前面有空格

  #platformDB='hgu95av2.db'

  platformDB=paste(platform,".db",sep="")

  if( platformDB  %in% rownames(installed.packages()) == FALSE) {

    BiocInstaller::biocLite(platformDB)

    #source("http://bioconductor.org/biocLite.R");

    #biocLite(platformDB )

  }

}

下载完了所有的包，就可以进行批量导出芯片探针与gene的对应关系！

for (i in 1:nrow(gpl_info)){

  print(i)

  platform=gpl_info[i,4]

  platform=gsub('^ ',"",platform)

  #platformDB='hgu95av2.db'

  platformDB=paste(platform,".db",sep="")

  if( platformDB  %in% rownames(installed.packages()) != FALSE) {

    library(platformDB,character.only = T)

    #tmp=paste('head(mappedkeys(',platform,'ENTREZID))',sep='')

    #eval(parse(text = tmp))

###重点在这里，把字符串当做命令运行

    all_probe=eval(parse(text = paste('mappedkeys(',platform,'ENTREZID)',sep='')))

    EGID <- as.numeric(lookUp(all_probe, platformDB, "ENTREZID"))

##自己把内容写出来即可

  }

}

参考：http://blog.sina.com.cn/s/blog_62b37bfe0101jbuq.html

一 28

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

Posted on 2016年1月28日 by ulwvfje

how to generate lollipop diagrams ?

问题来源：https://www.biostars.org/p/65347/

这个工具是网页版的，不需要下载，只需上传两列数据即可，就是基因上面的第几个氨基酸是如何突变的这样的信息

网页就可以把这些信息画在基因上面，而且注释到domain和cosmic数据库，挺好用的！

http://www.cbioportal.org/mutation_mapper.jsp

Hugo_Symbol	HUGO symbol for the gene	TP53
Protein_Change	Amino acid change	V600E

结果也很容易理解！

我输入了4个基因，网页就会对每个基因画一个图，其实ABCA1这个基因有6个突变，就都画在该基因上面了，基因突变分布在基因的哪个domain也看得一清二楚！

我的输入文件是这样的！

非常好用：

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)

甚至还有高手用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.

Live examples are found at BioJS: http://registry.biojs.net/client/#/detail/muts-needle-plot
Code is available at GitHub: https://github.com/bbglab/muts-needle-plot
And it is also a npm-package: https://www.npmjs.com/package/muts-needle-plot

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.

http://www.ensembl.org/Homo_sapiens/Transcript/ProteinSummary?db=core;g=ENSG00000141510;r=17:7565097-7590856;t=ENST00000269305

一 22

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

Posted on 2016年1月22日 by ulwvfje

是福建医科大学的学者开发的，

文章里面详细讲解了他们的这个差异分析的统计学原理

大意就是找到同一组织的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精确检验的结果是

这个基因在该疾病样本，相对于normal pool并没有差异表达！当然检验得到的P值最后可以做FDR校验。

依次这样，把所有的gene都分析完，就知道这个样本有哪些差异的gene了。

介绍完原理，我们拿一个具体的例子来看看吧：

首先我们下载一个2008年的一个人的肝脏表达数据（Gene Expression in Human Liver），都是正常组织，共427个样本。

不过这个芯片比较小众，是默克医药公司定制化的，需要下载探针对应gene的文件！

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

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

我们读取GSE9588这个数据，得到表达矩阵，然后计算rank矩阵，然后计算得到comp矩阵

> table(rank_comp)

rank_comp

down no up

58479 465752098 58479

不知道为什么这个数据，stable的那么少，不知道是不是出了什么问题！

其实我的程序都是对的了，只是因为这个数据集已经不是纯粹的表达量数据了，而是这427个样本的数据都减去了某个样本的表达量。

这样每个个体的基因之间的表达量排序就会被干扰，导致得到的稳定基因对非常少！！！

但是，我后来下载了GTEx的表达数据，拿那里面的normal组织样本表达量来做，可以得到非常多的稳定基因对。

实际代码大概是：

得到正常组织的表达矩阵：

然后计算表达矩阵的rank，得到各个样本自己的基因排序情况，得到排序矩阵！

处理排序矩阵，每个基因对之间都算一下是否稳定，得到稳定性描述矩阵！

然后根据每个疾病个体的基因表达情况，来循环每个基因，看看该基因是否差异！

一 14

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

Posted on 2016年1月14日 by ulwvfje

报错如下：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能系统性的学习一些基础知识！

参考：http://seqanswers.com/forums/showthread.php?t=17233

https://www.biostars.org/p/115819/

http://gatkforums.broadinstitute.org/gatk/discussion/2667/bam-is-malformed-depthofcoverage

一 14

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

Posted on 2016年1月14日 by ulwvfje

这个软件不仅仅能做QC，而且可以统计各个基因的RPKM值！尤其是TCGA计划里面的都是用它算的

一、程序安装

直接在官网下载java版本软件即可使用：http://www.broadinstitute.org/cancer/cga/tools/rnaseqc/RNA-SeQC_v1.1.8.jar

但是需要下载很多注释数据

二、输入数据

箭头所指的文件，一个都不少，只有那个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的信息

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.

软件的主页是：

http://www.broadinstitute.org/cancer/cga/rnaseqc_run

http://www.broadinstitute.org/cancer/cga/rnaseqc_download

帮助文件：http://www.broadinstitute.org/cancer/cga/sites/default/files/data/tools/rnaseqc/RNA-SeQC_Help_v1.1.2.pdf

一 14

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

Posted on 2016年1月14日 by ulwvfje

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

文献是：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

这些数据在很多时候非常有用，尤其是想跟自己得到的突变数据做交叉验证，或者做一下统计分析的时候！

人的基因组才300亿个位点，他们就注释了86亿！！！

所以有三百多G的压缩包数据，我想，一般的公司或者单位都不会去用这个数据了！

一 14

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

Posted on 2016年1月14日 by ulwvfje

最近遇到一个项目需要探究一个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:

Organisms: Homo sapiens (Human)
Availability: Free to all users
Standards: all

人类的六个主要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://openwetware.org/wiki/Protein-protein_interaction_databases

其中比较好用的是宾夕法尼亚州匹兹堡的大学的一个：http://severus.dbmi.pitt.edu/wiki-pi/

http://severus.dbmi.pitt.edu/wiki-pi/index.php/pair/view/3838/7157

(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数据库的整合）

2010年的综述：http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1000807

http://bib.oxfordjournals.org/content/early/2010/09/16/bib.bbq064.full

一 14

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

Posted on 2016年1月14日 by ulwvfje

它有两个ftp站点存储所有的数据！

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/

ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/

直接看最新版的数据，共有NA编号开头的1182个人，HG开头的1768个人！

ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/phase3/data/

也可以按照人种来查看这些数据：ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000_genomes_project/data/

每个人的目录下面都有四个数据文件夹

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文件，记录了这两千多人的所有变异位点信息！

可以直接看到所有的位点，具体到每个人在该位点是否变异！

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/

不过它的基因型信息是通过MVNcall+SHAPEIT这个程序call出来的，具体原理见：http://www.ncbi.nlm.nih.gov/pubmed/23093610

而且网站还提供一些教程：ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000_genomes_project/working/

它有两个ftp站点存储所有的数据！

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/

ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/

直接看最新版的数据，共有NA编号开头的1182个人，HG开头的1768个人！

ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/phase3/data/

也可以按照人种来查看这些数据：ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000_genomes_project/data/

每个人的目录下面都有四个数据文件夹

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文件，记录了这两千多人的所有变异位点信息！

可以直接看到所有的位点，具体到每个人在该位点是否变异！

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/

不过它的基因型信息是通过MVNcall+SHAPEIT这个程序call出来的，具体原理见：http://www.ncbi.nlm.nih.gov/pubmed/23093610

而且网站还提供一些教程：ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000_genomes_project/working/

还有Illumina的450K甲基化芯片数据：ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/hgsv_sv_discovery/working/20151214_450k_methylation/

还有一个小程序，ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/technical/browser/vcf_to_ped_converter/version_1.1/

还有Illumina的450K甲基化芯片数据：ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/hgsv_sv_discovery/working/20151214_450k_methylation/

还有一个小程序，ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/technical/browser/vcf_to_ped_converter/version_1.1/

一 11

批量运行GSEA，命令行版本

Posted on 2016年1月11日 by ulwvfje

之前用过有界面的那种，那样非常方便，只需要做好数据即可，但是如果有非常多的数据，每次都要点击文件，点击下一步，也很烦，不过，，它既然是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

三、运行命令

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

这个参数其实非常多，见文件http://software.broadinstitute.org/gsea/doc/linkedFiles/GSEAParameters.txt

但是一般我们都用默认的即可

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

四、输出数据

自己看官网去理解这些结果咯！

需要下载的数据如下：

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

http://software.broadinstitute.org/gsea/downloads.jsp

都是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表型文件~

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

然后就可以直接运行了

如果是芯片数据，第一列是芯片探针，那么还需要下载chip数据：ftp://ftp.broadinstitute.org/pub/gsea/annotations

一 06

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

Posted on 2016年1月6日 by ulwvfje

这种编程语言，我之前以为只有包，或者模块什么的才烦人，没想到最近还碰到了版本的问题。

一般来说，去R的官网，可以下载到基于各种操作系统的最新版R软件，但是在某些特别的情况下，我们可能需要下载以前的旧版本。这时候，想用百度这个破东西搜索出来下载地址简直是痴心妄想。

比如，我想搜索2.15.2这个版本，在Google里面很容易找到了

https://cran.r-project.org/bin/windows/base/old/2.15.2/

https://cran.r-project.org/bin/windows/base/old/2.15.2/R-2.15.2-win.exe

然后我看了一下域名，发现其实很有规律的

如果进入bin目录，可以看到各种操作系统的软件https://cran.r-project.org/bin/ ，里面都是可执行的文件，有下面几种操作系统

linux/	23-Jan-2008 19:47	-
macos/	19-Apr-2005 09:45	-
macosx/	12-Dec-2015 09:04	-
windows/	24-Feb-2012 18:41	-

大多数情况下面我们会用window这个平台的，毕竟，可视化嘛！

很多时候，我们也会需要linux平台的https://cran.r-project.org/bin/linux/ 但是linux平台本身就比较复杂

debian/	15-Dec-2015 02:06	-
redhat/	27-Jul-2014 21:12	-
suse/	16-Feb-2012 15:09	-
ubuntu/	06-Jan-2016 04:05	-

我用的Ubuntu比较多，即使选择了ubuntu，里面还有好几个选项，因为Ubuntu也有各种版本！


precise/	06-Jan-2016 04:03	-
trusty/	06-Jan-2016 04:04	-
vivid/	06-Jan-2016 04:04	-
wily/	06-Jan-2016 04:05	-

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

源码是在：https://cran.r-project.org/src/base/

R-0/	04-Oct-2004 10:20	-
R-1/	04-Oct-2004 19:02	-
R-2/	01-Mar-2013 09:11	-
R-3/	10-Dec-2015 09:13	-

windows系统里面一般是是exe的可执行文件，直接安装，不需要源码变异，只需要不停的下一步即可

所有的版本见：https://cran.r-project.org/bin/windows/base/old/

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)

一 06

picnic对拷贝数变异检测芯片数据进行分析

Posted on 2016年1月6日 by ulwvfje

这里的拷贝数变异检测芯片指的是Affymetrix Genome-Wide Human SNP Array 6.0

cel数据，处理成segment及genotype数据

一、程序安装

这本来是一个matlab程序，但是有linux版本，需要安装matlab编译环境

下载地址是 ftp://ftp.sanger.ac.uk/pub/cancer/picnic_software/picnic_src/

下载解压之后首先安装matlab环境：

./MCRInstaller.bin -console

因为我的服务器没有界面，所以我用了-console执行安装程序，很简单就安装好了

其中，cdf目录下面是我们从http://www.affymetrix.com 里面下载的Library FilesGenome-Wide Human SNP Array 6.0 (zip, 246 MB)

celConverter目录下面是一个java程序，可以把芯片出来的cel文件转为flat file

其中Matlab_running 是我执行安装的matlab环境

其余两个脚本run_HMM.sh run_preprocessing.sh是picnic程序的第二步和第三步

二、输入数据准备

我随便找了两个snp6.0芯片的raw data

-rw-rw-r-- 1 jmzeng jmzeng 66M Dec 30 06:30 GSM1949207.CEL

-rw-rw-r-- 1 jmzeng jmzeng 66M Sep 9 11:08 GSM887898.CEL

三、程序使用

程序主要分为三个步骤来实现

Step 1: Convert the binary *.cel file to a flat file
Step 2: Normalise and estimate the ploidy of the genome and the level of normal contamination
Step 3: Segment the data and produce the genotype

第一步：

暂时还不需要matlab工作环境

示例的code是：java -Xmx2G -jar CelFileConverter.jar -m Snp6FeatureMappings.csv -c 'cdf_file_including_path' -s 'directory name
of cel files' -t rootDir/outdir/raw

我的code是：

jmzeng@ubuntu:/home/jmzeng/bio-soft/picnic/c_code/celConverter$ java -jar CelFileConverter.jar -m Snp6FeatureMappings.csv -c ../cdf/GenomeWideSNP_6.Full.cdf -s ../celFiles/ -t ../result/
Loading CDF file: ../cdf/GenomeWideSNP_6.Full.cdf
Loading feature mapping file: Snp6FeatureMappings.csv
Processing CEL file: /home/jmzeng/bio-soft/picnic/c_code/celConverter/../celFiles/GSM1949207.CEL
Created feature intensity file: /home/jmzeng/bio-soft/picnic/c_code/celConverter/../result/GSM1949207.feature_intensity
Processing CEL file: /home/jmzeng/bio-soft/picnic/c_code/celConverter/../celFiles/GSM887898.CEL
Created feature intensity file: /home/jmzeng/bio-soft/picnic/c_code/celConverter/../result/GSM887898.feature_intensity

然后在输出目录就有了一个feature_intensity后缀的文本文件，约110M

-rw-rw-r-- 1 jmzeng jmzeng 105M Dec 30 06:35 GSM1949207.feature_intensity

-rw-rw-r-- 1 jmzeng jmzeng 109M Dec 30 06:35 GSM887898.feature_intensity

第二步：

这里需要用matlab啦，但是这里经常出现库的问题！！！

示例的code是

sh run_preprocessing.sh mcr_dir cell_name info_dir feature_int_dir normalised_outdir outdir sample_type in_pi
我的code是：

jmzeng@ubuntu:/home/jmzeng/bio-soft/picnic/c_code$ sh run_preprocessing.sh Matlab_running/v710/ GSM1949207 info/ result/ result/output result/
------------------------------------------
Setting up environment variables
---
LD_LIBRARY_PATH is .:Matlab_running/v710//runtime/glnxa64:Matlab_running/v710//bin/glnxa64:Matlab_running/v710//sys/os/glnxa64:Matlab_running/v710//sys/java/jre/glnxa64/jre/lib/amd64/native_threads:Matlab_running/v710//sys/java/jre/glnxa64/jre/lib/amd64/server:Matlab_running/v710//sys/java/jre/glnxa64/jre/lib/amd64/client:Matlab_running/v710//sys/java/jre/glnxa64/jre/lib/amd64
My Own Exception: Fatal error loading library /home/jmzeng/bio-soft/picnic/c_code/Matlab_running/v710/bin/glnxa64/libmwmclmcr.so Error: libXp.so.6: cannot open shared object file: No such file or directory

第三步：

也需要用matlab，库的问题必须解决！！！

我另外一台服务器的结局办法是用v714的matlab

示例的code是

sh run_HMM.sh /nfs/team78pc2/kwl_temp/segments/PICNIC/C/release/Matlab_Compiler_Runtime/v710 A01_CGP_PD3945a.feature_intensity '/nfs/team78pc3/KWL/segments/PICNIC/matlab/C/release/info/' '/nfs/team78pc2/kwl_temp/segments/PICNIC/data/normalized/' '/nfs/team78pc2/kwl_temp/segments/PICNIC/data/' '10' '0.33598' '1.9915' '0.40997'
我的code是：无所谓了，这个服务器不知道怎么回事，总是出现库文件的问题，而这个问题需要root权限，我懒得弄了，我在其它的服务器上面都木有问题的！！！

所以我就换了MATLAB，毕竟，这个软件本来就是matlab版本的，第一步照旧，第二步在matlab里面运行！

我这里只用GSM1949207做测试吧！！！

Genomic DNA was extracted from saliva, peripheral blood, or fibroblast cell lines using the QIAamp DNA Blood Mini Kit or QIAamp DNA Mini Kit. DNA quality and quantity was assessed using a Nanodrop Spectrophotometer and agarose gel electrophoresis.

打开matlab，进入picnic目录

重复第二步，输入：

preprocessing('GSM1949207.feature_intensity','info\','result\raw\','result\output\','result\')

需要7个参数，我这个是cell lines数据，所以后面两个参数省略不写！info文件夹自己下载放在picnic目录，其中result\raw 存放你第一步的结果文件！

这一步运行，会比较久！

同时可以看到程序在我的result目录里面新增加了两个目录用来存放结果，不过这一步的结果还是中间文件，就不解释了！

然后再重复第三步，输入：

HMM('GSM1949207.feature_intensity','info\','result\output\','result\',10,0,2.0221,0.40997)

这个软件的参数很没规律，最好不要像说明书那些用output做文件夹名字，不然，很容易出错。

这一步好像也很耗时间！

重要的结果有两个：

一 06

拷贝数变异检测芯片介绍

Posted on 2016年1月6日 by ulwvfje

这里的拷贝数变异检测芯片指的是Affymetrix Genome-Wide Human SNP Array 6.0

cel数据，需要处理成segment及genotype数据

这个芯片在TCGA计划里面用的非常多，是标配了。大家只要记住，这是一个跟拷贝数变异检测相关的芯片，而且还可以测一些genotype

Affymetrix Genome-Wide Human SNP Array 6.0是唯一可以真正将CNP(拷贝数多态性)转化成高分辨率的参考图谱的平台。主要应用领域包括全基因组SNP分型、全基因组CNV分型、全基因组关联分析、全基因组连锁分析。除了进行基因分型外，还为拷贝数研究和LOH研究提供帮助，从而能够进行：UPD检测、亲子鉴定、异常的亲代起源分析（针对 UPD和缺失）、纯合性分析、血缘关系鉴定。

参考：http://www.affymetrix.com/support/technical/byproduct.affx?product=genomewidesnp_6

SNP Array 6.0是昂飞公司继Mapping10k、100k、500k和SNP5.0芯片后推出的新一代SNP芯片。在一张芯片上可以分析一个样本906,600 个SNP的基因型, 大约有482，000个SNP来自于前代产品500K和SNP5.0芯片。剩下424，000个SNP包括了来源于国际HapMap计划中的标签 SNP，X，Y染色体和线粒体上更具代表性的SNP,以及来自于重组热点区域和500K芯片设计完成后新加入dbSNP数据库的SNP。该芯片同时含 946,000个非多态性CNV探针，用于检测拷贝数变异，其中202,000个用于检测5677个已知拷贝数变异区域的探针，这些区域来源于多伦多基因组变异体数据库。该数据库中每隔3,182个非重叠片段区域分别用61个探针来检测。除了检测这些已知的拷贝数多态区域，还有超过744,000个探针平均分配到整个基因组上，用来发现未知的拷贝数变异区域。SNP和CNV两种探针高密度且均匀地分布在整个基因组，作为拷贝数变异和杂合性缺失(LOH)检测的工具来发现微小的染色体增加和缺失。为广大生命科学研究者提高发现复杂疾病相关基因的可能提供了强有力的工具。
通过与哈佛大学合办的Broad研究所合作，SNP6.0芯片在数据准确性和一致性方面达到了新的高度。相应推出的Genotyping Console用来处理SNP6.0芯片数据和全基因组遗传分析及质量控制。

产品特点：

1.涵盖超过1,800,000个遗传变异标志物：包括超过906,600个SNP和超过946,000个用于检测拷贝数变化（CNV，Copy Number Variation）的探针；

2.SNP和CNV两种探针高密度且均匀地分布在整个基因组，不仅可以用于SNP基因精确分型，还可用于拷贝数变异CNV的研究；

3.744,000个探针平均分配到整个基因组上，用来发现未知的拷贝数变异区域；

4.可用于Copy-neutral LOH/UPD检测，亲子鉴定，纯合性分析、血缘关系鉴定、遗传病或其它疾病的研究。

参考：http://www.biomart.cn/specials/cnv2014/article/84169

在NCBI的GEO数据库里面可以查到这个芯片，已经有一万多个样本数据啦!

图中第一个是CCLE计划的近千个样本，可能是定制化了的snp6.0芯片吧

使用这个芯片数据来发文章的非常多，见列表：http://media.affymetrix.com/support/technical/other/snp6_array_publications.pdf

还有一篇2010-nature文章讲了如何用picnic来研究cnv，http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3145113/

也有一篇2010年的文章提出了新的软件来分析这个芯片cnv数据http://bioinformatics.oxfordjournals.org/content/26/11/1395.long

实现同样功能的软件，非常之多，还有一个R的bioconductor系列的包

http://www.bioconductor.org/help/search/index.html?q=cnv/

随便进去都可以找到很多raw data，可以自己进行分析的！

http://www.ncbi.nlm.nih.gov/geo/browse/?view=samples&platform=6801

比如：ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1949nnn/GSM1949207/suppl/GSM1949207%5FSB%5FCID0102B%5F071708%2ECEL%2Egz

十二 29

用firehose_get 来下载所有TCGA寄存在broad的数据

Posted on 2015年12月29日 by ulwvfje

该软件是broad institute写的一个数据接口，主要是供他人下载TCGA的所有寄存在broad institute的免费数据，主要是level3，level4的数据。(说错了，好像只有level4的数据，就是可以发文章的分析结果及图片)

软件下载地址：https://confluence.broadinstitute.org/display/GDAC/Download

懂它的使用规则，编码规则即可：

https://tcga-data.nci.nih.gov/datareports/codeTablesReport.htm

就是一个很简单的shell脚本而已，根据几个用户自定义参数来选择性的下载数据。

我们可以用-t这个参数来指定下载的数据类型，可以是mut/rna/mutsig/gistic等各种数据，至于这些单词代表什么意义，需要自己去看说明书啦

还可以指定时间，截止到什么时间的数据！

它支持的癌症种类：

ACC  BLCA  BRCA  CESC  COAD  COADREAD  DLBC  ESCA  
	GBM  HNSC  KICH  KIRC  KIRP  LAML  LGG  LIHC  
	LUAD  LUSC  OV  PAAD  PANCANCER  PANCAN8  PANCAN12  PRAD  
	READ  SARC  SKCM  STAD  THCA  UCEC  UCS

这些癌症种类的简称，也是可以去官网里面看到的！官网：http://gdac.broadinstitute.org

十二 29

所有TCGA收集的mRNA表达数据集数据集-GSE62944

Posted on 2015年12月29日 by ulwvfje

无意中看一篇文献，有提到这个数据集，我简单看了一下，居然还真的这么全面！！！

它处理了目前（大概是2015年6月）TCGA收集的所有癌症样本的mRNA表达数据，并且统一处理成了count和RPKM两种表达量形式。

GEO地址：http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62944

Title	Alternatively processed and compiled RNA-Sequencing and clinical data for thousands of samples from The Cancer Genome Atlas
Organism	Homo sapiens
Experiment type	Expression profiling by high throughput sequencing
Summary	We reprocessed RNA-Seq data for 9264 tumor samples and 741 normal samples across 24 cancer types from The Cancer Genome Atlas with "Rsubread". Rsubread is an open source R package that has shown high concordance with other existing methods of alignment and summarization, but is simple to use and takes significantly less time to process data. Additionally, we provide clinical variables publicly available as of May 20, 2015 for the tumor samples where the TCGA ids are matched.

这样就非常方便大家使用了。

可以直接拿来进行聚类，看看不同癌症种类如何聚集在一起，也可以直接拿来跟某个临床指标进行关联分析.

比如，如果我们想用DEseq来做差异分析，那么需要的数据就是基因的count，这里就有近一万个癌症样本的基因表达的count值，随便都可以分类看看差异情况，再富集分析分析。非常适合初学者练手。

癌症种类见官网：http://gdac.broadinstitute.org/

十二 29

自动化出网页报告的R语言包-Nozzle

Posted on 2015年12月29日 by ulwvfje

挺好用的：https://confluence.broadinstitute.org/display/GDAC/Nozzle

http://bioinformatics.oxfordjournals.org/content/29/8/1089.full.pdf+html

自动化出生物信息学分析报告：

一个简单的例子：https://raw.githubusercontent.com/parklab/Nozzle/develop/examples/basic.R

挺好玩的

开发者自己就用这个包写了一大堆的自动化报告；

http://gdac.broadinstitute.org/runs/analyses__latest/reports/cancer/OV/Mutation_Assessor/nozzle.html

很牛！

还有很多其它例子：https://github.com/parklab/Nozzle/tree/develop/examples

功能还是很多的，https://cran.r-project.org/web/packages/Nozzle.R1/Nozzle.R1.pdf

根据测试代码输出的报告如下：

这些报告里面的元素，都是测试代码添加的，很容易就能理解

十二 29

用Mutation-Assessor软件来看突变位点对基因或者蛋白功能的影响

Posted on 2015年12月29日 by ulwvfje

这是一个在线工具，非常好用，对snp位点进行注释，看看该突变是否影响蛋白功能，一定要收藏！！！

官网：http://mutationassessor.org/

也应该是有standalone版本，我没有去找，不过网页版就很好用了，只需要复制粘贴进去自己想分析的数据，按照一定的格式即可，比如：

该软件就能分析出该突变位点发生在BRCA2这个基因上面，对氨基酸的改变也能写出来，对蛋白功能的改变等选项都是可以自由定制化的。

输入数据非常广泛：

The server accepts list of variants, one variant per line, plus optional text describing your variants,
in genomic coordinates, "+" strand assumed :
<genome build>,<chromosome>,<position>,<reference allele>,<substituted allele>
Genome build is optional (build 18 assumed), accepted values: 'hg18' and 'hg19'
Examples:

hg19,13,32912555,G,T BRCA2
hg18,7,55178574,G,A GBM
7,55178574,G,A GBM

or in protein space: <protein ID> <variant> <text>, where <protein ID> can be :

1. Uniprot protein accession (i.e. EGFR_HUMAN)
2. NCBI Refseq protein ID (i.e. NP_005219)

examples:

EGFR_HUMAN R521K
EGFR_HUMAN R98Q Polymorphism
EGFR_HUMAN G719D disease
NP_000537 G356A
NP_000537 G360A dbSNP:rs35993958
NP_000537 S46A Abolishes phosphorylation

ID types can be mixed in one list in any way.

Page 6 of 10« First ‹ Previous 3 4 567 8 9 Next ›Last »

生信菜鸟团

欢迎去论坛biotrainee.com留言参与讨论，或者关注同名微信公众号biotrainee

Category Archives: 生信基础

在R里面操作SQLite

把bioconductor的gene mapping信息上传到mysql数据库

用crossmap代替liftover做基因组坐标转换

画基因结构图！

用R获取芯片探针与基因的对应关系三部曲-bioconductor

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

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

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

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

这个软件不仅仅能做QC，而且可以统计各个基因的RPKM值！尤其是TCGA计划里面的都是用它算的

一、程序安装

TCGA数据里面都会提供由RNA-SeQC软件处理得到的表达矩阵！

Expression

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

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

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

批量运行GSEA，命令行版本

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

picnic对拷贝数变异检测芯片数据进行分析

拷贝数变异检测芯片介绍

用firehose_get 来下载所有TCGA寄存在broad的数据

所有TCGA收集的mRNA表达数据集数据集-GSE62944

自动化出网页报告的R语言包-Nozzle

用Mutation-Assessor软件来看突变位点对基因或者蛋白功能的影响

2024年11月
一	二	三	四	五	六	日
« 九
				1	2	3
4	5	6	7	8	9	10
11	12	13	14	15	16	17
18	19	20	21	22	23	24
25	26	27	28	29	30