六 24

Genomemapper软件使用说明书

Posted on by

 我以前一直以为有了bwa跟bowtie，没什么必要用其它的alignment软件，直到我碰到了高插入删除的helicos三代测序数据，我才发现，这个古董软件genomemapper居然大有用武之地了。

一．下载并且安装该软件

这是最新版本了

Release 0.4.4 2012-10-30 source code including documentation

Wget http://1001genomes.org/data/software/genomemapper/genomemapper_0.4.4/genomemapper-0.4.4.tar.gz

这个软件安装很简单，解压进入目录，make一下即可

image001

看到make完了之后就会多了两个软件，其中一个是用来构建参考基因组索引，一个用来比对的！

二．准备数据

既然是比对软件，那么肯定是一个参考基因组，一个测序的fastq原始文件咯

当然这个软件比较奇葩，它还支持Multi-FASTA, FASTQ2 or SHORE flat file format,

三、比对命令

这里要分两步走，首先是构建参考基因组的索引，然后才是比对

/home/jmzeng/bio-soft/genomemapper-0.4.4/gmindex \

-i BRCA1.fa -x BRCA1.idx -t BRCA1.meta

首先构建索引，种子长度就用默认的12即可，然后构建完索引如下。

image002

然后进行比对即可

/home/jmzeng/bio-soft/genomemapper-0.4.4/genomemapper \

-i BRCA1.fa -q SRR258835.fastq -M 4 -G 2 -E 4 -o mapped_reads.fl -u unmapped_reads.fl

成功比对的都输出到了mapped_reads.fl -这个文件，未比对上的在unmapped_reads.fl

我有12344条序列，成功比对的只有5276条，但是如果我用精确比对的算法，只有一千五百条是可以比对的，所以用这个允许4个mismatch和2个gap的比对算法，大大提高了比对率。

然后我修改了比对参数可以达到5605,5654,5696的提升。但是没有质的飞跃，估计本身我的这种helicos测序数据错误率就太可怕了。

四，输出结果解读

image004

这个是很规则的tab键分割的文本字符，我就不解读了，大家看readme

六 08

探究各个步骤对snp-calling的影响

Posted on by

做snp-calling时很多标准流程都会提到去除PCR重复这个步骤，但是这个步骤对找snp的影响到底有多大呢？这里我们来探究一下

 

去除PCR重复前 样本名 去除PCR重复后
   106082 BC1-1.snp 103829
   101443 BC1-2.snp 99500
   103937 BC2-1.snp 101833
   102979 BC2-2.snp 101022
   105876 BC3-1.snp 103562
   109168 BC3-2.snp 107052
   107155 BC4-1.snp 104894
   108335 BC4-2.snp 106031
   100236 BC5-1.snp 98417
   102322 BC5-2.snp 100395
   103466 BC6-1.snp 101405
   112940 BC6-2.snp 110611
   113166 BC7-1.snp 110948
   114038 BC7-2.snp 116090
   123670 PC1-1.snp 121697
   111402 PC1-2.snp 109389
   106917 PC2-1.snp 105149
   108724 PC2-2.snp 106776

 

可以看到去除pcr重复这个脚本对snp-calling的结果影响甚小，就是少了那么一千多个snp，脚本如下，我是用picard-tools进行的去除PCR重复，当然也可以用samtools来进行同样的步骤

[shell]

<b>for i in *.sorted.bam</b>

<b>do</b>

<b>echo $i</b>

<b>java  -Xmx120g  -jar /home/jmzeng/snp-calling/resources/apps/picard-tools-1.119/MarkDuplicates.jar \</b>

<b>CREATE_INDEX=true REMOVE_DUPLICATES=True \</b>

<b>ASSUME_SORTED=True VALIDATION_STRINGENCY=LENIENT METRICS_FILE=/dev/null \</b>

<b>INPUT=$i OUTPUT=${i%%.*}.sort.dedup.bam</b>

<b>done</b>

[/shell]

然后我们首先看看没有产生变化的那些snp信息的改变

head -50  ../rmdup/out/snp/BC1-1.snp  |tail |cut -f 1,2,8

chr1 17222 ADP=428;WT=0;HET=1;HOM=0;NC=0

chr1 17999 ADP=185;WT=0;HET=1;HOM=0;NC=0

chr1 18091 ADP=147;WT=0;HET=1;HOM=0;NC=0

chr1 18200 ADP=278;WT=0;HET=1;HOM=0;NC=0

chr1 24786 ADP=238;WT=0;HET=1;HOM=0;NC=0

chr1 25072 ADP=24;WT=0;HET=1;HOM=0;NC=0

chr1 29256 ADP=44;WT=0;HET=1;HOM=0;NC=0

chr1 29265 ADP=44;WT=0;HET=1;HOM=0;NC=0

chr1 29790 ADP=351;WT=0;HET=1;HOM=0;NC=0

chr1 29939 ADP=109;WT=0;HET=1;HOM=0;NC=0

head -50   BC1-1.snp  |tail |cut -f 1,2,8

chr1 17222 ADP=457;WT=0;HET=1;HOM=0;NC=0

chr1 17999 ADP=196;WT=0;HET=1;HOM=0;NC=0

chr1 18091 ADP=155;WT=0;HET=1;HOM=0;NC=0

chr1 18200 ADP=313;WT=0;HET=1;HOM=0;NC=0

chr1 24786 ADP=254;WT=0;HET=1;HOM=0;NC=0

chr1 25072 ADP=25;WT=0;HET=1;HOM=0;NC=0

chr1 29256 ADP=46;WT=0;HET=1;HOM=0;NC=0

chr1 29265 ADP=46;WT=0;HET=1;HOM=0;NC=0

chr1 29790 ADP=440;WT=0;HET=1;HOM=0;NC=0

chr1 29939 ADP=123;WT=0;HET=1;HOM=0;NC=0

可以看到，同一位点的snp仍然可以找到，仅仅是对测序深度产生了影响

 
然后我们再看看去除PCR重复这个步骤减少了的snp，在原snp里面是怎么样的

perl -alne '{$file++ if eof(ARGV);unless ($file){$hash{"$F[0]_$F[1]"}=1} else {print if not exists $hash{"$F[0]_$F[1]"} } }' ../rmdup/out/snp/BC1-1.snp BC1-1.snp |less

这个脚本就可以把去除PCR重复找到的snp位点在没有去除PCR重复的找到的snp文件里面过滤掉，查看那些去除PCR重复之前独有的snp

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

8.00    8.00   11.00   44.26   25.00 7966.00

图片1

 

可以看到被过滤的snp大多都是测序深度太低了的，如下面的例子

chr1 726325 a 9 CCC.ccc,^:, IEHGHHG/9

chr1 726325 a 5 C.c,^:, IGH/9

 

chr1 726338 g 16 TTT.ttt,,....,,, IHGI:9<HIIFIHC5H

chr1 726338 g 10 T.t,,...,, II:HIIFH5H

 

可以看到这一步还是很有用的，但是怎么说呢，因为最后对snp的过滤本来就包含了一个步骤是对snp的测序深度小于20的给过滤掉

 

但是也有个别的测序深度非常高的snp居然也是被去除PCR重复这个步骤给搞没了！很奇怪，我还在探索之中.

grep 13777 BC1-1.mpileup  |head

chr1 13777 G 263 ........,.C,,,,,.,,,.......,,,..,....,,......,.....c,........,,,,,,,..,...,,,,,.........,......C.......,,,,,,,,,,.....,,,,,,,.,,,..C,,,,,,CC,c,,,...C..,,,,cC.C..CC.CC,,cc,.C...C,,,,CCc,c,,,,,,,c,C.C.CC...C.cc,c...,C.CCcc...,CCC.C.CC..CCC..CC.c,cc,cc,,cc,C.,,^!.^6.^6.^6.^!, HIHIIIIEIEIHGIIIFIHIG?IIIIHIIHIFHIIHICIIIHIIGIEIIGIIIGHIIIIIIHIIHIHIIIIIIIHII1I?GHHHEHHIIEIEHIIEIIHHIIFIIIFHIHIIIIHIHIIHIIHHIIEIIIIIIHIIIIIIIIIG1HIIIIHIHIEHIHIHIIIIIIIIIIIHICIHIIIIIEIIIIHICIHGGIIIIIIIIEHIHIIIIIIHFIGGIIIIGIIIGICIIIHIIIIIIIIIIIHHHIIIIIHIIHDDII>>>>>

grep 13777 BC1-1.rmdup.mpileup  |head

chr1 13777 G 240 ........,.C,,,,,.,,,.......,,,..,....,,......,....c,......,,,,,,,..,...,,,,,.........,......C......,,,,,,,,,,.....,,,,,,,.,,,..C,,,,,,CC,c,,,...C..,,,,cC..CC.CC,cc,.C...C,,,,CCc,c,,,,,,,cC.C.C..C.c,c...,C.CCcc...,CC.C.CCC..C.c,cc,,c,.,,^!.^6.^6.^!, HIHIIIIEIEIHGIIIFIHIG?IIIIHIIHIFHIIHICIIIHIIGIEIIIIIHIIIIIHIIHIHIIIIIIIHII1I?GHHHEHHIIEIEHIIEIHHIIFIIIFHIHIIIIHIHIIHIIHHIIEIIIIIIHIIIIIIIIIG1HIIIIHIHIEHIHIIIIIIIIIIHICIHIIIIIEIIIIHICIHGGIIIIIIHIHIIIIIHFIGGIIIIGIIIGCIIIIIIIIIIHHIIIHIHDII>>>>

 

然后我再搜索了一些

