Monthly Archives: 12月 2015

十二 30

用GSEA来做基因集富集分析

Posted on by
how to use GSEA?
这个有点类似于pathway（GO,KEGG等）的富集分析，区别在于gene set（矫正好的基于文献的数据库）的概念更广泛一点，包括了

how to download GSEA ?

what's the input for the GSEA?

说明书上写的输入数据是：GSEA supported data files are simply tab delimited ASCII text files, which have special file extensions that identify them. For example, expression data usually has the extension *.gct, phenotypes *.cls, gene sets *.gmt, and chip annotations *.chip. Click the More on file formats help button to view detailed descriptions of all the data file formats.
实际上没那么复杂，一个表达矩阵即可！然后做一个分组说明的cls文件即可。
主要是自己看说明书，做出要求的数据格式：http://www.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats
表达矩阵我这里下载GSE1009数据集做测试吧！
cls的样本说明文件，就随便搞一搞吧，下面这个是例子：
6 2 1
# good bad
good good good bad bad bad
文件如下，六个样本，根据探针来的表达数据，分组前后各三个一组。
clipboard
现在开始运行GSEA！

start to run the GSEA !

首先载入数据
clipboard
确定无误，就开始运行，运行需要设置一定的参数！
clipboard

what's the output ?

输出的数据非常多，对你选择的gene set数据集里面的每个set都会分析看看是否符合富集的标准，富集就出来一个报告。
点击success就能进入报告主页，里面的链接可以进入任意一个分报告。
最大的特色是提供了大量的数据集：You can browse the MSigDB from the Molecular Signatures Database page of the GSEA web site or the Browse MSigDB page of the GSEA application. To browse the MSigDB from the GSEA application:
 
有些文献是基于GSEA的：

 

十二 29

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

Posted on by
该软件是broad institute写的一个数据接口，主要是供他人下载TCGA的所有寄存在broad institute的免费数据，主要是level3，level4的数据。(说错了，好像只有level4的数据，就是可以发文章的分析结果及图片)
软件下载地址：https://confluence.broadinstitute.org/display/GDAC/Download

懂它的使用规则，编码规则即可：
就是一个很简单的shell脚本而已，根据几个用户自定义参数来选择性的下载数据。
clipboard
我们可以用-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 by
无意中看一篇文献，有提到这个数据集，我简单看了一下，居然还真的这么全面！！！
它处理了目前（大概是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/
clipboard

 

十二 29

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

Posted on by
根据测试代码输出的报告如下：
clipboard
这些报告里面的元素，都是测试代码添加的，很容易就能理解
十二 29

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

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

官网：http://mutationassessor.org/

也应该是有standalone版本，我没有去找，不过网页版就很好用了，只需要复制粘贴进去自己想分析的数据，按照一定的格式即可，比如：
clipboard
该软件就能分析出该突变位点发生在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.

 

十二 25

使用进度条和并行计算来模拟计算赌博输赢概率

Posted on by

http://www.zhihu.com/question/19605599

以前看到过一个赌徒「必胜方法」
有一个广泛流传于赌徒中的“必胜方法”，是关于赌徒悖论极好的说明:一个赌场设有“猜大小”的赌局，玩家下注后猜“大”或者猜“小”，如果输了，则失去赌注，如果赢的话，则获得本金以及0.9倍的利润。 “必胜法”是这么玩的：
（1）押100猜“大”；
（2）如果赢的话，返回（1）继续；
（3）如果输的话则将赌注翻倍后继续猜“大”，因为不可能连续出现“大”，总会有获胜的时候，而且由于赌注一直是翻倍，只要赢一次，就会把所有输掉的钱赢回；
（4）只要赢了，就继续返回（1）。

