28

用R的shiny包写一个基因的ID转换小程序

基因的ID现在可能有一千多种,但是我们一般只用NCBI的entrez ID,和HUGO组装规定的gene symbol和gene name
再者有人会用ensembl的ID,至于UCSC的ID很少人用了,更别说剩余的一大堆稀奇古怪的ID了。
ID_example
其实没必要自己写一个软件来转换,因为非常多的现成软件都可以做到!
不过可以拿这个简单的任务来练手!
首先自己做好数据库文件:
我是把数据写到了一个sqlite文件里面了!结构很简单
然后设计软件的界面,俗称UI端,也是很容易,在shiny里面,记住几个控件的函数即可!
UI_1 UI_2
具体代码见我的github代码:https://github.com/jmzeng1314/my-R/tree/master/4-ID-convert-using-shiny
28

在R里面操作SQLite

我前面写到过如何把数据库写到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存储在文件里面!
 2
最后这些数据都保存在了当前工作目录下的hg19_bioconductor.sqlite文件里面!
在其它程序里面就可以直接调用这个文件,而不需要加载一大堆的包了!
1
library(KEGG.db)
library(GO.db)

library(org.Hs.eg.db )

 

28

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

在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"
1
因为我用的merge的all=TRUE,所以最后记录会越来越多
最后的结果就是下面这样,各种ID转换都储存到了my_gene_mapping这个数据对象里面!
可以看到有些基因是没有ensemble数据库对应的ID的,非常多的基因没有被注释到KEGG通路数据库!
我这里如果一个基因对应多个通路,我用冒号连接起来了,成了一个字符串!在后面会有用!

2

现在需要把它上传到数据库里面,这样我就可以用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做基因组坐标转换

其实国际三大主流生物信息学数据库运营单位都出了自己的基因组坐标转换,它们分别是 (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

22

没必要学shell进阶语法

因为大部分生物信息学软件都是linux版本的,所以生物信息学数据分析工作者必备技能就是linux,但是大部分人只是拿他当个中转站,我以前也是,直到接触了大批量的任务,自动化流程,才明白这里面的水太深了,不过无所谓,凭我个人的观点,其实shell的进阶语法真的不必要!

当然,只是我一家之言!
我实在是不想去背诵大括号,小括号,中括号以及双重括号到底区别是什么!
http://www.bio-info-trainee.com/?p=1018  [],[[]],(),(()),{},{{}},以及在前面加上$的区别,以及它们互相杂交组合的区别!!!

我也不想去搞明白操作符两边是否加空格的区别是什么了。

if((i%5==0)) 来判断变量是否被一个数整除
i=$((i+1))来表示变量自增。
这些东西真的很诡异!
如果你有qsub,condor等任务提交系统,那么你只需要熟悉他们就可以了,但大部分散兵游勇的生物信息学家并没有集群,所以压根不会接触任务提交系统,就需要些自动化脚本了!
受限制与机器的cpu以及内存数,需要判断提交了多少任务,等待多久再执行,所以会把一个简单的自动化脚本写的很复杂!
比如下面这个脚本:cat >download_hg38_from_UCSC.sh

for i in $(seq 1 22) X Y M;
do echo $i;
done
gunzip *.gz
for i in $(seq 1 22) X Y M;
do cat chr${i}.fa >> hg38.fa;
done
rm -fr chr*.fa
可以下载hg38基因组的fasta文件,但是是分染色体一个个下载的!
再比如下面这个,批量做GSEA分析的脚本:
while read id
do
echo $id
gene=`echo $id |awk '{print $1}'`
probe=`echo $id |awk '{print $2}'`
echo $i
do_GSEA $probe $gene; ##这里是我自己定义的一个function,就不贴出来了
if((i%5==0))
then
sleep 10  ##重点就在这里,每次提交的任务有限制,所以需要休息,不然机器的cpu负载太高!
fi
i=$((i+1))
done <$1

如果,还有其它功能需要实现,我们可以把脚本写的更负载,纯粹的用shell,需要搜索更多的shell技巧。

但是事实上并没有这个必要,我们现在有了更方便的脚本语言,比如我所擅长的perl
我写一个nohup提交任务的脚本!
## perl nohup.pl   deep_count.sh  0
## perl nohup.pl   deep_count.sh  1
## perl nohup.pl   deep_count.sh  2

[perl]
## perl nohup.pl   deep_count.sh  0
## perl nohup.pl   deep_count.sh  1
## perl nohup.pl   deep_count.sh  2
$i=1;
open FH,$ARGV[0];
while(<FH>){
    chomp;
    next unless $.%3==$ARGV[1];
    $cmd="nohup  $_  &";
    print "$cmd\n";
    system($cmd);
    sleep(10800) if $i%5==4;
    $i++;
    #exit;
}
[/perl]

我尝试过用shell,写了很久,总是报错,但是用perl,一分钟我就写完了,所以,最好是用自己熟悉的一种语法最好!
21

画基因结构图!

一个基因有一个甚至多个转录本!根据转录及翻译现象可以把一个基因人为的定义成多种结构!
首先自己查资料搞明白转录本,外显子和内含子,5端UTR区域和CDS区域,还有3端UTR区域的概念。
真核生物的基因含有外显子和内含子,是区别原核生物的特征之一,能转录的就是外显子,不能转录的就是内含子!
内含子(英语:Intron)是一个基因中非编码DNA片段,它分开相邻的外显子。更精确的定义是:内含子是阻断基因线性表达的序列。DNA上的内含子会被转录到前体RNA中,但RNA上的内含子会在RNA离开细胞核进行转译前被剪除。在成熟mRNA被保留下来的基因部分被称为外显子。
通常我们看到基因结构示意图会像下面这样!
gene-structure
但是对于这个图,我之前有非常多的疑问,直到我自己写程序去仔细统计,画图探究后才真正搞明白!
从图中可以看到5端UTR区域后面接着的是第一外显子,但是事实上不是!
外显子和内含子是转录本上面的概念,是基于转录这个行为的定义。
而5端UTR区域和CDS区域,还有3端UTR区域,是基于翻译这种行为的定义!
把它们画成这样,是严重的干扰信息。
实际上,我们需要更正很多概念,我检验了一下,得到下面几个正确的:
1,如果基因有多个转录本,基因的起始坐标,就是该基因所有转录本的第一个外显子的起始坐标的最小值,同理基因的终止坐标,就是该基因的所有转录本的最后一个外显子的终止坐标的最大值。
2,通过这个概念,可以把基因分成闭合基因和非闭合基因。 闭合基因:有一个最长转录本使得基因起始终止坐标等于该最长转录本的起始终止坐标。(这个是我乱说的,并没有这个定义)
3,如果基因只有一个转录本,那么基因的起始终止坐标,就是转录本的起始终止坐标!
4,一个基因的一个转录本的5'utr区域可以包括多个外显子区域,前者是翻译行为,后者是转录行为
5,起始密码子和终止密码子是CDS的起止处,是基于翻译的概念
6,一个基因的多个转录本的外显子坐标不一定会排列整齐,每个转录本的剪切位点并不一定要比其它转录本一致!
比如对这个ANXA1基因来说,非常多的转录本,但是基因的起始终止坐标,是所有转录本起始终止坐标的极大值和极小值!同时,它是一个闭合基因,因为它存在一个转录本的起始终止坐标等于该基因的起始终止坐标。可以看到它的外显子并不是非常整齐的,虽然多个转录本会共享某些外显子,但是也存在某些外显子比同区域其它外显子长的现象!
 anxa1-gene-structure
再比如下面这个例子:对DDX11L11这个基因来说,前两个外显子都不会翻译,直到第三个外显子才开始翻译,构成CDS。
有些转录本是没有utr的,所以该转录本的起始坐标,就是CDS的起始坐标
5UTR包括多个exon行为
首先把gtf文件格式化导入mysql数据库
我用的是broadinstitute提供的GENCODE的gtf文件,是hg19版本的
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
##format: gtf
##date: 2011-03-23
很简单,下载数据,parse好之后,用R导入mysql即可!
下面是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]

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

 