chr8 43092928 . A T . PASS ADP=7966;WT=0;HET=1;HOM=0;NC=0 GT:GQ:SDP:DP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 0/1:255:7967:7966:6261:1663:20.9%:0E0:39:39:3647:2614:1224:439

chr8 43092908 . T C . PASS ADP=6968;WT=0;HET=1;HOM=0;NC=0 GT:GQ:SDP:DP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 0/1:255:7002:6968:5315:1537:22.06%:0E0:37:38:3022:2293:890:647

chr8 43092898 . T G . PASS ADP=6517;WT=0;HET=1;HOM=0;NC=0 GT:GQ:SDP:DP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 0/1:255:6517:6517:4580:1587:24.35%:0E0:38:38:2533:2047:920:667

chr7 100642950 . T C . PASS ADP=770;WT=0;HET=1;HOM=0;NC=0 GT:GQ:SDP:DP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 0/1:255:771:770:615:155:20.13%:3.9035E-51:38:38:277:338:65:90

终于发现规律啦！！！原来它们的突变率都略高于20%，在没有去处PCR重复之前，是高于snp的阈值的，但是去除PCR重复对该位点的突变率产生了影响，使之未能通过筛选。

 

六 01

Samtools无法同时得到mpileup格式的数据和bcftools格式的数据

Posted on by

 来自于： https://www.biostars.org/p/63429/

I'm using samtools mpileup and would like to generate both a pileup file and a vcf file as output. I can see how to generate one or the other, but not both (unless I run mpileup twice). I suspect I am missing something simple.

Specifically, calling mpileup with the -g or -u flag causes it to compute genotype likelihoods and output a bcf. Leaving these flags off just gives a pileup. Is there any way to get both, without redoing the work of producing the pileup file? Can I get samtools to generate the bcf _from_ the pileup file in some way? Generating the bcf from the bam file, when I already have the pileup, seems wasteful.

Thanks for any help!

我写了脚本来运行，才发现我居然需要两个重复的步骤来得到mpileup格式的数据和bcftools格式的数据，而这很明显的重复并且浪费时间的工作

for i in *sam

do

echo $i

samtools view -bS $i >${i%.*}.bam

samtools sort ${i%.*}.bam ${i%.*}.sorted

samtools index ${i%.*}.sorted.bam

samtools mpileup -f /home/jmzeng/ref-database/hg19.fa  ${i%.*}.sorted.bam  >${i%.*}.mpileup

samtools mpileup -guSDf  /home/jmzeng/ref-database/hg19.fa  ${i%.*}.sorted.bam  | bcftools view -cvNg - > ${i%.*}.vcf

Done

我想得到mpileup格式，是因为后续的varscan等软件需要这个文件来call snp

而得到bcftools格式可以直接用bcftools进行snp-calling

samtools mpileup 命令只有用了-g或者-u那么就只会输出bcf文件

如果想得到mpileup格式的数据，就只能用-f参数。

  • bcftools doesn't work on pileup format data. It works on bcf/vcf files.
  • samtools provides a script called sam2vcf.pl, which works on the output of "samtools pileup". However, this command is deserted in newer versions. The output of "samtools mpileup" does not satisfy the requirement of sam2vcf.pl. You can check the required pileup format on lines 95-99, which is different from output of "samtools mpileup".

 

五 29

用R语言的RCurl包结合XML包批量下载生信课件

Posted on by

首先是宾夕法尼亚州立大学（The Pennsylvania State University缩写PSU)的生信课件下载，这个生信不仅有课件，而且在中国的优酷视频网站里面还有全套授课视频，非常棒！

image001

课程主页是http://www.personal.psu.edu/iua1/courses/2013-BMMB-597D.html

可以看出所有的课件pdf链接都在这一个页面，所以是非常简单的代码！

下面是R代码：

library(XML)

library(RCurl)

library(dplyr)

psu_edu_url='http://www.personal.psu.edu/iua1/courses/2013-BMMB-597D.html';

wp=getURL(psu_edu_url)

base='http://www.personal.psu.edu/iua1/courses/file';

#pse_edu_links=getHTMLLinks(psu_edu_url)

psu_edu_links=getHTMLLinks(wp)

psu_edu_pdf=psu_edu_links[grepl(".pdf$",psu_edu_links,perl=T)]

for (pdf in psu_edu_pdf){

down_url=getRelativeURL(pdf,base)

filename=last(strsplit(pdf,"/")[[1]])

cat("Now we down the ",filename,"\n")

#pdf_file=getBinaryURL(down_url)

#FH=file(filename,"wb")

#writeBin(pdf_file,FH)

#close(FH)

download.file(down_url,filename)

}

因为这三十个课件都是接近于10M，所以下载还是蛮耗时间的

image003

其实R语言里面有这个down_url函数，可以直接下载download.file(down_url,filename)

然后我开始下载德国自由大学的生信课件，这次不同于宾夕法尼亚州立大学的区别是，课程主页里面是各个课题的链接，而pdf讲义在各个课题里面，所以我把pdf下载写成了一个函数对我们的课题进行批量处理

library(XML)

library(RCurl)

library(dplyr)

base="http://www.mi.fu-berlin.de/w/ABI/Genomics12";

down_pdf=function(url){

links=getHTMLLinks(url)

pdf_links=links[grepl(".pdf$",links,perl=T)]

for (pdf in pdf_links){

down_url=getRelativeURL(pdf,base)

filename=last(strsplit(pdf,"/")[[1]])

cat("Now we down the ",filename,"\n")

#pdf_file=getBinaryURL(down_url)

#FH=file(filename,"wb")

#writeBin(pdf_file,FH)

#close(FH)

download.file(down_url,filename)

}

}

down_pdf(base)

list_lecture= paste("http://www.mi.fu-berlin.de/w/ABI/GenomicsLecture",1:15,"Materials",sep="")

for ( url in list_lecture ){

cat("Now we process the ",url ,"\n")

try(down_pdf(url))

}

image005

同样也是很多pdf需要下载

接下来下载Minnesota大学的关于生物信息的教程的ppt合集

主页是： https://www.msi.umn.edu/tutorial-materials

 

这个网页里面有64篇pdf格式的ppt，还有几个压缩包，本来是准备写爬虫来爬去的，但是后来想了想有点麻烦，而且还不一定会看，反正也是玩玩

就用linux的命令行简单实现了这个爬虫功能。

curl https://www.msi.umn.edu/tutorial-materials >tmp.txt

perl -alne '{/(https.*?pdf)/;print $1 if $1}' tmp.txt >pdf.address

perl -alne '{/(https.*?txt)/;print $1 if $1}' tmp.txt

perl -alne '{/(https.*?zip)/;print $1 if $1}' tmp.txt >zip.address

wget -i pdf.address

wget -i pdf.zip

这样就可以啦！

 

用爬虫也就是几句话的事情，因为我已经写好啦下载函数，只需要换一个主页即可下载页面所有的pdf文件啦！

 

五 29

R语言对Ozone数据处理笔记

Posted on by

一．首先加载一些包，这样才能获得该ozone数据

数据集介绍：

Ozone数据集是一个三维数组，记录了24×24个空间网格内，从1995年1月到2000年12月，共72个时间点上，中美洲每月的平均臭氧水平。

前两维分别表示纬度和经度，第三维表示时间。

加载包的代码如下：

library("MASS", lib.loc="C:/Program Files/R/R-3.1.1/library")

library("ggplot2", lib.loc="C:/Program Files/R/R-3.1.1/library")

library("plyr", lib.loc="C:/Program Files/R/R-3.1.1/library")

library("dplyr", lib.loc="C:/Program Files/R/R-3.1.1/library")

library("reshape2", lib.loc="C:/Program Files/R/R-3.1.1/library")

二．我们首先简单看看第一个地点的72个月的臭氧水平变化图。

plot(1:72,ozone[1,1,],type="l")

box(lty = '1373', col = 'red')

grid(nx=NA,ny=NULL,lty=1,lwd=1,col="gray")

看起来还算是蛮有规律的。

image001

三．然后把这72个月的数据分成年份来画图

绘图第一种方式如下：

value <-ozone[1, 1, ]

plot(value[1:12],type="b",pch=19,lwd=2,xaxt="n",col="black",

xlab="month",ylab="value")