一个朋友说写了程序来验证正确与否，我感觉蛮有趣的，想想自己也学了不少，也应该实践一下啦！
写出程序不难，就一个递归而已，不过是应该用并行计算，不然算的太慢了
另外一个重点是，如何随机？？？
[R]
judge=function(l=left,c=count,n=number,w=win){
c=c+1;
if(rnorm(1)>0){
l=l+2^n
n=1 # win
w=w+1
}else{
l=l-2^n
n=n+1 #lose
}
return(c(l,c,n,w))
}
my_fun=function(x){
tmp=judge(0,1,1,0)
for(i in 1:x){
tmp=judge(tmp[1],tmp[2],tmp[3],tmp[4])
print(tmp)
}
return(tmp[1])
}
my_fun(10)
[/R]
第一个函数是根据现在的状态，接受你的本金，第几次赌了，以及第几次连续下注，最后，你总共赢了多少次，这四个参数，然后随机给一个大小概率，这样你的本金会变化，赌的次数会加1
第二个函数是传入我们需要赌多少次，看结果：
> my_fun(10)
[1] -2  2  2  0
[1] 2 3 1 1
[1] 0 4 2 1
[1] -4  5  3  1
[1] 4 6 1 2
[1] 6 7 1 3
[1] 8 8 1 4
[1] 6 9 2 4
[1]  2 10  3  4
[1] 10 11  1  5
[1]  8 12  2  5
[1] 8
如果我们赌11次，可以看到，我最后会剩余8块钱，每次输赢的情况都反应在里面了，可以自己模拟多次看看！
因为我只赌了11次，所以很快，如果我赌1000次，而且还想检验一下10000次模拟结果，就会比较慢了！
我首先使用进度条模拟一下结果，代码如下：
2
还是比较慢的##Time difference of 1.861802 mins
我用了apply，好像时间是节省了一些，不过聊胜于无！
> library(pbapply)
> start=Sys.time()
> results=pbsapply(1:10000,function(i){my_fun(1000)})
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
> Sys.time()-start  ##9.909524 secs
Time difference of 1.301539 mins
那接下来我们分析一下模拟的结果吧
> summary(results)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -64580     974     998     985    1020    1110
可以看到，我们平均下来是会赢钱的，而且赢面很大，赢的次数非常多！！！
但是，我们看下面的图就知道，一个很诡异的结果！
而且，我这里用的模拟胜负，并不是很好
3
这其实还只是小批量模拟，如果要模拟百亿次，首先我的笔记本肯定不行，cpu太破了，如果用服务器，就需要用并行计算啦！
##下面是用多核并行计算的代码，大家有兴趣可以自己去玩一下
library(parallel)
cl.cores <- detectCores() 
cl <- makeCluster(cl.cores)
clusterExport(cl, "judge") 
start=Sys.time()
results=parSapply( 
  cl=cl, 
  1:10000,
  my_fun(1000)
)
Sys.time()-start  ##4.260994 secs
stopCluster(cl)
十二 25

做癌症研究一定要把这几十篇TCGA的大文章看完

Posted on by

都是发表在nature,cell还有新英格兰医学杂志上面的超级文章！每个文章附件都有一百多页，比博士论文还长，但是它们的分析套路其实都一样，都是那几种数据，包括WGS，WES，RNA-Seq，芯片表达量，miRNA表达量，甲基化数据，蛋白数据。分析过程也差不多，无法就是对癌症进行进一步的分类，癌症亚型，或者看看driver mutation，进一步解释癌症病变，转移，扩散机理，或者找标记物signature，辅助治疗等等，具体的要等我把这些文献看完了才能再进一步讲解，请做癌症研究方向的一定要把它们看完。

1

我已经下载完了，大家如果没有权限下载，就需要自己想办法啦！

image

非常值得大家阅读！！！

 

十二 24

使用R包cgdsr来下载TCGA的数据

Posted on by

前面我讲到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 →

十二 24

TCGA数据下载大全

Posted on by
并不是所有的数据都能下载，很多数据需要有权限才能下载的！！！
首先我们可以根据TCGA的文章来下载数据：

总共也就几十篇文章，都是发表在大杂志上面的。
每篇文章都会提供数据的打包下载，例如：

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

那个portal链接点击进去，就可以看到所有的下载链接了！
这是根据文章来分类打包好的数据！

同时也可以通过其它数据接口来下载

Tools for Exploring Data and Analyses

TCGA Data Portal

