前面我讲到TCGA的数据可以在5个组织机构可以获取，他们都提供了类似的接口来供用户下载数据
每个接口都有使用教程，比如http://firebrowse.org/tutorial/FireBrowse-Tutorial.pdf
非常详细！！！
有的还专门写了软件接口：https://confluence.broadinstitute.org/display/GDAC/Download
或者写了R的接口：http://www.cbioportal.org/cgds_r.jsp
接下来我们要讲的就是cbioportal网站提供的一个R接口，非常好用，只需记住4个函数即可！！！ Continue reading
The molecular taxonomy of primary prostate cancer
Cell Volume 163 Issue 4: p1011-1025 Read the full article
Portal Publication Site and Associated Data Files
Comprehensive Molecular Characterization of Papillary Renal Cell Carcinoma
NEJM. Published on line on Nov 4th, 2015 Read the full article
Portal Publication Site and Associated Data Files
[perl]
open FH,"/home/jmzeng/hg19/chrY.gene.special.position" or die "file error !!!";
while(<FH>){
chomp;
@F=split;
foreach ($F[1]..$F[2]){
$h{$_}=$F[3];
}
$length{$F[3]}=$F[2]-$F[1]+1;
}
close FH;
open FH,$ARGV[0];
while(<FH>){
chomp;
@F=split;
next unless $F[0] eq 'chrY';
next if $F[2]<20;
if (exists $h{$F[1]}){
$count{$h{$F[1]}}++ ;
}else{
$count{'other'}++ ;
}
}
close FH;
print "$_\t$length{$_}\t$count{$_}\n" foreach sort keys %count;</pre>
</div>
<div>[/perl]
| AMELY | 8111 | 1269 |
| BCORP1 | 47724 | 689 |
| CSPG4P1Y | 3799 | 538 |
| DAZ1 | 69739 | 762 |
| DAZ2 | 71901 | 228 |
| DAZ3 | 73222 | 233 |
| DAZ4 | 73222 | 540 |
| DDX3Y | 12825 | 3654 |
| EIF1AY | 17445 | 929 |
| FAM224A | 4295 | 82 |
| FAM224B | 4293 | 85 |
| GOLGA2P3Y | 4866 | 68 |
| GYG2P1 | 15476 | 547 |
| HSFY2 | 42277 | 3950 |
| KDM5D | 39526 | 7425 |
| NLGN4Y | 319396 | 3872 |
| PCDH11Y | 105374 | 6627 |
| PRKY | 107577 | 1390 |
| PRORY | 3388 | 735 |
| RBMY1B | 14451 | 232 |
| RBMY1D | 14411 | 117 |
| RBMY1E | 14410 | 157 |
| RBMY1J | 14407 | 65 |
| RBMY2EP | 6416 | 27 |
| RBMY2FP | 7348 | 419 |
| RPS4Y1 | 25376 | 1856 |
| RPS4Y2 | 24966 | 1831 |
| SRY | 888 | 703 |
| TBL1Y | 180999 | 3231 |
| TGIF2LY | 958 | 808 |
| TMSB4Y | 2457 | 534 |
| TSPY4 | 132211 | 1525 |
| TTTY14 | 205048 | 394 |
| TTTY4C | 36811 | 39 |
| TTTY9A | 9317 | 580 |
| TXLNGY | 23067 | 1968 |
| USP9Y | 159610 | 10508 |
| UTY | 232293 | 6670 |
| VCY | 742 | 291 |
| XKRY2 | 1582 | 980 |
| ZFY | 47437 | 3125 |
| other | 100328 |
| NLGN4Y | 319396 | 575 |
| PCDH11Y | 105374 | 1643 |
| PRKY | 107577 | 82 |
| TGIF2LY | 958 | 191 |
| TTTY14 | 205048 | 139 |
| other | 54297 |
cat Sample5.gatk.UG.vcf |perl -alne '{next unless $F[0] eq "chrX";next unless $F[5]>30;$h{(split/:/,$F[9])[0]}++}END{print "$_\t$h{$_}" foreach keys %h}'
之前我们用超几何分布检验的方法做了富集分析，使用的是GSE63067.diffexp.NASH-normal.txt的logFC的绝对值大于0.5，并且P-value小雨0.05的基因作为差异基因来检验kegg的pathway的富集情况
结果是这样的
我们接下来用另外一种方法来做富集分析，顺便检验一下，是不是超几何分布统计检验的富集分析方法就是最好的呢？
这种方法是-重抽样+主成分分析
大概的原理是，比如对上图中的，04380这条pathway来说，总共有128个基因，那么我从原来的表达矩阵里面随机抽取128个基因的表达矩阵做主成分分析，并且抽取一千次，每次主成分分析都可以得到第一主成分的贡献度值。那么，当我并不是随机抽取的时候，我就抽04380这条pathway的128个基因，也做主成分分析，并且计算得到第一主成分分析的重要性值。我们看看这个值，跟随机抽1000次得到的值差别大不大。
这时候就需要用到表达矩阵啦！
setwd("D:\\my_tutorial\\补\\用limma包对芯片数据做差异分析")
exprSet=read.table("GSE63067_series_matrix.txt.gz",comment.char = "!",stringsAsFactors=F,header=T)
rownames(exprSet)=exprSet[,1]
exprSet=exprSet[,-1]
我们根据ncbi里面对GSE63067的介绍可以知道，对应NASH和normal的样本的ID号，就可以提取我们需要的表达矩阵
把前面两属于Steatosis的样本去掉即可，exprSet=exprSet[,-c(1:2)]
然后再把芯片探针的id转换成entrez id
exprSet=exprSet[,-c(1:2)]
library(hgu133plus2.db)
library(annotate)
platformDB="hgu133plus2.db";
probeset <- rownames(exprSet)
rowMeans <- rowMeans(exprSet)
EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))
match_row=aggregate(rowMeans,by=list(EGID),max)
colnames(match_row)=c("EGID","rowMeans")
dat=data.frame(EGID,rowMeans,probeset)
tmp_prob=merge(dat,match_row,by=c("EGID","rowMeans"))
relevantProbesets=as.character(tmp_prob$probeset)
length(relevantProbesets) #hgu133plus2.db 20156
exprSet=exprSet[relevantProbesets,]
EGID_name=as.numeric(lookUp(relevantProbesets, platformDB, "ENTREZID"))
rownames(exprSet)=as.character(EGID_name)
d=exprSet
最后得到表达矩阵表格
我们首先得到1000次随机挑选128个基因的表达矩阵的主成分分析，第一主成分贡献度值。
gene128=sapply(1:1000,function(y) {
dat=t(d[sample(row.names(d), 128, replace=TRUE), ]);
round(100*summary(fast.prcomp(dat))$importance[2,1],2)
}
)
很快就能得到结果，可以看到数据如下
> summary(gene128)
Min. 1st Qu. Median Mean 3rd Qu. Max.
19.1 25.8 28.8 29.8 32.5 59.7
那么接下来我们挑选这个04380这条pathway特有128个基因来算第一主成分贡献度值
path_04380_gene=intersect(rownames(d),as.character(Path2GeneID[['04380']]))
dat=t(d[path_04380_gene,]);
round(100*summary(fast.prcomp(dat))$importance[2,1],2)
得到的值是38.83，然后看看我们的这个38.83在之前随机得到的1000个数里面是否正常，就按照正态分布检验来算
1-pnorm((38.83-mean(gene128))/sd(gene128))
[1] 0.0625
可以看到已经非常显著的不正常了，可以说明这条通路被富集了。
至少说明超几何分布检验的方法得到的富集分析结果跟我们这次的重抽样+主成分分析结果是一致的，当然，也有不一致的，不然就不用发明一种新的方法了。
如果写一个循环同样可以检验所有的通路，但是这样就不需要事先准备好差异基因啦！！！这是这个分析方法的特点！
主成分分析是为了简化变量的个数。
我这里不涉及到任何高级统计知识来简单讲解一下主成分分析，首先我们用下面的代码随机创造一个矩阵：
options(digits = 2)
x=c(rnorm(5),rnorm(5)+4)
y=3*c(rnorm(5),rnorm(5)+4)
dat=rbind(x,y,a=0.1*x,b=0.2*x,c=0.3*x,o=0.1*y,p=0.2*y,q=0.3*y)
colnames(dat)=paste('s',1:10,sep="")
dat
library(gmodels)
pca=fast.prcomp(t(dat))
pca
summary(pca)$importance
biplot(pca, cex=c(1.3, 1.2));
之前我提到过kegg数据库里面的有些pathway下面没有对应任何基因，当时我还在奇怪，怎么会有这样的通路呢？
然后，我随机挑选了其中一条通路（hsa01000），进行查看，发现正好是所有的酶的信息。
好奇怪，不明白为什么kegg要列出所有的酶信息
http://www.genome.jp/kegg-bin/get_htext?hsa01000
下载htext格式的酶的信息
798K Dec 15 2015 hsa01000.keg
查看文件，发现也是层级非常清楚的结构，D已经是最底级别的酶了，而E对应的基因是属于该酶的。
简单统计一下，发现跟酶相关的基因有3688个，而最底级别的酶有5811个，应该会持续更新的
如果想做成kegg那样的基因与酶对应表格，也是非常简单的！
国际系统分类法按酶促反应类型,将酶分成六个大类:
1、氧化还原酶类(oxidoreductases) 催化底物进行氧化还原反应的酶类。包括电子或氢的转移以及分子氧参加的反应。常见的有脱氢酶、氧化酶、还原酶和过氧化物酶等。
2、转移酶类(transferases) 催化底物进行某些基团转移或交换的酶类,如甲基转移酶、氨基转移酶、转硫酶等。
3、水解酶类(hydrolases)催化底物进行水解反应的酶类。如淀粉酶、粮糖苷酶、蛋白酶等。
4、裂解酶类(lyases)或裂合酶类(synthases) 催化底物通过非水解途径移去一个基团形成双键或其逆反应的酶类,如脱水酶、脱羧酸酶、醛缩酶等。如果催化底物进行逆反应,使其中一底物失去双键,两底物间形成新的化学键,此时为裂合酶类。
5、异构酶类(isomerases)催化各种同分异构体、几何异构体或光学异构体间相互转换的酶类。如异构酶、消旋酶等。
6、连接酶类（ligases）或合成酶类(synthetases)催化两分子底物连接成一个分子化合物的酶类。
上述六大类酶用EC(enzyme commission)加1.2.3.4.5.6编号表示,再按酶所催化的化学键和参加反应的基团,将酶大类再进一步分成亚类和亚-亚类,最后为该酶在这亚-亚类中的排序。如α淀粉酶的国际系统分类纩号为:EC3.2.1.1
EC3——Hydrolases 水解酶类
EC3.2——Glycosylases 转葡糖基酶亚类
EC3.2.l——Glycosidases 糖苷酶亚亚类i.e.enzymes hydmlyzing O-and S-glycosyl compound即能水解O-和S-糖基化合物
EC3.2.1.1 Alpha-amylase, α-淀粉酶
值得注意的是,即使是同一名称和EC编号,但来自不同的物种或不同的组织和细胞的同一种酸,如来自动物胰脏、麦芽等和枯草杆菌BF7658的α-淀粉酶等,它们的一级结构或反应机制可解不同,它们虽然都能催化淀扮的水解反应,但有不同的活力和最适合反应条件。
可以按照酶在国际分类编号或其推荐名,从酶手册(Enzyme Handbook)、酶数据库中检索到该酶的结构、特性、活力测定和Km值等有用信息。著名的手册和数据库有:
手册:
1、Schomburg,M.Salzmann and D.Stephan:Enzyme Handbook 10 Volumes
2、美国Worthington Biochemical Corporation:Enzyme Manual
(http://www.worthington-biochem.com/index/manual.htm/)
数据库:
l、德国BRENDA：Enzyme Database(http://www.brenda wnzymes.org)
2、Swissprot:EXPASYENZYME Enzyme nomenclature database (http://www.expasy.org/enzyme/)
3、IntEnz:Integrated relational Enzyme database (http://www.ebi.ac.uk/mtenz)
首先要明白，需要下载什么？
要下载四万多条GO记录的详细信息（http://purl.obolibrary.org/obo/go/go-basic.obo）
要下载GO与GO之间的关系（http://archive.geneontology.org/latest-termdb/go_daily-termdb-tables.tar.gz）
要下载GO与基因之间的对应关系！（物种）（ftp://ftp.ncbi.nlm.nih.gov/gene/DATA）
去官网！
http://geneontology.org/page/download-ontology
grep '\[Term\]' go-basic.obo |wc
43992 43992 307944
版本的区别!刚才我们下载的GO共有43992条，而以前的版本才38804条
GO与GO之间的关系
对应关系也在更新
> as.list(GOBPPARENTS['GO:0000042'])
$`GO:0000042`
is_a is_a is_a is_a
"GO:0000301" "GO:0006605" "GO:0016482" "GO:0072600"
library(org.Hs.eg.db)
library(GO.db)
> tmp=toTable(org.Hs.egGO) ##这个只包括基因与最直接的go的对应关系
> dim(tmp)
[1] 213101 4
> tmp2=toTable(org.Hs.egGO2ALLEGS) #这个是所有的基因与go的对应关系
> dim(tmp2)
[1] 2218968 4
基因与GO的对应关系也在更新
grep '^9606' gene2go |wc -l ### ##这个只包括基因与最直接的go的对应关系
269063
ftp://ftp.informatics.jax.org/pub/reports/index.html#go
ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
ftp://ftp.informatics.jax.org/pub/reports/index.html#go
其实主要要讲的不是用excel来做差异分析，只是想讲清楚差异分析的原理，用excel可视化的操作可能会更方便理解，而且想告诉大家，其实生物信息学分析，本来就很简单的，那么多软件，只有你理解了原理，你自己就能写出来的！
首先，还是得到表达矩阵，下面绿色的样本是NASH组，蓝色的样本是normal组
我们进行差异分析，很简单，就是看两组的表达值，是否差异，而检验的方法就是T检验。
=AVERAGE(D2:L2) ##求NASH组的平均表达量
=AVERAGE(M2:S2) ###求normal的平均表达量
=T2-U2 ##计算得到logFOLDchange值
=AVERAGE(D2:S2) ###得到所有样本的平均表达量
=T.TEST(D2:L2,M2:T2,2,3) ###用T检验得到两个组的表达量的差异显著程度。
简单检查几个值就可以看到跟limma包得到的结果差不多。
下载该R语言包，然后看说明书，需要自己做好三个数据（表达矩阵，分组矩阵，差异比较矩阵），总共三个步骤（lmFit,eBayes,topTable）就可以啦
首先做第一个数据，基因表达矩阵！
自己在NCBI里面可以查到下载地址，然后用R语言读取即可
exprSet=read.table("GSE63067_series_matrix.txt.gz",comment.char = "!",stringsAsFactors=F,header=T)
rownames(exprSet)=exprSet[,1]
exprSet=exprSet[,-1]
然后做好分组矩阵，如下
然后做好，差异比较矩阵，就是说明你想把那些组拿起来做差异分析，如下
最后输出结果：
我进行了6次比较，所以会输出6次比较结果
最后打开差异结果，解读，说明书如下！
在我的github有完整代码
打开官网：http://www.genome.jp/kegg-bin/get_htext?hsa00001+3101
http://www.genome.jp/kegg-bin/get_htext#A1 （这个好像打不开）
可以在里面找到下载链接
下载得到文本文件，可以看到里面的结构层次非常清楚，
C开头的就是kegg的pathway的ID所在行，D开头的就是属于它的kegg的所有的基因
A,B是kegg的分类，总共是6个大类，42个小类
grep ^A hsa00001.keg
A<b>Metabolism</b>
A<b>Genetic Information Processing</b>
A<b>Environmental Information Processing</b>
A<b>Cellular Processes</b>
A<b>Organismal Systems</b>
A<b>Human Diseases</b>
也可以看到，到目前为止（2015年12月8日20:26:57），共有343个kegg的pathway信息啦
接下来我们就把这个信息解析一下：
perl -alne '{if(/^C/){/PATH:hsa(\d+)/;$kegg=$1}else{print "$kegg\t$F[1]" if /^D/ and $kegg;}}' hsa00001.keg >kegg2gene.txt
这样就得到了
但是我发现了一个问题，有些通路竟然是没有基因的，我不是很明白为什么？
C 04030 G protein-coupled receptors [BR:hsa04030]
C 01020 Enzyme-linked receptors [BR:hsa01020]
C 04050 Cytokine receptors [BR:hsa04050]
C 03310 Nuclear receptors [BR:hsa03310]
C 04040 Ion channels [BR:hsa04040]
C 04031 GTP-binding proteins [BR:hsa04031]
那我们来看看kegg数据库更新的情况吧。
首先我们看org.Hs.eg.db这个R包里面自带的数据
Date for KEGG data: 2011-Mar15
org.Hs.egPATH has 5869 entrez genes and 229 pathways
2015年八月我用的时候是 6901 entrez genes and 295 pathways
现在是299个通路，6992个基因
所以这个更新其实很缓慢的，所以大家还在用DAVID这种网络工具做kegg的富集分析结果也差不大！
matlab毕竟是收费软件，而且是有界面的。所以搞生物信息的都用R和linux替代了，但是很多高大上的单位，比如大名鼎鼎的broadinstitute，是用matlab的，所以他们开发的程序也会以matlab代码的形式发布。但是考虑到大多研究者用不起matlab，或者不会用，所以就用linux系统里面安装matlab运行环境来解决这个问题，我们仍然可以把人家写的matlab程序，在linux命令行下面，当做一个脚本来运行！
比如，这次我就需要用broadinstitute的一个软件：Mutsig，找cancer driver gene的，http://www.broadinstitute.org/cancer/cga/mutsig_run，但是我看了说明才发现，它是用matlab写的，所以我要想在我的服务器用，就必须按照安装matlab运行环境，在官网可以下载：http://www.mathworks.com/products/compiler/mcr/
我这里选择的是R2013a (8.1)，下载之后解压是这样的，压缩包约四百多M
然后直接在解压后的目录里面运行那个install即可，然后如果你的linux可以传送图像，那么就会想安装windows软件一样方便！如果你的linux是纯粹的命令行，那么，就需要一步步的命令行交互，选择安装地址，等等来安装了。
记住你安装之后，会显示一些环境变量给你，请千万要记住，然后自己去修改自己的环境变量，如果你忘记了，就需要搜索来解决环境变量的问题啦！安装之后是这样的：
请记住你的安装目录，以后你运行其它matlab相关的程序，都需要把这个安装目录，当做一个参数传给你的其它程序的！！！
如果你没有设置环境变量，就会出各种各样的错误，用下面这个脚本可以设置
其中MCRROOT一般是$path/biosoft/matlab_running/v81/ 这样的东西，请务必注意，LD_LIBRARY_PATH非常重要，非常重要，非常重要！！！！
MCRROOT=$1
echo ---
LD_LIBRARY_PATH=.:${MCRROOT}/runtime/glnxa64 ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/bin/glnxa64 ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/sys/os/glnxa64;
MCRJRE=${MCRROOT}/sys/java/jre/glnxa64/jre/lib/amd64 ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/native_threads ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/server ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/client ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE} ;
XAPPLRESDIR=${MCRROOT}/X11/app-defaults ;
export LD_LIBRARY_PATH;
export XAPPLRESDIR;
echo LD_LIBRARY_PATH is ${LD_LIBRARY_PATH};
如果你没有设置正确，那么会报一下的错误！
error while loading shared libraries: libmwmclmcrrt.so.8.1: cannot o pen shared object file: No such file or directory
error while loading shared libraries: libmwlaunchermain.so: cannot o pen shared object file: No such file or directory
等等！！！
很多文献都提到外显子组的序列仅占全基因组序列的1%左右，但大多数与疾病相关的变异位于外显子区。通过外显子组测序可鉴定约8万个变异，全基因组测序可鉴定300万个变异，因此与全基因组测序相比，外显子组测序不仅费用较低，数据阐释也更为简单。外显子组测序技术以其经济有效的优势广泛应用于孟德尔遗传病、罕见综合征及复杂疾病的研究，并于2010年被Science杂志评为十大突破之一。所以我一直想验证一下外显子到底占基因组什么样的百分比，是哪些位点。
首先我们拿到外显子记录和基因信息，下面是通过R来从bioconductor中心提取数据
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
txdb
exon_txdb=exons(txdb)
genes_txdb=genes(txdb)
tmp=as.data.frame(exon_txdb)
write.table(tmp,'exon.pos',row.names=F)
tmp=as.data.frame(genes_txdb)
write.table(tmp,'gene.pos',row.names=F)
首先我们看看我刚才从TxDb.Hsapiens.UCSC.hg19.knownGene包里面提取的外显子记录文件exon.pos
我简单用脚本统计了一下：perl -alne '{$tmp+=$F[3]} END{print $tmp}' exon.pos
共289970个外显子记录，长度共98759283bp，也就是98.5Mb的区域都是外显子，感觉有点不大对，毕竟真正的外显子也就三十多M的，所以必须有些外显子是并列关系，也就是有的外显子包含着另外一些外显子，然后我写脚本检查了一下，的确太多了，这样的重复！
也可以去NCBI里面拿到 consensus coding sequence (CCDS)记录：ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/
那么我再看看NCBI里面拿到 consensus coding sequence (CCDS)记录CCDS.20150512.txt文件。
共33140行，一个基因一行，所以需要格式化成外显子文件，简单看了一下文件的列，很容易格式化。每一行的信息如下
1 NC_000001.11 SAMD11 148398 CCDS2.2 Public + 925941 944152 [925941-926012, 930154-930335, 931038-931088, 935771-935895, 939039-939128, 939274-939459, 941143-941305, 942135-942250, 942409-942487, 942558-943057, 943252-943376, 943697-943807, 943907-944152] Identical
用下面这个脚本格式化
格式化后的文件如下
外显子记录增加到了346493个，我再用单行命令计算长度，是56321786bp，这次只有56M了，还是有点大，所以应该是还有重复！
所以，我写了一个脚本来精确计算外显子的总长度，不能那样马马虎虎的用单行命令了。
这次才是34729283bp，也就是约35M，根很多文献里面说的都一样！
同理，我统计了一下外显子的侧翼长度
54692160 外显子加上前后50bp
73066288 外显子加上前后100bp
90362533 外显子加上前后150bp
而hg19版本基因组里面有着entrez gene ID号的基因是23056个基因，所以我接下来探究一下这些基因的信息！
我们首先看看基因与基因之间的交叉情况
其中有12454266bp的位点，是多个外显子共有的，可能是一个基因的多个外显子，或者是不同基因的外显子
我们主要看看不同基因的交叉情况
不同基因交叉情况共516种，共涉及498种基因！
ZNF341 NFS1
ZNF397 ZSCAN30
ZNF428 SRRM5
ZNF436-AS1 RPL11
ZNF578 ZNF137P
ZNF619 ZNF621
ZNF816 ZNF816-ZNF321P
ZSCAN26 ZSCAN31
ZSCAN5A ZNF542P
ZZEF1 ANKFY1
我随便在数据库里面搜索一下验证，发现这些基因的确是交叉重叠的
X-DAS-Version: DAS/0.95 X-DAS-Status: 200 Content-Type:text Access-Control-Allow-Origin: * Access-Control-Expose-Headers: X-DAS-Version X-DAS-Status X-DAS-Capabilities UCSC DAS Server. See http://www.biodas.org for more info on DAS. Try http://genome.ucsc.edu/cgi-bin/das/dsn for a list of databases. See our DAS FAQ (http://genome.ucsc.edu/FAQ/FAQdownloads#download23) for more information. Alternatively, we also provide query capability through our MySQL server; please see our FAQ for details (http://genome.ucsc.edu/FAQ/FAQdownloads#download29). Note that DAS is an inefficient protocol which does not support all types of annotation in our database. We recommend you access the UCSC database by downloading the tab-separated files in the downloads section (http://hgdownload.cse.ucsc.edu/downloads.html) or by using the Table Browser (http://genome.ucsc.edu/cgi-bin/hgTables) instead of DAS in most circumstances.
wget http://clavius.bc.edu/~erik/freebayes/freebayes-5d5b8ac0.tar.gz
tar xzvf freebayes-5d5b8ac0.tar.gz
cd freebayes
make
一个小插曲，安装的过程报错：/bin/sh: 1: cmake: not found
所以我需要自己下载安装cmake，然后把cmake添加到环境变量
首先下载源码包http://www.cmake.org/cmake/resources/software.html
wget http://cmake.org/files/v3.3/cmake-3.3.2.tar.gz
解压进去，然后源码安装三部曲,首先 ./configu 然后make 最后make install
cmake 会默认安装在 /usr/local/bin 下面，这里需要修改，因为你可能没有 /usr/local/bin 权限,安装到自己的目录，然后把它添加到环境变量！
step2：准备数据
wget http://bioinformatics.bc.edu/marthlab/download/gkno-cshl-2013/chr20.fawget ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/data/NA12878/alignment/NA12878.chrom20.ILLUMINA.bwa.CEU.low_coverage.20121211.bam
wget ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/data/NA12878/alignment/NA12878.chrom20.ILLUMINA.bwa.CEU.low_coverage.20121211.bam.bai用这个代码就可以下载千人基因组计划的NA12878样本的第20号染色体相关数据啦
freebayes -f chr20.fa \
NA12878.chrom20.ILLUMINA.bwa.CEU.low_coverage.20121211.bam >NA12878.chr20.freebayes.vcf一般就用默认参数即可
chr7 127471196 127472363 Pos1 0 + 127471196 127472363 255,0,0 chr7 127472363 127473530 Pos2 0 + 127472363 127473530 255,0,0
我的数据如下，需要自己创建成一个GRanges对象
library(GenomicRanges)
pfam.hg38 <- GRanges(seqnames=Rle(a[,2]),
ranges=IRanges(a[,3], a[,4]),
strand=a[,5])
library(rtracklayer)
ch = import.chain("hg38ToHg19.over.chain")
pfam.hg19 = liftOver(pfam.hg38, ch)
BLCA, BRCA, COADREAD, GBM, HNSC, KIRC, LAML, LGG, LUAD, LUSC, OV, PRAD, SKCM, STAD, THCA, UCEC.
Acute myeloid leukaemia
Bladder
Breast
Carcinoid
Chronic lymphocytic leukaemia
Colorectal
Diffuse large B-cell lymphoma
Endometrial
Oesophageal adenocarcinoma
Glioblastoma multiforme
Head and neck
Kidney clear cell
Lung adenocarcinoma
Lung squamous cell carcinoma
Medulloblastoma
Melanoma
Multiple myeloma
Neuroblastoma
Ovarian
Prostate
Rhabdoid tumour
HCD features: download
这是高置信度的癌症驱动基因列表：共280多个基因
Cancer5000 features: download
这是一篇对接近5000个癌症样本的研究得到的癌症相关基因列表：共230多个基因
参考：http://bg.upf.edu/oncodrive-role/
http://bioinformatics.oxfordjournals.org/content/30/17/i549.full
http://www.nature.com/nature/journal/v505/n7484/full/nature12912.html?WT.ec_id=NATURE-20140123
As of May 30, 2015, 201 different HPV types had been completely sequenced and officially recognized and divided into five PV-genera: Alpha-, Beta-, Gamma-, Mu-, and Nupapillomavirus.
根据文献，我找到了hpv所有已知测序种类的参考基因组网站：
http://www.hpvcenter.se/html/refclones.html
到目前（2015年7月31日15:17:59）已经有了205种，我爬取它们的genebank ID号，然后用python程序批量下载了它们的序列，能下载的序列共179条，都是8K左右的碱基序列。
根据genebank ID或者其它ID号批量下载核酸序列的脚本如下：
[python]</pre>
import sys
import time
import random
from Bio import Entrez
ids=[]
infile=sys.argv[1]
for line in open(infile,'r'):
line=line.strip()
ids.append(line)
for i in range(1,len(ids)):
# t = random.randrange(0,5)
handle =
Entrez.efetch(db="nucleotide", id=ids[i],rettype="fasta",email="jmzeng1314@163.com")
# time.sleep(t)
print handle.read()
[/python]
脚本使用很简单，保持输入文件是一行一个ID号即可。
同时，根据文献我们也能得到hbv病毒提取方法
当然，我是看不懂的。
同样的拿到下载的178条序列我们可以做一个进化树，当然，这个文章已经做好了，我就不做了，进化树其实蛮简单的。
下载179条hpv序列，每条序列都是8KB左右
我还用了R脚本批量下载
library(ape)
a=read.table("hpv_all.ID") #输入文件是一行一个ID号即可
for (i in 1:nrow(a)){
tmp=read.GenBank(a[i,1],seq.names = a[1,1],as.character = T)
write.dna(tmp,"tmp.fa",format="fasta", append=T,colsep = "")
}
然后用muscle做比对，参照我之前的笔记
http://www.bio-info-trainee.com/?p=659
http://www.bio-info-trainee.com/?p=660
http://www.bio-info-trainee.com/?p=626
muscle -in mouse_J.pro -out mouse_J.pro.a
muscle -maketree -in mouse_J.pro.a -out mouse_J.phy
貌似时间有点长呀，最后还莫名其妙的挂掉了，可能是我的服务器配置有点低。
进化树如下所示：
强烈推荐xshell安装，及破解网上找的一个注册码：101210-450789-147200
强烈editplus或者notepad++，搜索注册码是
注册名：Free User
注册码：6AC8D-784D8-DDZ95-B8W3A-45TFA
强烈建议大家玩一玩双系统制作ubuntu12.4版本的安装U盘，下载ultraISO，注册码为：王涛7C81-1689-4046-626F
肯定有很多人玩过了双系统就烦了，麻烦的是如何删除呢？
重点是需要修复MBR
1.进入win7，下载个软件MbrFix，放在C:\windows\system32文件夹中
2.点击开始》所有程序》附件》命令提示符
3.在命令提示符中输入MbrFix /drive 0 fixmbr /yes
4.此时已经修复了MBR，重启电脑后会发现已经没有了linux启动选项，直接进入windows了
这个需求很常见，因为双系统linux是可以看window的内容，但是window却无法直接查看linux内容。
下载ext2fsd可以把双系统里面linux分区的ext硬盘里面的文件通过window平台读取出来，软件下载地址如下
http://sourceforge.net/projects/ext2fsd/files/latest/download