axis(1,at=1:12,labels=c("Jan", "Feb", "Mar", "Apr", "May",

"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))

lines(value[13:24],col="red",type="b",pch=19,lwd=2)

lines(value[25:36],col="orange",type="b",pch=19,lwd=2)

lines(value[37:48],col="purple",type="b",pch=19,lwd=2)

lines(value[49:60],col="blue",type="b",pch=19,lwd=2)

lines(value[61:72],col="green",type="b",pch=19,lwd=2)

legend("bottomright",legend=c("1995","1996","1997","1998","1999","2000"),

lty=1,lwd=2,pch=rep(19,6),col=c("black","red","orange","purple","blue","green"),

ncol=1,bty="n",cex=1.2,

text.col=c("black","red","orange","purple","blue","green"),inset=0.01)

是首先画第一年的，然后逐年添加一条线，然后画图例，算是蛮复杂的。

image003

还有一个简单的方法，就是用ggplot这个包来画。

values <-ozone[1, 1, ]

months=c("Jan", "Feb", "Mar", "Apr", "May",

"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

years=c("1995","1996","1997","1998","1999","2000")

dat=data.frame(month=rep(months,6),value=values,year=rep(years,each=12))

ggplot(dat,aes(x=month,y=value,group=year,colour=year))+geom_line()

image005

四：测试一下稳健回归模型

稳健回归是加权最小二乘回归，或称文艺最小二乘回归。

MASS 包中的 rlm命令提供了不同形式的稳健回归拟合方式。

回归分析就是用数理统计的方法，研究自然界中变量之间存在的非确定的相互依赖和制约关系，并把这种非确定的相互依赖和制约关系用数学表达式表达出来。其目的在于利用这些数学表达式以及对这些表达式的精度估计，对未知变量作出预测或检验其变化，为决策服务。

介绍几个线性回归（linear regression）中的术语：

残差（Residual）: 基于回归方程的预测值与观测值的差。

离群点（Outlier）: 线性回归（linear regression）中的离群点是指对应残差较大的观测值。也就是说，当某个观测值与基于回归方程的预测值相差较大时，该观测值即可视为离群点。 离群点的出现一般是因为样本自身较为特殊或者数据录入错误导致的，当然也可能是其他问题。

杠杆率（Leverage）: 当某个观测值所对应的预测值为极端值时，该观测值称为高杠杆率点。杠杆率衡量的是独立变量对自身均值的偏异程度。高杠杆率的观测值对于回归方程的参数有重大影响。

影响力点：（Influence）: 若某观测值的剔除与否，对回归方程的系数估计有显著相应，则该观测值是具有影响力的，称为影响力点。影响力是高杠杆率和离群情况引起的。

Cook距离（Cook's distance）: 综合了杠杆率信息和残差信息的统计量。

values <-ozone[1, 1, ]

month.abbr=c("Jan", "Feb", "Mar", "Apr", "May",

"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

month_72 <-factor(rep(month.abbr, length = 72), levels = month.abbr)

deseas1 <-rlm(value ~ month_72 -1)

summary(deseas1)

image007

五．对我们的24*24个地理位置的数据都做稳健回归分析

deseasf<-function(value) rlm(value ~ month_72 -1, maxit= 50)

models <-alply(ozone, 1:2, deseasf) #model是一个list，储存着24*24次稳健回归的结果

failed <-laply(models, function(x) !x$converged)

coefs<-laply(models, coef)#coefs是一个三维数组，记录了所有24×24个位置中每个位置的12个系数

dimnames(coefs)[[3]] <-month.abbr

names(dimnames(coefs))[3] <-"month"

deseas<-laply(models, resid) #deseas是一个三维数组，记录了所有24×24个位置中每个位置的72个残差

dimnames(deseas)[[3]] <-1:72

names(dimnames(deseas))[3] <-"time"

六．对我们的稳健回归系数的三维矩阵进行降维处理，方便画图

通过reshape包可以对三维数组进行降维

coefs_df<-melt(coefs)

head(coefs_df)

lat   long month   value

1 -21.2 -113.8   Jan 264.3964

2 -18.7 -113.8   Jan 261.3284

3 -16.2 -113.8   Jan 260.9643

4 -13.7 -113.8   Jan 258.9999

5 -11.2 -113.8   Jan 255.9999

6 -8.7 -113.8   Jan 254.9999

可以看到第三维的month成功被降维了

还可以通过plyr这个数据工厂包来进行降维

coefs_df<-ddply(coefs_df, .(lat, long), transform, avg= mean(value), std= value / max(value))

>head(coefs_df)   lat   long month   value     avg       std1 -21.2 -113.8   Jan 264.3964 268.6604 0.92770682 -21.2 -113.8   Feb 259.2036 268.6604 0.90948623 -21.2 -113.8   Mar 255.0000 268.6604 0.89473684 -21.2 -113.8   Apr 252.0052 268.6604 0.88422885 -21.2 -113.8   May 258.5089 268.6604 0.90704866 -21.2 -113.8   Jun 265.3387 268.6604 0.9310129

可以看到，不仅成功降维了，还添加了几个属性变量

 

七．最后对降维的coef系数数据画热图

 

coef_limits<-range(coefs_df$value)

coef_mid<-mean(coefs_df$value)

monthsurface<-function(mon)

{

df<-subset(coefs_df, month == mon)

qplot(long, lat, data = df, fill = value, geom="tile") +

scale_fill_gradient(limits = coef_limits,

low = "lightskyblue", high = "yellow")

}

pdf("ozone-animation.pdf", width = 8, height = 8)

l_ply(month.abbr, monthsurface, .print = TRUE)

dev.off()

会在当前R的工作目录下面看到一个pdf文件，里面储存着12个月的在24*24个地理位置的系数热图。

image009

五 28

org.Xx.eg.db系列包概述

Posted on by

在bioconductor的官网里面可以查到共有111个系列包，基本上跨越了我们常见的物种啦！

org.Xx.eg.db系列包介绍65

斑马鱼：Bioconductor - org.Dr.eg.db - /packages/release/data/annotation/html/org.Dr.eg.db.html

Details biocViews AnnotationData , Danio_rerio , OrgDb Version 3

拟南芥：Bioconductor - org.At.tair.db - /packages/release/data/annotation/html/org.At.tair.db.html

Details biocViews AnnotationData , Arabidopsis_thaliana , OrgDb Version 3

小鼠：Bioconductor - org.Mm.eg.db - /packages/release/data/annotation/html/org.Mm.eg.db.html

Details biocViews AnnotationData , Mus_musculus , OrgDb , mouseLLMappings Version 3

人类：Bioconductor - org.Hs.eg.db - /packages/release/data/annotation/html/org.Hs.eg.db.html

Details biocViews AnnotationData , Homo_sapiens , OrgDb , humanLLMappings Version 3

对这些系列包的函数都一样，包括以下几个：

columns(x)  keytypes(x)  keys(x, keytype, ...)  select(x, keys, columns, keytype, ...)  saveDb(x, file)  loadDb(file, dbType, dbPackage, ...)

 

 

这些包就是bioconductor已经做好的数据库，我们可以根据定义好的ID号来进行任意的基因转换，现在支持的信息有一下几种！

keytypes(org.Hs.eg.db)

[1] "ENTREZID"     "PFAM"         "IPI"          "PROSITE"      "ACCNUM"       "ALIAS"        "CHR"

[8] "CHRLOC"       "CHRLOCEND"    "ENZYME"       "MAP"          "PATH"         "PMID"         "REFSEQ"

[15] "SYMBOL"       "UNIGENE"      "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS" "GENENAME"     "UNIPROT"

[22] "GO"           "EVIDENCE"     "ONTOLOGY"     "GOALL"        "EVIDENCEALL"  "ONTOLOGYALL"  "OMIM"

[29] "UCSCKG"

这些包的确非常有用，大家可以看我博客里面关于它们的介绍！！！

 

五 28

菜鸟团第二次作业的部分答案

Posted on by

> library(org.Hs.eg.db)

载入需要的程辑包：AnnotationDbi载入需要的程辑包：stats4载入需要的程辑包：GenomeInfoDb载入需要的程辑包：S4Vectors载入需要的程辑包：IRanges载入程辑包：‘AnnotationDbi’The following object is masked from ‘package:GenomeInfoDb’:     species载入需要的程辑包：DBI

 

1、人共有多少个entrez id的基因呢？

x <- org.Hs.egENSEMBLTRANS

# Get the entrez gene IDs that are mapped to an Ensembl ID

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

length(x)

[1] 47721

可知共有47721个基因都是有entrez ID号的

2、能对应转录本ID的基因有多少个呢？

length(xx)

[1] 20592

可以看到共有20592个基因都是有转录本的！

2、能对应ensembl的gene ID的基因有多少个呢？

x <- org.Hs.egENSEMBL

# Get the entrez gene IDs that are mapped to an Ensembl ID

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

> length(x)

[1] 47721

> length(xx)

[1] 26019

可以看到只有26019是有ensembl的gene ID的

3、那么基因对应的转录本分布情况如何呢？

table(unlist(lapply(xx,length)))

菜鸟团第二次作业的部分答案863

可以看出绝大部分的基因都是20个转录本一下的，但也有极个别基因居然有高达两百个转录本，很可怕！

4、那么基因在染色体的分布情况如何呢？

x <- org.Hs.egCHR

# Get the entrez gene identifiers that are mapped to a chromosome

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

> length(x)

[1] 47721

> length(xx)

[1] 47232

可以看到有接近五百个基因居然是没有染色体定位信息的！！！

table(unlist(xx))

用barplot函数可视化一下，如图

 

菜鸟团第二次作业的部分答案1209

6、那么有多多少基因是有GO注释的呢？

x <- org.Hs.egGO

# Get the entrez gene identifiers that are mapped to a GO ID

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

length(xx)

[1] 18229

> length(x)

[1] 47721

可以看到只有18229个基因是有go注释信息的。

那么基因被注释的go的分布如何呢？

菜鸟团第二次作业的部分答案1477

可以看到大部分的基因都是只有30个go的，但是某些基因特别活跃，高达197个go注释。

还有kegg和omin数据库的我就不写了！

五 28

实战R语言bioconductor的seqinr包探究人的所有转录本的性质

Posted on by

首先安装这个包

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

biocLite("seqinr")

然后加载包，并读取我们的CDS.fa文件

library("seqinr")

human_cds=read.fasta("CDS.fa")

#这一个步骤非常耗时间，可能是因为我们的转录本文件有十万多个转录本的原因吧

str(human_cds) #查看可知读入了一个list，其中每个转录本都是list的一个元素

List of 100778

$ ENST00000415118:Class 'SeqFastadna'  atomic [1:8] g a a a ...

.. ..- attr(*, "name")= chr "ENST00000415118"

.. ..- attr(*, "Annot")= chr ">ENST00000415118 havana_ig_gene:known chromosome:GRCh38:14:22438547:22438554:1 gene:ENSG00000223997 gene_biotype:TR_D_gene tran"| __truncated__

$ ENST00000448914:Class 'SeqFastadna'  atomic [1:13] a c t g ...

.. ..- attr(*, "name")= chr "ENST00000448914"

.. ..- attr(*, "Annot")= chr ">ENST00000448914 havana_ig_gene:known chromosome:GRCh38:14:22449113:22449125:1 gene:ENSG00000228985 gene_biotype:TR_D_gene tran"| __truncated__

对list的每个元素都有几种函数可以处理得到信息：

Length,table,GC,count

其中count函数很有趣，数一数序列里面的这些组合出现的次数

count(dengueseq, 1)

count(dengueseq, 2)接下来我们随机取human_cds这个list的一个元素用这几个函数对它处理一下

> tmp=human_cds[[1]]

> tmp

[1] "g" "a" "a" "a" "t" "a" "g" "t"

attr(,"name")

[1] "ENST00000415118"

attr(,"Annot")

[1] ">ENST00000415118 havana_ig_gene:known chromosome:GRCh38:14:22438547:22438554:1 gene:ENSG00000223997 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene"

attr(,"class")

[1] "SeqFastadna"

再看看函数的结果

> length(tmp)

[1] 8

> table(tmp)

tmp

a g t

4 2 2

> GC(tmp)

[1] 0.25

> count(tmp,1)

 

a c g t

4 0 2 2

> count(tmp,2)

 

aa ac ag at ca cc cg ct ga gc gg gt ta tc tg tt

2  0  1  1  0  0  0  0  1  0  0  1  1  0  0  0

>

还是挺好用的，接下来我们应用R的知识来对着十万多个转录本进行一些简单的总结

human_cds_length=unlist(lapply(human_cds,length))

human_cds_gc=unlist(lapply(human_cds,GC))

这样就得到了所有转录本的长度和GC含量信息

然后我们简单统计一下，并画几个图表吧！

> summary(human_cds_length)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

3     366     699    1132    1425  108000

> summary(human_cds_gc)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

0.1467  0.4577  0.5285  0.5264  0.5932  0.8917

可以看到还是有很多很极端的转录本的存在！

最长的转录本也不过10k，而我记得最长的基因高达8M，看了内含子远大于外显子呀。

但是GC含量有很多高于80%，这些基因在二代测序的研究中是一个盲区。

实战R语言bioconductor的seqinr包探究人的所有转录本的性质2075

这些极端基因可以通过biomaRt等包进行注释，得到gene名和功能信息。

 

hist(human_cds_gc)

hist(log10(human_cds_length))

GC含量分布如图

实战R语言bioconductor的seqinr包探究人的所有转录本的性质2177

长度分布如图

实战R语言bioconductor的seqinr包探究人的所有转录本的性质2186

附表：

http://www.bioinformatics.org/sms/iupac.html 所有字符的碱基氨基酸意义表格

 

五 25

用R画GO注释二级分类统计图

Posted on by

群里有朋友问这个图怎么画，我想了想，这肯定是ggplot完成的，非常简单，但是菜鸟们缺乏实践，可能会困惑，所以我模拟数据画了一个！

图片1

首先构造数据

dat=data.frame(name=LETTERS[1:21],
 number=abs(rnorm(21)*10),
 type=c(rep("BP",7),rep("CC",7),rep("MF",7))
)
# 请务必自己查看dat是一个什么数据，print出来即可
# 然后对这个数据画图，一行代码即可！！！
library(ggplot2)
ggplot(dat,aes(x=name,y=number,fill=type))+geom_bar(stat="identity")+coord_flip()

看起来是不是很像回事啦！细节我就懒得调控啦！

图片2

其实自己搜索即可！坐标轴和主题都是可以控制的

http://rstudio-pubs-static.s3.amazonaws.com/3364_d1a578f521174152b46b19d0c83cbe7e.html

http://docs.ggplot2.org/0.9.3.1/coord_flip.html

​
五 21

Biostrings包简介

Posted on by

首先讲讲它的对象

有下面几个字符串对象BString, DNAString, RNAString and AAString可以通过以下代码构造它们：

b <- BString("I am a BString object")

d <- DNAString("TTGAAAA-CTC-N")

这两个对象的区别是DNAstring对象对字符串的要求严格一些，只有IUPAC字符和+-字符可以。

对构造好的对象可以通过下标来取子字符串对象，也可以通过subseq来取，但是子字符串仍然是数据对象，只有通过toString函数才能把它们转化成字符串。

用length(dd2)和nchar(toString(dd2))都可以找到我们Biostrings对象的长度。但是后者速度会很慢。

Views(RNAString("AU"), start=0, end=2)这个函数可以把string对象任意截取成list

start, end and width可以作用于我们截取的list，判断list里面的元素在原来的string对象上面的起始终止及长度信息。

 

接下来讲这个包带有的一个比对函数！

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede")

Global PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] succ--eed subject: [1] supersede score: -33.99738

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",type = "local")

Local PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] su subject: [1] su score: 5.578203

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",gapOpening = 0, gapExtension = 1)

Global PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] su-cce--ed subject: [1] sup--

可以看出这个比对函数可以调整的参数实在是太多了，而且改变参数之后比对情况大不一样，还有很多参数就不一一细讲了。

这个比对结果可以赋值给一个变量，保存比对的对象。

psa1 <- pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede")

class(psa1)

summary(psa1)

class(pattern(psa1))

class(summary(psa1))

score(psa2)

还可以自己构建打分矩阵来进行比对。

submat <-

+ matrix(-1, nrow = 26, ncol = 26, dimnames = list(letters, letters))

diag(submat) <- 0

Biostrings包简介1454

psa2 <-pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",substitutionMatrix = submat, gapOpening = 0, gapExtension = 1)

我们的包还自带了两个非常流行的氨基酸比对矩阵PAM和BLOSUM

ls("package:Biostrings")可以查看这个包所有的对象。

data(package="Biostrings")可以查看这个包所有的数据对象

还有很多其它函数

还可以去除adaptor，挺好玩的

既然有配对比对函数，那么就有多重比对函数！

我们可以读取clustaW, Phylip and Stolkholm这几种不同的比对结果文件来构造多重比对对象。

library(Biostrings)这个包里面自带了两个文件，我们可以示范一下构建对象。

origMAlign <- readDNAMultipleAlignment(filepath = system.file("extdata", "msx2_mRNA.aln", package="Biostrings"), format="clustal")

phylipMAlign <- readAAMultipleAlignment(filepath = system.file("extdata","Phylip.txt", package="Biostrings"),format="phylip")

 

对构造好的多重比对对象就可以构建进化树啦，代码如下！

sdist <- stringDist(as(origMAlign,"DNAStringSet"), method="hamming")

> clust <- hclust(sdist, method = "single")

> pdf(file="badTree.pdf")

> plot(clust)

> dev.off()

Biostrings包简介2345

五 21

Bioconductor的DO.db包介绍

Posted on by

Bioconductor的包都是同样的安装方法：

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

还有GO.bd包是完全一模一样的规则！！！

加载这个包可以发现它依赖于好几个其它的包，这也是我比较喜欢R的原因，它会自动把它需要的包全部安装加载进来，不需要自己一个个调试！

> library(DO.db)

载入需要的程辑包：AnnotationDbi

载入需要的程辑包：stats4

载入需要的程辑包：GenomeInfoDb

载入需要的程辑包：S4Vectors

载入需要的程辑包：IRanges

载入需要的程辑包：DBI

> help(DO.db)

> ls("package:DO.db")

[1] "DO"          "DO_dbconn"   "DO_dbfile"   "DO_dbInfo"   "DO_dbschema" "DOANCESTOR"  "DOCHILDREN"  "DOID"        "DOMAPCOUNTS"

[10] "DOOBSOLETE"  "DOOFFSPRING" "DOPARENTS"   "DOSYNONYM"   "DOTERM"      "DOTerms"     "Secondary"   "show"        "Synonym"

[19] "Term"

这个包里面有19个数据对象！都是比较高级的S4对象。

比如我们可以拿DOTERM[1:10]这个小的数据对象来做例子！example=DOTERM[1:10]

因为example是一个高级对象，所以无法直接查看，需要用as.list方法来查看

> as.list(example)

$`DOID:0001816`DOID: DOID:0001816Term: angiosarcomaSynonym: DOID:267Synonym: DOID:4508Synonym: "hemangiosarcoma" EXACT []Secondary: DOID:267Secondary: DOID:4508

~~~~~~~~~~~~共十个DO条目

对每一个DO条目来说都有DOID,Term,Synony这些函数可以取对应的值。

下面是对DO的有向无环图的数据解读

xx <- as.list(DOANCESTOR)可以查看每个DO与它所对应的上级条目DO，每个DO都会有不止一个的上级DO。

xx <- as.list(DOPARENTS)可以查看每个DO与它所对应的父条目DO，每个DO都有且只有一个父DO。

xx <- as.list(DOOFFSPRING)可以查看每个DO与它所对应的下级DO的关系列表，大多数DO都不止一个子条目DO，所有的下级DO都会列出。

xx <- as.list(DOCHILDREN)以查看每个DO与它所对应的子条目DO的关系列表，大多数DO都不止一个子条目DO。

还有Lkeys(DOTERM)可以查看数据库里面的所有的DO条目的ID号

> head(keys(DOTERM))

[1] "DOID:0000000" "DOID:0001816" "DOID:0002116" "DOID:0014667" "DOID:0050004" "DOID:0050012"

dbmeta(GO_dbconn(), "GOSOURCEDATE")

可以查看这个DO库的制备时间

> dbmeta(DO_dbconn(), "DOSOURCEDATE")

[1] "20140417"

五 21

RNA-seq流程对基因和转录本的表达量的计算

Posted on by

bedtools multicov和htseq-count都可以用来对基因和转录本的表达量的计算！！！

我们总共有四个样本，已经比对到小鼠的mm9基因组上面了，数据大小如下

RNA-seq流程对基因和转录本的表达量的计算111

然后对基因和转录本计数需要一些额外的信息，即各个基因及转录本的位置信息，gtf文件需要在UCSC等各大数据库下载

RNA-seq流程对基因和转录本的表达量的计算170

然后我们制作一个config文件配置我们的数据地址

cat sample_bam.config 可以看到文件内容如下

/data/mouse/ptan/740WT1.bam

/data/mouse/ptan/741WT2.bam

/data/mouse/ptan/742KO1.bam

/data/mouse/ptan/743KO2.bam

几个批处理文件名及内容分别如下

bedtools_multicov.sh  bedtools_multicov_transcript.sh  htseq.sh  htseq_transcript.sh

 

while read id

do

echo $id

new=`echo $id |cut -d"/" -f 5`

echo $new

bedtools multicov -bams $id -bed /data/mouse/mouse_mm9_gene.bed  > $new.gene.bedtools_multicov.count

done <$1

 

while read id

do

echo $id

new=`echo $id |cut -d"/" -f 5`

echo $new

bedtools multicov -bams $id -bed /data/mouse/mouse_mm9_transcript.bed  > $new.transcript.bedtools_multicov.count

done <$1

 

while read id

do

echo $id

new=`echo $id |cut -d"/" -f 5`

htseq-count -f bam $id /data/mouse/Mus_musculus.NCBIM37.67.gtf.chr  > $new.gene.htseq.count

done <$1

 

while read id

do

echo $id

new=`echo $id |cut -d"/" -f 5`

htseq-count -f bam --idattr transcript_id $id /data/mouse/Mus_musculus.NCBIM37.67.gtf.chr  > $new.transcript.htseq.count

done <$1

 

批量运行这些程序后就能对它们分别分情况进行计数，也能比较这两种计数方法的区别！

RNA-seq流程对基因和转录本的表达量的计算1201

 

可以看出区别还是很大的！！！

RNA-seq流程对基因和转录本的表达量的计算1219

我肯定没搞懂它们的原理，这完全就不一样，已经不是区别的问题了！！！

对于每个个体输出的计数文件，接下来就可以用DESeq等包来进行差异基因分析啦！

五 21

SAMStat软件使用说明书

Posted on by

这个软件是对我们的比对结果（通常是bwa,bowtie,tophat,hisat,star）bam或者sam来进行一个可视化的总结，类似于fastqc对我们的fastq测序结果做一个可视化总结，非常好用。

一．下载并安装该软件

软件主页是http://samstat.sourceforge.net/ 里面对该软件进行非常详细的说明

包括installation和usage，我这里简单的翻译一下。

Wget http://liquidtelecom.dl.sourceforge.net/project/samstat/samstat-1.5.tar.gz

解压开看里面的readme有介绍如何安装这个软件

Unpack the tarball:

bash-3.1$ tar -zxvf samstat-XXX.tar.gz

bash-3.1$ cd samstat

bash-3.1$ ./configure

bash-3.1$ make

bash-3.1$ make check

bash-3.1$ make install

如果用root命令就可以直接用samstat啦

如果没有root权限，安装的时候稍微有点不同

./configure  --prefix=/home/jmzeng/my-bin/

make

make install

很简单的

二，数据，就是我们的bam文件啦

三，运行命令

SAMStat软件使用说明书710

四，结果

简单看看samtools flagstat 740WT1.bam  的结果

19232378 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 duplicates

18846845 + 0 mapped (98.00%:-nan%)

0 + 0 paired in sequencing

0 + 0 read1

0 + 0 read2

0 + 0 properly paired (-nan%:-nan%)

0 + 0 with itself and mate mapped

0 + 0 singletons (-nan%:-nan%)

0 + 0 with mate mapped to a different chr

0 + 0 with mate mapped to a different chr (mapQ>=5)

然后再看看我们的samstat的结果！

740WT1.bam.samstat.html

一个网页，非常丰富的内容

SAMStat软件使用说明书1168

内容太多了，我懒得解释了

见软件说明书http://davetang.org/wiki/tiki-index.php?page=SAMStat

五 21

使用Bedtools对RNA-seq进行基因计数

Posted on by

以前是没有想过用这个软件的，直到有一个我的htseq无法对比对的bam文件进行基因计数（后来我才发现htseq无法计数的原因是gtf版本不同导致坐标不同，而且gtf对染色体编号没有加上chr），我简单搜索了一下，发现bedtools multicov也有类似的功能，所以我选择它来试试看！

首先注意它需要sort的bam文件及bam的index

bedtools multicov depends upon index BAM files in order to count the number of overlaps in each BAM file. As such, each BAM file should be position sorted (samtool sort aln.bam aln.sort) and indexed (samtools index aln.sort.bam) with either samtools or bamtools.

首先安装它：

wget https://github.com/arq5x/bedtools2/releases/download/v2.23.0/bedtools-2.23.0.tar.gz

解压开后

Make clean

Make all

然后就可以看到它的bin目录下全部是程序啦

Bedtools使用笔记639

命令很简单的

bedtools multicov [OPTIONS] -bams BAM1 BAM2 BAM3 ... BAMn -bed  <BED/GFF/VCF>

By default, multicov will report the count of alignments in each input BAM file that overlap.

例子：

bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed ivls-of-interest.bed

ivls-of-interest.bed这个文件是必须的，可能需要自己制作，其实用gtf文件也可以的

chr1 0   10000   ivl1

chr1 10000   20000   ivl2

chr1 20000   30000   ivl3

chr1 30000   40000   ivl4

输出结果前三列是坐标，第四列是基因名，跟我们的bed文件一样，只是最后三列是三个样本的计数，是添加上来的！

chr1 0       10000   ivl1    100 2234    0

chr1 10000   20000   ivl2    123 3245    1000

chr1 20000   30000   ivl3    213 2332    2034

chr1 30000   40000   ivl4    335 7654    0

 

同样是对gene的reads计数，bedtools的multicov工具与htseq-count的区别是

i'd guess it's due to the fact that htseq-count only reports one hit per aligned read assuming that read is aligned uniquely and does not overlap multiple features listed in your GTF. if an aligned read hits more than one feature in your GTF then it doesn't report that hit. bedtools gives you raw hits which includes every 1 hit for every intersection of every alignment with any features in the GTF no matter how many times it aligned or how many features it hit. you might think, "wow, htseq-count is dropping a lot of information". yes, it is! i've moved to using other tools to count hits to genes (RSEM/eXpress) since they disambiguate those ambiguous alignments and as a result you get counts for all of your aligned reads. in a genome with alternative splicing you lose too much data using htseq-count, in my opinion.

而且专门有个文献在讨论这个问题

http://barcwiki.wi.mit.edu/wiki/SOPs/rna-seq-diff-expressions

http://www.nature.com/nbt/journal/v31/n1/abs/nbt.2450.html

Differential analysis of gene regulation at transcript resolution with RNA-seq

下面我讲一个实际的例子

我的bam文件如下

Bedtools使用笔记2406

bedtools multicov -bams 740WT1.bam 741WT2.bam 742KO1.bam 743KO2.bam -bed mm9.bed

Bedtools使用笔记2491

得到的这个矩阵就可以去用DESeq包来进行差异分析啦！

五 18

R语言DESeq找差异基因

Posted on by

一：安装并加装该R包

安装就用source("http://bioconductor.org/biocLite.R") ;biocLite("DESeq")即可，如果安装失败，就需要自己下载源码包，然后安装R模块。

 

二．所需要数据

它的说明书指定了我们一个数据

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

安装了pasilla这个包之后，在这个包的安装目录就可以找到一个表格文件，就是我们的DESeq需要的文件。

C:\Program Files\R\R-3.2.0\library\pasilla\extdata\pasilla_gene_counts.tsv

说明书原话是这样的

The table cell in the i-th row and the j-th column of the table tells how many reads have been mapped to gene i in sample j.

一般我们需要用htseq-count这个程序对我们的每个样本的sam文件做处理计数，并合并这样的数据

下面这个是示例数据，第一列是基因ID号，后面的每一列都是一个样本。

图片1

de = newCountDataSet( pasillaCountTable, condition )  #根据我们的样本因子把基因计数表格读入成一个cds对象，这个newCountDataSet函数就是为了构建对象！

对我们构建好的de对象就可以直接开始找差异啦！非常简单的几步即可

de=estimateSizeFactors(de)

de=estimateDispersions(de)

res=nbinomTest(de,"K","W") #最重要的就是这个res表格啦！

uniq=na.omit(res)

我这里是对4个样本用htseq计数后的文件来做的，贴出完整代码吧

library(DESeq)

#首先读取htseq对bam或者sam比对文件的计数结果

K1=read.table("742KO1.count",row.names=1)

K2=read.table("743KO2.count",row.names=1)

W1=read.table("740WT1.count",row.names=1)

W2=read.table("741WT2.count",row.names=1)

data=cbind(K1,K2,W1,W2)

data=data[-c(43630:43634),]

#把我们的多个样本计数结果合并起来成数据框，列是不同样本，行是不同基因

colnames(data)=c("K1","K2","W1","W2")

type=rep(c("K","W"),c(2,2))

#构造成DESeq的对象，并对分组样本进行基因表达量检验

de=newCountDataSet(data,type)

de=estimateSizeFactors(de)

de=estimateDispersions(de)

res=nbinomTest(de,"K","W")

#res就是我们的表达量检验结果

library(org.Mm.eg.db)

tmp=select(org.Mm.eg.db, keys=res$id, columns="GO", keytype="ENSEMBL")

ensembl_go=unlist(tapply(tmp[,2],as.factor(tmp[,1]),function(x) paste(x,collapse ="|"),simplify =F))

#首先输出所有的计数数据，加上go注释信息

all_res=res

res$go=ensembl_go[res$id]

write.csv(res,file="all_data.csv",row.names =F)

#然后输出有意义的数据，即剔除那些没有检测到表达的基因

uniq=na.omit(res)

sort_uniq=uniq[order(uniq$padj),]

write.csv(sort_uniq,file="sort_uniq.csv",row.names =F)

#然后挑选出padj值小于0.05的差异基因数据来做富集，富集用的YGC的两个包，在我前面的博客已经详细说明了！

tmp=select(org.Mm.eg.db, keys=sort_uniq[sort_uniq$padj<0.05,1], columns="ENTREZID", keytype="ENSEMBL")

diff_ENTREZID=tmp$ENTREZID

require(DOSE)

require(clusterProfiler)

diff_ENTREZID=na.omit(diff_ENTREZID)

ego <- enrichGO(gene=diff_ENTREZID,organism="mouse",ont="CC",pvalueCutoff=0.01,readable=TRUE)

ekk <- enrichKEGG(gene=diff_ENTREZID,organism="mouse",pvalueCutoff=0.01,readable=TRUE)

write.csv(summary(ekk),"KEGG-enrich.csv",row.names =F)

write.csv(summary(ego),"GO-enrich.csv",row.names =F)

 

五 10

Perl的电子书打包下载！

Posted on by

我看到总是有人问我perl在生信方面的书，这里推荐一下我以前看到的一个帖子，因为书籍比较大，我就没有下载，这样它们在115网盘里面，下载可能会有一点麻烦的，但是它们都是mobi格式的，可以转换成epub格式，这样所有的电子书阅读器都可以阅读，绝对不是扫描版，全文字版！

我的网盘里面是精选的，可能会更有用一点

http://club.topsage.com/thread-2818785-1-1.html

里面有很多生物信息学的perl书籍

打包下载地址：
http://115.com/file/c2h0mgg2#
mobi_collection_2012-04-16(86_books).7z

以下是书目录

Perl/Advanced Perl Programming, 2nd Edition.mobi
Perl/Automating System Administration with Perl, 2nd Edition.mobi
Perl/Beginning Perl for Bioinformatics.mobi
Perl/Beginning Perl Web Development - From Novice to Professional.tpz
Perl/CGI Programming with Perl.mobi
Perl/Effective Perl Programming.mobi
Perl/Higher-Order Perl.mobi
Perl/Intermediate Perl.mobi
Perl/Learning Perl, 5th Edition.mobi
Perl/Learning Perl, 6th Edition.mobi
Perl/Mastering Algorithms with Perl.mobi
Perl/Mastering Perl for Bioinformatics.mobi
Perl/Mastering Perl.mobi
Perl/Mastering PerlTk - Graphical User Interfaces in Perl.mobi
Perl/Perl & LWP.mobi
Perl/Perl Best Practices.mobi
Perl/Perl by Example, 4th Edition.mobi
Perl/Perl Cookbook.mobi
Perl/Perl For Dummies.mobi
Perl/Perl Hacks - Tips & Tools for Programming, Debugging, and Surviving.mobi
Perl/Perl Pocket Reference.mobi
Perl/Perl Testing - A Developer's Notebook.mobi
Perl/Practical Text Mining with Perl.mobi
Perl/Pro Perl Parsing.tpz
Perl/Pro Perl.tpz
Perl/Programming Perl - Unmatched power for text processing and scripting, 4th Edition.mobi
Perl/Programming Perl, 3rd Edition.mobi
Perl/Programming the Perl DBI - Database programming with Perl.mobi
Perl/Randal Schwartz's Perls of Wisdom.tpz
Perl/The Definitive Guide to Catalyst - Writing Extensible, Scalable and Maintainable Perl Based Web Applications.mobi
Advanced Programming in the UNIX Environment, 2nd Edition.mobi
Agile Project Management with Scrum.mobi
Apache Security.mobi
Backup & Recovery.mobi
Beginning iOS 5 Development - Exploring the iOS SDK.mobi
Beginning Rails 3 (Expert's Voice in Web Development).mobi
Clean Code - A Handbook of Agile Software Craftsmanship.mobi
Core Java, Volume I--Fundamentals, 8th Edition.mobi
Core Java, Volume II--Advanced Features, 8th Edition.mobi
Cpp Primer, 4th Edition.mobi
Dreamweaver CS5.5 Mobile and Web Development with HTML5, CSS3, and jQuery.mobi
Erlang Programming.mobi
Essential System Administration, 3rd Edition.mobi
Even Faster Web Sites.mobi
Hacking Vim 7.2.mobi
Hacking-The Art of Exploitation, 2nd Edition.mobi
Hacking-The Next Generation.mobi
High Performance MySQL, 3rd Edition.mobi
High Performance Web Sites.mobi
IPv6 Essentials.mobi
Java The Complete Reference, 8th Edition.mobi
JavaScript - The Definitive Guide, 6th Edition.mobi
Joomla! Programming.mobi
Land of Lisp.mobi
Learning GNU Emacs, 3rd Edition.mobi
Learning PHP, MySQL, and JavaScript.mobi
Learning Python, 3rd Edition.mobi
Learning the bash Shell, 3rd Edition.mobi
Learning the vi and Vim Editors.mobi
Linux Kernel Development, 3rd Edition.mobi
Mac OS X for Unix Geeks.mobi
Managing Projects with GNU Make, 3rd Edition.mobi
Mastering Regular Expressions.mobi
Mercurial The Definitive Guide.mobi
MySQL High Availability.mobi
MySQL Stored Procedure Programming.mobi
Nagios System and Network Monitoring, 2nd Edition.mobi
PhoneGap Beginner's Guide.mobi
PHP and MySQL Web Development, 4th Edition.mobi
Pro jQuery Mobile.mobi
Programming in C, 3rd Edition.mobi
Programming in Objective-C, 4th Edition.mobi
Programming iOS 4 - Fundamentals of iPhone, iPad, and iPod touch Development.mobi
Python for Unix and Linux System Administration.mobi
Real World Haskell.mobi
Ruby on Rails 3 Tutorial - Learn Rails by Example.mobi
Sed & awk, 2nd Edition.mobi
Steve Jobs.mobi
The Art of SEO, 2nd Edition.mobi
The Rails 3 Way, 2nd Edition.mobi
The Ruby Programming Language.mobi
Time Management for System Administrators.mobi
Unix and Linux System Administration Handbook, 4th Edition.mobi
UNIX Power Tools, 3rd Edition.mobi
vi and Vim Editors Pocket Reference.mobi
Website Optimization.mobi

五 10

RNA-seq比对软件HISAT说明书

Posted on by

取代bowtie+tophat进行RNA-seq比对

HISAT全称为Hierarchical Indexing for Spliced Alignment of Transcripts，由约翰霍普金斯大学开发。它取代Bowtie/TopHat程序，能够将RNA-Seq的读取与基因组进行快速比对。这项成果发表在3月9日的《Nature Methods》上。

HISAT利用大量FM索引，以覆盖整个基因组。以人类基因组为例，它需要48,000个索引，每个索引代表~64,000 bp的基因组区域。这些小的索引结合几种比对策略，实现了RNA-Seq读取的高效比对，特别是那些跨越多个外显子的读取。尽管它利用大量索引，但HISAT只需要4.3 GB的内存。这种应用程序支持任何规模的基因组，包括那些超过40亿个碱基的。

HISAT软件可从以下地址获取：http://ccb.jhu.edu/software/hisat/index.shtml。

首先，我们安装这个软件！

Wget http://ccb.jhu.edu/software/hisat/downloads/hisat-0.1.5-beta-source.zip

官网下载的是源码包，需要make一下，make之后目录下面就多了很多程序，绿色的那些都是，看起来是不是很眼熟呀！！！

哈哈，这完全就是bowtie的模拟版本！！！

HISAT取代bowtie+tophat进行RNA-seq比对1222

也可以从github里面下载，wget https://codeload.github.com/infphilo/hisat/zip/master

下载后直接解压即可使用啦。当然这个软件本身也有着详尽的说明书

http://ccb.jhu.edu/software/hisat/manual.shtml

然后就是准备数据，它跟tophat一样的功能。就是把用RNA-seq方法测序得到的fastq文件比对到参考基因组上面，所以就准这两个文件了哦

接下来是运行程序！

说明书上面写着分成两个步骤，构建索引和比对。

这个软件包模仿bowtie自带了一个example数据，而且它的说明书也是针对于那个example来的，我也简单运行一下。

$HISAT_HOME/hisat-build $HISAT_HOME/example/reference/22_20-21M.fa 22_20-21M_hisat

构建索引的命令如上，跟bowtie一样我修改了一下

/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat-build 22_20-21M.fa  my_hisat_index

连日志都跟bowtie一模一样，哈哈，可以看到我们的这个参考fasta文件 22_20-21M.fa 就变成索引文件啦，索引还是很多的！

HISAT取代bowtie+tophat进行RNA-seq比对1871

然后就是比对咯，还是跟bowtie一样

$HISAT_HOME/hisat -x 22_20-21M_hisat -U $HISAT_HOME/example/reads/reads_1.fq -S eg1.sam

我的命令是

/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat -x  my_hisat_index -U ../reads/reads_1.fq  -S reads1.sam

1000 reads; of these:

1000 (100.00%) were unpaired; of these:

0 (0.00%) aligned 0 times

1000 (100.00%) aligned exactly 1 time

0 (0.00%) aligned >1 times

100.00% overall alignment rate

哈哈，到这里。这个软件就运行完毕啦！！！是不是非常简单，只有你会用bowtie，这个就没有问题。当然啦，软件还是有很多细节是需要调整的。我下面就简单讲一个实际的例子哈！

首先，我用了1.5小时把4.6G的小鼠基因组构建了索引

/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat-build  Mus_musculus.GRCm38.fa.fa mouse_hisat_index

HISAT取代bowtie+tophat进行RNA-seq比对2512

然后对我的四个测序文件进行比对。

for i in *fq

do

/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat  -x  /home/jmzeng/hoston/mouse/mouse_hisat_index  \

-p 30 -U  $i.trimmed.single  -S ./hisat_out/${i%.*}.sam

done

它运行的速度的确要比tophat快好多，太可怕的速度！！！！至于是否多消耗了内存我就没有看了

4.6G的小鼠，5G的测序数据，我只用了五个核，居然十分钟就跑完了！

然后听群友说是因为没有加 --known-splicesite-infile <path>这个参数的原因，没有用gtf文件来指导我们的RNA数据的比对，这样是不对的！

需要用下面这个脚本把gtf文件处理一下，然后导入什么那个参数来指导RNA比对。

extract_splice_sites.py genes.gtf > splicesites.txt

但是我报错了，错误很奇怪，没解决，但是我换了个 extract_splice_sites.py  程序，就可以运行啦！之前是HISAT 0.1.5-beta release 2/25/2015里面的python程序，后来我换做了github里面的就可以啦！

/home/jmzeng/hoston/RNA-soft/hisat-master/extract_splice_sites.py Mus_musculus.GRCm38.79.gtf >mouse_splicesites.txt

HISAT取代bowtie+tophat进行RNA-seq比对3218

21192819 reads; of these:
21192819 (100.00%) were unpaired; of these:
14236834 (67.18%) aligned 0 times
5437800 (25.66%) aligned exactly 1 time
1518185 (7.16%) aligned >1 times

感觉没有变化，不知道为什么？

21192819 reads; of these:

21192819 (100.00%) were unpaired; of these:

14236838 (67.18%) aligned 0 times

5437793 (25.66%) aligned exactly 1 time

1518188 (7.16%) aligned >1 times

32.82% overall alignment rate

发表这个软件的文献本身也把这个软件跟其它软件做了详尽的对比

http://www.nature.com/nmeth/journal/v12/n4/full/nmeth.3317.html

Program Run time (min) Memory usage (GB)
Run times and memory usage for HISAT and other spliced aligners to align 109 million 101-bp RNA-seq reads from a lung fibroblast data set. We used three CPU cores to run the programs on a Mac Pro with a 3.7 GHz Quad-Core Intel Xeon E5 processor and 64 GB of RAM.
HISATx1 22.7 4.3
HISATx2 47.7 4.3
HISAT 26.7 4.3
STAR 25 28
STARx2 50.5 28
GSNAP 291.9 20.2
OLego 989.5 3.7
TopHat2 1,170 4.3

 

 

 

参考：http://www.plob.org/2015/03/20/8980.html

http://nextgenseek.com/2015/03/hisat-a-fast-and-memory-lean-rna-seq-aligner/

 

五 10

RNA-seq的比对软件star说明书

Posted on by

类似于tophat的软件

首先当然是下载软件啦！

两个地方可以下载，一个是谷歌code中心，被墙啦，另一个是github，我的最爱。

wget https://codeload.github.com/alexdobin/STAR/zip/master

2013-RNA-seq的star-类似于tophat的软件说明书217

解压即可使用啦，其中程序在bin目录下面，根据自己的平台调用即可！

然后doc里面还有个pdf的说明文档，写的非常清楚，我也是看着那个文档学的这个软件！

接下来就是准备数据啦！

既然是类似于tophat一样的比对软件，当然是准备参考基因组和测序数据咯，毫无悬念。

然后 该软件也给出了一些测试数据

ftp://ftp2.cshl.edu/gingeraslab/tracks/STARrelease/2.1.4/

然后就是运行程序的命令！

分为两个步骤：首先构建索引，然后比对即可，中间的参数根据具体需要可以细调！

构建索引时候，软件说明书给的例子是

The basic options to generate genome indices are as follows:
--runThreadN NumberOfThreads
--runMode genomeGenerate
--genomeDir /path/to/genomeDir
--genomeFastaFiles /path/to/genome/fasta1 /path/to/genome/fasta2 ...
--sjdbGTFfile /path/to/annotations.gtf
--sjdbOverhang ReadLength-1

我模仿了一下。对我从ensembl ftp里面下载的老鼠基因组构建了索引

/home/jmzeng/hoston/RNA-soft/STAR-master/bin/Linux_x86_64/STAR  \

--runThreadN 30 #我的服务器还比较大，可以使用30个CPU  \

--runMode genomeGenerate \

--genomeDir  /home/jmzeng/hoston/mouse/STAR-mouse #构建好的索引放在这个目录 \

--genomeFastaFiles  /home/jmzeng/hoston/mouse/Mus_musculus.GRCm38.fa.fa \

--sjdbGTFfile /home/jmzeng/hoston/mouse/Mus_musculus.GRCm38.79.gtf \

--sjdbOverhang 284   #我的测序数据长短不一，最长的是285bp

当然注释的地方你要删除掉才行呀，因为cpu用的比较多。

2013-RNA-seq的star-类似于tophat的软件说明书1302

算一算时间，对4.6G的小鼠基因组来说，半个小时算是非常快的了！Bowtie2的index要搞两个多小时。

然后就是比对咯。这也是很简单的，软件说明书给的例子是

The basic options to run a mapping job are as follows:
--runThreadN NumberOfThreads
--genomeDir /path/to/genomeDir
--readFilesIn /path/to/read1 [/path/to/read2]

我稍微理解了一下参数，然后写出了自己的命令。

fq=740WT1.fq.trimmed.single

mkdir  740WT1_star

/home/jmzeng/hoston/RNA-soft/STAR-master/bin/Linux_x86_64/STAR  \

--runThreadN 20  \

--genomeDir  /home/jmzeng/hoston/mouse/STAR-mouse   \

--readFilesIn  $fq \

--outFileNamePrefix  ./740WT1_star/740WT1

如果输出文件需要被cufflinks套装软件继续使用。就需要用一下参数

Cufflinks/Cuffdiff require spliced alignments with XS strand attribute, which STAR will generate with --outSAMstrandField intronMotif option.

还有--outSAMtype参数可以修改输出比对文件格式，可以是sam也可以是bam，可以是sort好的，也可以是不sort的。

最后是输出文件解读咯！

其实没什么好解读的，输出反正就是sam类似的比对文件咯，如果还有其它文件，需要自己好好解读说明书啦。我就不废话了！

 

值得一提的是，该程序提供了2次map的建议

The basic idea is to run 1st pass of STAR mapping with the usual parameters , then collect the junctions detected in the first pass, and use them as ”annotated” junctions for the 2nd pass mapping.

在对RNA-seq做snp-calling的时候可以用到，尤其是GATK官方还给出了教程，大家可以好好学习学习！

http://www.broadinstitute.org/gatk/guide/article?id=3891

 

 

五 08

GWAS研究现状及资源下载

Posted on by

 GWAS研究是非常火的，NHGIR还专门为它开辟了专栏来介绍，下面这个图片也是来自于NHGIR组织，是GWAS近年来发表文章的状况。

图片1

可以在该文章上面下载这个所有的数据

wget http://www.genome.gov/admin/gwascatalog.txt

截至目前为止。2015年5月8日21:08:34

这个文档有19603行的数据，但是只有2113篇pubmed文献，共涉及到七千多个基因

有293种杂志都发过GWAS的文章，总共有2113篇文献，发表关联分析突变位点最多的是这篇文献23251661 在PLoS One杂志上面，共 949个rs突变

杂志排序
cut -f 2,5 gwascatalog.txt |perl -alne '{$hash{$_}++}END{print "$_" foreach sort {$hash{$a} <=> $hash{$b}} keys %hash}' |cut -f 2 |perl -alne '{$hash{$_}++}END{print "$_\t$hash{$_}" foreach sort {$hash{$a} <=> $hash{$b}} keys %hash}'
Hum Genet 41
Am J Hum Genet 62
Mol Psychiatry 64
PLoS One 132
PLoS Genet 145
Hum Mol Genet 168
Nat Genet 397
文章的rs突变点排序
cut -f 2,5 gwascatalog.txt |perl -alne '{$hash{$_}++}END{print "$_ $hash{$_}" foreach sort {$hash{$a} <=> $hash{$b}} keys %hash}'
24324551 PLoS One 241
24097068 Nat Genet 245
24816252 Nat Genet 299
23382691 PLoS Genet 699
23251661 PLoS One 949

数据打开如下：

我取了表头和第一行数据，然后把它转置了，这样方便查看

Date Added to Catalog 10/22/2014
PUBMEDID 24528284
First Author Ji Y
Date 08/01/2014
Journal Br J Clin Pharmacol
Link http://www.ncbi.nlm.nih.gov/pubmed/24528284
Study Citalopram and escitalopram plasma drug and metabolite concentrations: genome-wide associations.
Disease/Trait Response to serotonin reuptake inhibitors in major depressive disorder (plasma drug and metabolite levels)
Initial Sample Description 300 European ancestry Escitalpram treated individuals, 130 European ancestry Citalopram treated individuals
Replication Sample Description NA
Region 17q25.3
Chr_id 17
Chr_pos 79831041
Reported Gene(s) CBX4
Mapped_gene CBX8 - CBX4
Upstream_gene_id 57332
Downstream_gene_id 8535
Snp_gene_ids
Upstream_gene_distance 33.93
Downstream_gene_distance 2.12
Strongest SNP-Risk Allele rs9747992-?
SNPs rs9747992
Merged 0
Snp_id_current 9747992
Context Intergenic
Intergenic 1
Risk Allele Frequency 0.086
p-Value 2.00E-07
Pvalue_mlog 6.698970004
p-Value (text) (S-DCT concentration)
OR or beta NR
95% CI (text) NR
Platform [SNPs passing QC] Illumina [7,537,437] (Imputed)
CNV N

 

 

上面这个文件是由tab键分割的，每一列的意义如下！

Note: The SNP data in the catalog has been mapped to dbSNP Build 142 and Genome Assembly,

GRCh38/hg37.p13.

DATE ADDED TO CATALOG: Date added to catalog

PUBMEDID: PubMed identification number

FIRST AUTHOR: Last name of first author

DATE: Publication date (online (epub) date if available)

JOURNAL: Abbreviated journal name

LINK: PubMed URL

STUDY: Title of paper (linked to PubMed abstract)

DISEASE/TRAIT: Disease or trait examined in study

INITIAL SAMPLE SIZE: Sample size for Stage 1 of GWAS

REPLICATION SAMPLE SIZE: Sample size for subsequent replication(s)

REGION: Cytogenetic region associated with rs number (NCBI)

CHR_ID: Chromosome number associated with rs number (NCBI)

CHR_POS: Chromosomal position associated with rs number (dbSNP Build 132,

NCBI)

REPORTED GENE (S): Gene(s) reported by author

MAPPED GENE(S): Gene(s) mapped to the strongest SNP (NCBI). If the SNP is

located within a gene, that gene is listed. If the SNP is intergenic, the upstream and

downstream genes are listed, separated by a hyphen. UPSTREAM_GENE_ID:

Entrez Gene ID for nearest upstream gene to rs number, if not within gene (NCBI)

DOWNSTREAM_GENE_ID: Entrez Gene ID for nearest downstream gene to rs

number, if not within gene (NCBI)

SNP_GENE_IDS: Entrez Gene ID, if rs number within gene; multiple genes

denotes overlapping transcripts (NCBI)

UPSTREAM_GENE_DISTANCE: distance in kb for nearest upstream gene to rs

number, if not within gene (NCBI)

DOWNSTREAM_GENE_DISTANCE: distance in kb for nearest downstream

gene to rs number, if not within gene (NCBI)

STRONGEST SNP-RISK ALLELE: SNP(s) most strongly associated with trait +

risk allele (? for unknown risk allele). May also refer to a haplotype.

SNPS: Strongest SNP; if a haplotype is reported above, may include more than one

rs number (multiple SNPs comprising the haplotype)

MERGED: denotes whether the SNP has been merged into a subsequent rs record

(0 = no; 1 = yes; NCBI)

SNP_ID_CURRENT: current rs number (will differ from strongest SNP when

merged = 1)

CONTEXT: SNP functional class (NCBI)

INTERGENIC: denotes whether SNP is in intergenic region (0 = no; 1 = yes;

NCBI)

RISK ALLELE FREQUENCY: Reported risk allele frequency associated with

strongest SNP

P-VALUE: Reported p-value for strongest SNP risk allele (linked to dbGaP

Association Browser)

PVALUE_MLOG: -log(p-value)

P-VALUE (TEXT): Information describing context of p-value (e.g. females,

smokers).

Note that p-values are rounded to 1 significant digit (for example, a published pvalue of 4.8 x 10-7 is rounded to 5 x 10-7).

OR or BETA: Reported odds ratio or beta-coefficient associated with strongest

SNP risk allele

95% CI (TEXT): Reported 95% confidence interval associated with strongest SNP

risk allele

PLATFORM (SNPS PASSING QC): Genotyping platform manufacturer used in

Stage 1; also includes notation of pooled DNA study design or imputation of

SNPs, where applicable

CNV: Study of copy number variation (yes/no)

Updated: January 13, 2015

 

 

 

 

五 05

Bioconductor的数据包library(biomaRt)简介

Posted on by

 

这是发布在bioconductor平台上面的一个数据库文件，可以通过R里面下载安装并使用，非常方便。其实在ensembl数据库里面也有一个biomart，我之前也讲过这个平台，非常好用，可以把任意的数据库之间的ID号进行转换。

为了更好的理解和掌握biomaRt，我们可以先通过在线资源来了解一下它的原型biomart （http://www.biomart.org）。 biomart是为生物科研提供数据服务的免费软件，它为数据下载提供打包方案。它有许多成功的应用实例，比如欧洲生物信息学中心（The European Bioinformatics Institute ，EBI）维护的Ensembl数据库（http://www.ensembl.org/）就使用biomart提供数据批量下载服务， 还有COSMIC, Uniprot, HGNC, Gramene, Wormbase以及dbSNP等。

这个就是一个R平台的biomart而已，但是非常好用！

> library(biomaRt)

> head(listMarts(), 3)

biomart                           version

1    ensembl      ENSEMBL GENES 79 (SANGER UK)

2        snp  ENSEMBL VARIATION 79 (SANGER UK)

3 regulation ENSEMBL REGULATION 79 (SANGER UK)

这是这个biomart最具有代表性的三个数据库，用listMarts()可以查看得知，它总共有58个数据库。

ensembl <-  useMart("ensembl", dataset = "hsapiens_gene_ensembl")

这是创建了人的ensembl数据库对象

> head(listFilters(ensembl), 3)

name     description

1 chromosome_name Chromosome name

2           start Gene Start (bp)

3             end   Gene End (bp)

可以看到对人的数据库ensembl来说，有多种字段可以来挑选自己感兴趣的东西，最常用的的当然是染色体号及起始终止坐标啦，用listFilters(ensembl),以查看得知，它总共有284中挑选感兴趣数据的方式。

既然 chromosome_name是其中一个挑选字段，那么我们就可以看看，是如何进行挑选的

用filterOptions(myFilter, ensembl)可以看到它挑选参数非常之多，远不止我们所认为的染色体号码。

染色体号一般是1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y,MT

还有一堆稀奇古怪的标志，LRG_101,LRG_102,LRG_103,LRG_104，因为我们组装好的人的标准基因组还有很多小的片段不被计入染色体中。

然后还可以看到人的ensembl数据库对象，有很多的属性，最常见的当然是基因ID和转录本ID和蛋白的ID号啦！

> head(listAttributes(ensembl), 3)

name           description

1       ensembl_gene_id       Ensembl Gene ID

2 ensembl_transcript_id Ensembl Transcript ID

3    ensembl_peptide_id    Ensembl Protein ID

用listAttributes(ensembl),,以查看得知，它总共有1166个ID号，太恐怖了，我实在是没有想到！

 

那么接下来我简单讲讲这个包的几个应用吧

首先是根据entrez ID号来找

ensembl <-  useMart("ensembl", dataset = "hsapiens_gene_ensembl")

这样就得到了人的信息，然后我探究以下两个基因的其它信息。

entrzID=c("672","1")

getBM(attributes=c("entrezgene","hgnc_symbol","ensembl_gene_id"), filters = "entrezgene", values =entrzID, mart=ensembl )

entrezgene hgnc_symbol ensembl_gene_id

1          1        A1BG ENSG00000121410

2        672       BRCA1 ENSG00000012048

3        672       BRCA1         LRG_292

 

其实这个函数很简单，就是根据自定义的entrzID这个变量来找到一些数据，数据的属性是我自己定义的entrezgene","hgnc_symbol","ensembl_gene_id"，所以它就显示这个信息给我，在我之前弄好的人的数据库里面寻找！listAttributes(ensembl),,以查看得知，它总共有1166个ID号，就是说，你可以挑选你想要的基因的1166种信息，包罗万象！！！

其它功能也是很简单的啦，自己多看帮助文档！

 

从上面的操作来看，使用biomaRt只需要两步，1，指定mart数据库，2，使用getBM获得注释。但是首先，我们如何知道有哪些服务器，以及这些服务器上哪些数据库呢？其次，我们如何获阳getBM中attributes，filters的正确设置呢？

关于第一个问题，我们可以使用biomaRt中的listMarts以及listDatasets两个函数来解决。

> marts <- listMarts(); head(marts) #查看当前可用的数据源 ，总共有58个数据源。

> ensembl <- useMart("ensembl") #使用ensembl数据源

> datasets <- listDatasets(ensembl); datasets[1:10,] #查看ensembl中可用数据库，共有69个物种的数据库！

对于第二个问题，我们使用biomaRt中的listFilters以及listAttributes两个函数来解决。

> mart <- useMart("ensembl", "hsapiens_gene_ensembl")  #首先使用人的数据库

>listAttributes(ensembl) #,以查看得知，它总共有1166个ID号，就是说，人的数据库可供挑选的信息多达1166种。

> filters <- listFilters(mart); filters[grepl("entrez", filters[,1]),] #总共有284中挑选感兴趣数据的方式。

最后的问题是，biomaRt会被如何使用呢？我们做注释的时候，怎么就想到要使用biomaRt呢？因为在注释上，各种ID，symbol, name之间的转换都可以考虑使用biomaRt来做。更重要的是，biomaRt还会有很多SNP, alternative splicing, exon, intron, 5’utr, 3’utr等等信息。当然，只要能做也数据库并使用SQL访问的数据都可以使用biomaRt来获取。所以我们的思路可以更加发散一些。