这几个接口都挺好用的：
非常详细！！！而且还专门写了软件接口：https://confluence.broadinstitute.org/display/GDAC/Download
或者写了R的接口：http://www.cbioportal.org/cgds_r.jsp
一般都推荐用TCGA自己的数据接口：https://tcga-data.nci.nih.gov/tcga/
里面对所有的样本都进行了统计
通过https://tcga-data.nci.nih.gov/tcga/dataAccessMatrix.htm可以进行定制化的数据下载！
里面有很多TCGA自定义的名词：

Data Levels and Data Types

https://tcga-data.nci.nih.gov/docs/dictionary/ 可以看到所有名词的解释：
数据的种类如下：
还记得以前看到一篇TCGA自己的关于胃癌的文章，发表在nature上面，文章涉及到了TCGA的各个方面的分析，所以附件PDF竟然有133页！！！

包含的其它数据有：
十二 23

看看Y染色体上面的基因在测序数据里的覆盖度和测序深度

Posted on by
作用：可以检测别人是否把自己的样本搞混，也可以看看测序是否分布均匀！
首先，我们要拿到Y染色体上面的基因的坐标信息！
因为我们的是hg19，所以我们要下载hg19的基因信息！
我们首先解析refGene文件，找到chrY的unique基因！
这四列分别是：chromosome/start/end/gene_symbol
clipboard
4程序如下：

[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]

对一个男性样本，结果会如下：
gene/length/pos
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
从结果可以看出来，很多基因都是y染色体特有的，这个结果是表明我们的测序非常棒

 

十二 22

根据chrY独有区域的覆盖度及测序深度来判断性别

Posted on by
这个也是基于bam文件来的，判断chrY独有区域的覆盖度及测序深度
首先下载chrY独有区域的记录文件，https://www.familytreedna.com/documents/bigy_targets.txt
然后用samtools depth来统计测序深度，samtools  depth $i |grep 'chr[XY]'
depth统计结果文件如下：
mzeng@ubuntu:/home/jmzeng/gender_determination$ head Sample3.depth
chrX    60085    1
chrX    60086    1
chrX    60087    1
chrX    60088    1
chrX    60089    1
chrX    60090    1
chrX    60091    1
chrX    60092    1
chrX    60093    1
chrX    60094    1
然后我随便写了一个脚本来对测序深度文件进行再统计，统计覆盖度及测序深度

[perl]
open FH,"bigy_targets.txt";
while(<FH>){
 chomp;
 @F=split;
 $all+=$F[2]-$F[1]+1;
 foreach ($F[1]..$F[2]){
  $h{$_}=1;
 }
}
close FH;
open FH,$ARGV[0];
while(<FH>){
 chomp;
 @F=split;
 next unless $F[0] eq 'chrY';
 if (exists $h{$F[1]}){
  $pos++;
  $depth+=$F[2];
 }
}
close FH;
$average=$depth/$pos;
$coverage=$pos/$all;
print "$pos\t$average\t$coverage\n" ;

[/perl]

 

 
这样对那三个样本结果如下：
clipboard
可以看到只有sample4，是覆盖率极低的，而且记录到的pos位点也特别少，所以她是女性！
这里测序深度没有意义。
十二 22

根据X,Y染色体比对上的reads数来判断性别

Posted on by
针对高通量测序数据，包括WGS，WES
我这里主要讲根据bam文件里面的chrX和chrY的reas比例来判断性别，大家可以自己处理数据得到bam文件！
主要是读取bam文件，选择chrX上面的记录，统计genotype即可：
也可以统计测序深度，samtools  depth $i |grep 'chr[XY]'
如果chrX和chrY的reads比例超过一定值，比如50倍，就判定为女性！
脚本很简单，就统计bam文件的第三列就可以啦，第三列就是染色体信息
 samtools  view $i |  perl -alne '{$h{$F[2]}++}END{print "$_\t$h{$_}" foreach sort keys %h}' >$out.chromosome.stat
对于Sample3，
chrX 1233061
chrY 140506
 
对于Sample4，很明显，这个是女性！
chrX 2734815
chrY 51860
 
对于Sample5
chrX 1329302
chrY 156663
十二 22