18

生物信息学分析过程中常见文件格式

刚开始接触生物信息学的时候我也很纠结什么fastq,fastq,sam,bam,vcf,maf,gtf,bed,psl等等,甚至还有过时了的NCBI,ENSEMBL格式,如果是我刚开始 学的时候,我倒是很愿意把他们全部搞透彻,写详细的说明书,但是现在成长了,这些东西感觉很low了,正好我看到了一篇帖子讲数据格式的收集大全,分享给大家,希望初学者能多花点时间好好钻研!

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

每种文件格式的定义,都是有它的道理的,大部分是因为一个比较流行的软件,少量的数据格式是因为国际组织广泛认可而流行的
15

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

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

但是重要而且常用的芯片并不多!
一般分析芯片数据都需要把探针的ID切换成基因的ID,我一般喜欢用基因的entrez ID。
一般有三种方法可以得到芯片探针与gene的对应关系。
金标准当然是去基因芯片的厂商的官网直接去下载啦!!!
一种是直接用bioconductor的包

一种是从NCBI里面下载文件来解析好!
首先,我们说官网,肯定可以找到,不然这种芯片出来就没有意义了!
然后,我们看看NCBI下载的,会比较大
这两种方法都比较麻烦,需要一个个的来!
所以我接下来要讲的是用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)
    #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"))
##自己把内容写出来即可
  }
}

 

 

31

画基因的外显子覆盖度图

一般情况下,我们得到了测序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来看看基因上面突变的分布

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并不是绝对的一一对应的

很多时候,我们都无法确定到底是基于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的思想来做差异基因分析

是福建医科大学的学者开发的,
文章里面详细讲解了他们的这个差异分析的统计学原理
大意就是找到同一组织的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生存分析的风险因子(比例风险模型)

再次强调一下,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

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

其中算法就是,点的坐标该如何确定?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

首先,我们需要了解一下基础知识:
网络图是为了展示数据与数据之间的联系,在生物信息学领域,一般是基因直接的相互作用关系,或者通路之间的联系!
通俗点说,就是我有一些点,它们之间并不是两两相互联系,而是有部分是有连接的,那么我应该如何把这些点画在图片上面呢?因为这些都并没有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文件几次遇到文件头错误

报错如下: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值

这个软件不仅仅能做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

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

华盛顿大学把所有的变异数据都用自己的方法注释了一遍,然后提供下载:
文献是: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)数据库大全

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

它有两个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