根据X染色体的snp的纯和率来判断性别

Posted on by
针对高通量测序数据，包括WGS，WES，甚至snp6芯片也行。
我这里主要讲根据vcf文件里面的chrX的纯合率来判断性别，大家可以自己处理数据得到vcf文件！
主要是读取vcf文件，选择chrX上面的记录，统计genotype即可：
clipboard
我这里拿之前的自闭症项目数据来举例子：
根据数据提供者的信息，3-4-5分别就是孩子、父亲、母亲，统计chrX的snp的的纯合和杂合的比例，代码很简单
vcf文件一般第一列是染色体，第6列是质量，第10列是基因型已经测序深度相关信息

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}' 

如果纯合的snp是杂合的倍数超过一定阈值，比如4倍，就可以判断是男性。

对于Sample3来说，很明显，是男孩，因为X染色体都是纯合突变

0/1 391
1/1 2463
2/2 1
1/2 32
0/2 1
对于Sample4来说，很明显，应该是母亲，证明之前别人给我的信息有误
1/1 3559
1/2 27
0/1 1835
0/2 5
那么Sample5很明显就是父亲咯
1/1 2626
0/1 356
1/2 22
十二 21

用重抽样+主成分方法来做富集分析

Posted on by

之前我们用超几何分布检验的方法做了富集分析，使用的是GSE63067.diffexp.NASH-normal.txt的logFC的绝对值大于0.5，并且P-value小雨0.05的基因作为差异基因来检验kegg的pathway的富集情况

结果是这样的

image001

我们接下来用另外一种方法来做富集分析，顺便检验一下，是不是超几何分布统计检验的富集分析方法就是最好的呢？

这种方法是-重抽样+主成分分析

大概的原理是，比如对上图中的，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

最后得到表达矩阵表格

image002

我们首先得到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

image003

 

那么接下来我们挑选这个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)

image004

得到的值是38.83，然后看看我们的这个38.83在之前随机得到的1000个数里面是否正常，就按照正态分布检验来算

1-pnorm((38.83-mean(gene128))/sd(gene128))

[1] 0.0625

可以看到已经非常显著的不正常了，可以说明这条通路被富集了。

至少说明超几何分布检验的方法得到的富集分析结果跟我们这次的重抽样+主成分分析结果是一致的，当然，也有不一致的，不然就不用发明一种新的方法了。

如果写一个循环同样可以检验所有的通路，但是这样就不需要事先准备好差异基因啦！！！这是这个分析方法的特点！

十二 21

主成分分析略讲

Posted on by

主成分分析是为了简化变量的个数。
我这里不涉及到任何高级统计知识来简单讲解一下主成分分析，首先我们用下面的代码随机创造一个矩阵：

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));

Continue reading →

十二 15

用超几何分布检验做富集分析

Posted on by

我们可以直接使用R的bioconductor里面的一个包，GOstats里面的函数来做超几何分布检验，看看每条pathway是否会富集

我们直接读取用limma包做好的差异分析结果

setwd("D:\\my_tutorial\\补\\用limma包对芯片数据做差异分析")

DEG=read.table("GSE63067.diffexp.NASH-normal.txt",stringsAsFactors = F)

View(DEG)

image001

我们挑选logFC的绝对值大于0.5，并且P-value小雨0.05的基因作为差异基因，并且转换成entrezID

probeset=rownames(DEG[abs(DEG[,1])>0.5 & DEG[,4]<0.05,])

library(hgu133plus2.db)

library(annotate)

platformDB="hgu133plus2.db";

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

length(unique(EGID))

#[1] 775

diff_gene_list <- unique(EGID)

这样我们的到来775个差异基因的一个list

首先我们直接使用R的bioconductor里面的一个包，GOstats里面的函数来做超几何分布检验，看看每条pathway是否会富集

library(GOstats)

library(org.Hs.eg.db)

#then do kegg pathway enrichment !

hyperG.params = new("KEGGHyperGParams", geneIds=diff_gene_list, universeGeneIds=NULL, annotation="org.Hs.eg.db",

categoryName="KEGG", pvalueCutoff=1, testDirection = "over")

KEGG.hyperG.results = hyperGTest(hyperG.params);

htmlReport(KEGG.hyperG.results, file="kegg.enrichment.html", summary.args=list("htmlLinks"=TRUE))

结果如下：

image002

但是这样我们就忽略了其中的原理，我们不知道这些数据是如何算出来的，只是由别人写好的包得到了结果罢了。

事实上，这个包的这个hyperGTest函数无法就是包装了一个超几何分布检验而已。

如果我们了解了其中的统计学原理，我们完全可以写成一个自建的函数来实现同样的功能。

超几何分布很简单，球分成黑白两色，数量已知，那么你随机抽有限个球，应该抽多少白球的问题！
公式就是 exp_count=n*M/N
然后你实际上抽了多少白球，就可以计算一个概率值！
换算成通路的富集概念就是，总共有多少基因，你的通路有多少基因，你的通路被抽中了多少基因（在差异基因里面属于你的通路的基因），这样的数据就足够算出上面表格里面所有的数据啦！
tmp=toTable(org.Hs.egPATH)
GeneID2Path=tapply(tmp[,2],as.factor(tmp[,1]),function(x) x)
Path2GeneID=tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)
#phyper(k-1,M, N-M, n, lower.tail=F)
#n*M/N
diff_gene_has_path=intersect(diff_gene_list,names(GeneID2Path))
n=length(diff_gene_has_path) #321 # 这里算出你总共抽取了多少个球
N=length(GeneID2Path) #5870  ##这里算出你总共有多少个球（这里是错的，有多少个球取决于背景基因！一般是两万个）
options(digits = 4)
for (i in names(Path2GeneID)){
 M=length(Path2GeneID[[i]])  ##这个算出你的所有的球里面，白球有多少个
 exp_count=n*M/N  ###这里算出你抽取的球里面应该多多少个是白色
 k=0         ##这个k是你实际上抽取了多少个白球
 for (j in diff_gene_has_path){
 if (i %in% GeneID2Path[[j]]) k=k+1
 }
 OddsRatio=k/exp_count
 p=phyper(k-1,M, N-M, n, lower.tail=F)  ##根据你实际上抽取的白球个数，就能算出富集概率啦！
 print(paste(i,p,OddsRatio,exp_count,k,M,sep="    "))
}
随便检查一下，就知道结果是一模一样的！

 

十二 15

下载所有的酶的信息，并且解析好

Posted on by

之前我提到过kegg数据库里面的有些pathway下面没有对应任何基因，当时我还在奇怪，怎么会有这样的通路呢？

然后，我随机挑选了其中一条通路（hsa01000），进行查看，发现正好是所有的酶的信息。

好奇怪，不明白为什么kegg要列出所有的酶信息

http://www.genome.jp/kegg-bin/get_htext?hsa01000

image001

下载htext格式的酶的信息

798K Dec 15  2015 hsa01000.keg

查看文件，发现也是层级非常清楚的结构，D已经是最底级别的酶了，而E对应的基因是属于该酶的。

简单统计一下，发现跟酶相关的基因有3688个，而最底级别的酶有5811个，应该会持续更新的

image002

如果想做成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)

 

十二 13

可怕的中国教育（转）

Posted on by
无意中看到这篇文章，感觉挺有趣的！
我不知道国外的教育是咋样的，但是文中描述的中国教育的试题都是千真万确的，我也经历过。
我没看过有钱人家的小孩，不知道那些花了几百万买个几平米的学区房的那些人的小孩是不是也要经受这样的教育摧残？
我还有个问题，我是如何出淤泥而不染的呢？
随便瞎扯，大家自己看这篇文章吧！
可怕的中国教育，聪明伶俐进去呆若木鸡出来
       作者：葛小琴
　　来源：杂文月刊（文摘版）
　　侄子在读高二，考了一道历史题：成吉思汗的继承人窝阔台，公元哪一年死?最远打到哪里?答不出来，我帮他查找资料，所以到现在我都记得，是打到现在的匈牙利附近。
　　在一次偶然的机会，我发现美国世界史这道题目不是这样考的。它的题目是这样的：成吉思汗的继承人窝阔台，当初如果没有死，欧洲会发生什么变化?试从经济、政治、社会三方面分析。
　 　有个学生是这样回答的：这位蒙古领导人如果当初没有死，那么可怕的黑死病，就不会被带到欧洲去，后来才知道那个东西是老鼠身上的跳蚤引起的鼠疫。但是六 百多年前，黑死病在欧洲猖獗的时候，谁晓得这个叫做鼠疫?如果没有黑死病，神父跟修女就不会死亡。神父跟修女如果没有死亡，就不会怀疑上帝的存在。如果没 有怀疑上帝的存在，就不会有意大利弗罗伦斯的文艺复兴?
　　如果没有文艺复兴，西班牙、南欧就不会强大，西班牙无敌舰队就不可能建立。如果西班牙、意大利不够强大，盎格鲁—撒克逊会提早200年强大，日耳曼会控制中欧，奥匈帝国就不可能存在。
　　教师一看“棒，分析得好。”但他们没有分数，只有等级A。其实这种题目老师是没有标准答案的，可是大家都要思考。
　　不久前，我去了趟日本，日本总是和我们在历史问题上产生纠葛，所以我在日本很注意高中生的教科书。
　 　他们的教师给高中生布置了这样一道题：日本跟中国100年打一次仗，19世纪打了日清战争(即甲午战争)，20世纪打了一场日中战争(即抗日战 争)，21世纪如果日本跟中国开火，你认为大概是什么时候?可能的远因和近因在哪里?如果日本赢了，是赢在什么地方?输了是输在什么条件上?分析之。
　 　其中有个高中生是这样分析的：我们跟中国很可能在台湾回到中国以后，有一场激战。台湾如果回到中国，中国会把基隆与高雄封锁，台湾海峡就会变成中国的内 海，我们的油轮就统统走右边，走基隆和高雄的右边。这样，会增加日本的运油成本。我们的石油从波斯湾出来跨过印度洋，穿过马六甲海峡，上中国南海，跨台湾 海峡进东海到日本海，这是石油生命线，中国政府如果把台湾海峡封锁起来，我们的货轮一定要从那里经过，我们的主力舰和驱逐舰就会出动，中国海军一看到日本 出兵，马上就会上场，就开打。
　　按照判断，公元2015年至2020年之间，这场战争可能爆发。所以，我们现在就要做对华抗战的准备。
　　我看其他学生的判断，也都是中国跟日本的摩擦会从东海从台湾海峡开始，时间判断是2015年至2020年之间。
　　这种题目和答案都太可怕了。
　 　撇开政治因素来看这道题，我们的历史教育就很有问题。翻开我们的教科书，题目是这样出的：甲午战争是哪一年爆发的？签订的叫什么条约？割让多少土地？赔 偿多少银两？每个学生都努力做答案。结果我们一天到晚研究什么时候割让辽东半岛，什么时候丢了台湾、澎湖、赔偿二万银两，1894年爆发甲午战争、 1895年签订马关条约，背得滚瓜烂熟，都是一大堆枯燥无味的数字。
　　那又怎么样，反正都赔了嘛，银两都给了嘛，最主要的是将来可能会怎样。
　　人家是在培养能力，而我们是在灌输知识，这是值得深思的部分！
　　看外面的教育，再看我们的教育？
　　老妈去参加我侄子的家长会，带回了两套侄子的考试试卷，我很好奇，拿过来看了现在小学生的试卷后，我震惊了！这是什么狗屁教育？这样的教育有希望吗？下面给大家详细说说我看到了什么。
　　侄子在本市某著名小学读书，有这么几道题。
　　一个春天的夜晚，一个久别家乡的人，望着皎洁的月光不禁思念起了故乡，于是吟起了一首诗：( )，( )。
　 　我看到侄子答的是：举头望明月，低头思故乡。但后面是一把大大的×，我就奇怪了，我也是想到的这两句。好奇地问侄子，这个不对??那答案是什么?侄子说 标准答案是：春风又绿江南岸，明月何时照我还。哎，这就奇怪了，因为是个春天的夜晚，就要是这句有春风的？？？要这个思念故乡的人不是江南的，是不可能说 出春风又绿江南岸这句话的!!!举头望明月，低头思故乡应该更准确。再扯远点，思念故乡，一千个人可以吟一千句不一样的诗，这个也可以有标准答案的么?
　　接下来是默写，题目是：我们学过《桂林山水》一文，请将下面句子默写下来。然后就是整段的要默写，这有什么用？死记硬背别人的文字有什么用?
　 　还有个题目，《匆匆》这篇课文，是现代着名作家朱自清先生写的，同学们都很喜欢这篇散文，你能把自己最喜欢，印象最深刻的一句写下来吗？我侄子写的是： 我的日子滴在时间的流里，没有声音，也没有影子。后面一把好大的×。标准答案竟然是：但是，聪明的你告诉我，我们的日子为什么一去不复返呢？这就更奇怪 了，一篇文章，你可以喜欢这句，我可以喜欢那句，难道最喜欢的一句话也要统一么?为什么“我的日子滴在时间的流里，没有声音，也没有影子。”这句不能喜 欢?就一定要喜欢“但是，聪明的你告诉我，我们的日子为什么一去不复返呢?”这句？我觉得这个题目应该是“你能把老师最喜欢，印象最深刻的一句写下来 吗?”才对!
　　再看别的试卷，更莫名其妙了，比如请说出阿拉伯数字的来历，是哪个国家创造的？侄子不知道，问我，我也不知道。我只好去搜一下，才知道是古印度人发明的。莫非我吃块猪肉，还一定得知道它是哪个养猪场养出来的?
　 　最后有个题目让我彻底崩溃了：请用一句话说明π的含义。侄子回答π的含义是圆周率。竟然打的是×，这就奇怪了，正好我老婆大学说读的是理科，我马上问 她，π是什么意思，她说圆周率啊。两个人狂汗，问了侄子半天，标准答案大概是，π是一个在数学及物理学领域普遍存在的数学常数……
　　如果你也觉得这种教育很无耻，就请转发吧，让更多的人来参与呼吁改变，为了孩子为了国家的未来 …… 别让孩子聪明伶俐地进去呆若木鸡地出来！
十二 12

下载最新版GO，并且解析好

Posted on by

首先要明白，需要下载什么？

要下载四万多条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

image001

grep '\[Term\]' go-basic.obo |wc

43992   43992  307944

版本的区别!刚才我们下载的GO共有43992条，而以前的版本才38804条

image002

GO与GO之间的关系

image003

对应关系也在更新

> as.list(GOBPPARENTS['GO:0000042'])

$`GO:0000042`

is_a         is_a         is_a         is_a

"GO:0000301" "GO:0006605" "GO:0016482" "GO:0072600"

image004

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

http://www.ebi.ac.uk/QuickGO

ftp://ftp.ncbi.nlm.nih.gov/gene/DATA

ftp://ftp.informatics.jax.org/pub/reports/index.html#go

 

 

十二 11

关于limma包差异分析结果的logFC解释

Posted on by

首先，我们要明白，limma接受的输入参数就是一个表达矩阵，而且是log后的表达矩阵（以2为底）。

那么最后计算得到的logFC这一列的值，其实就是输入的表达矩阵中case一组的平均表达量减去control一组的平均表达量的值，那么就会有正负之分，代表了case相当于control组来说，该基因是上调还是下调。

我之前总是有疑问，明明是case一组的平均表达量和control一组的平均表达量差值呀，跟log foldchange没有什么关系呀。

后来，我终于想通了，因为我们输入的是log后的表达矩阵，那么case一组的平均表达量和control一组的平均表达量都是log了的，那么它们的差值其实就是log的foldchange

首先，我们要理解foldchange的意义，如果case是平均表达量是8，control是2，那么foldchange就是4，logFC就是2咯

那么在limma包里面，输入的时候case的平均表达量被log后是3，control是1，那么差值是2，就是说logFC就是2。

这不是巧合，只是一个很简单的数学公式log(x/y)=log(x)-log(y)