26

一个基因的生信之旅

感觉大家对很多生物信息学的术语都不甚了解,我这里简单的从一个基因开始,扩展开来讲一讲生信数据库,及它相关的一些术语!

我要讲的基因是BRCA1,这是一个与乳腺癌以及卵巢癌都息息相关的基因。而BRCA1是它的英文缩写简称,也是通常学者们进行交流十它的名字。它的全称是breast cancer 1,每个基因都会有一个简称,比如下面这些,在human里面这些简称多大47732个,正常人都不会认识它们所有,只需要碰到了去数据库搜索即可,但是搞医疗健康的,必须熟悉癌症50基因。

一个基因的生信之旅247

这样的缩写简称其实弊端很多,单词毕竟是有限的,而且缩写也没有语义。所以NCBI给每个基因都定义了一个entrez ID号,是整数的排序,具体大家可以去看NCBI发的一篇文献,专门讲解了entrez ID号的好处。

1 A1BG

2 A2M

3 A2MP1

9 NAT1

10 NAT2

11 NATP

12 SERPINA3

13 AADAC

14 AAMP

这里我们来找一下我们的BRCA1这个基因在生物信息数据库里面的其它信息,在NCBI的ftp里面有一个文件是Homo_sapiens.gene_info里面包含着人类所以基因的全部信息

9606  首先这个基因在human上面的,而human被NCBI定义的taxid是9606

672  然后这个基因的被NCBI定义的entrez ID号是672

BRCA1  这个当然就是这个基因的英文缩写名称啦

-      这个表明这个基因在负链什么

BRCAI|BRCC1|BROVCA1|FANCS|IRIS|PNCA4|PPP1R53|PSCP|RNF53

这个可能是基因以前的名称,或者是在其它研究领域的一些名称。MIM:113705|HGNC:HGNC:1100|Ensembl:ENSG00000012048|HPRD:00218|Vega:OTTHUMG00000157426

这里面包含在它在其它数据库的信息,我们的NCBI用entrez ID号672来标识它,相应的ensembl数据用ensembl ID号ENSG00000012048来标识它,还有什么MIM数据库,HGNC数据库,Vega数据库我就不详细讲啦

17 17q21 这个说明它在human的17号染色体的位置信息

下面一堆都是这个基因的描述,它的功能等等。

breast cancer 1, early onset protein-coding BRCA1 breast cancer 1, early onset

O BRCA1/BRCA2-containing complex, subunit 1|Fanconi anemia, complementation group S|RING finger protein 53|breast and ovarian cancer susceptibility protein 1|breast and ovarian cancer sususceptibility protein 1|breast cancer type 1 susceptibility protein|protein phosphatase 1, regulatory subunit 53

20150201

这样我们就把好几个数据库给串起来了,也大致了解了一个基因的各种信息,但是,这样肯定是不够的。

接下来我们就不用BRCA1来称呼这个基因了,我们统一用NCBI定义entrez ID号672来称呼这个基因,当然用ensembl ID号ENSG00000012048也可以,它们都是比较通用的。

ENSG00000012048 672 这个基因在GO数据库里面可以找到67个功能信息,分别是以下

GO:0000151 GO:0000724 GO:0000724 GO:0000794 GO:0003677 GO:0003684 GO:0003713 GO:0003723 GO:0004842 GO:0005515 GO:0005634 GO:0005654 GO:0005694 GO:0005737 GO:0005886 GO:0006260 GO:0006281 GO:0006301 GO:0006302 GO:0006302 GO:0006349 GO:0006357 GO:0006359 GO:0006633 GO:0006915 GO:0006974 GO:0006978 GO:0007059 GO:0007098 GO:0008270 GO:0008274 GO:0008630 GO:0009048 GO:0010212 GO:0010575 GO:0010628 GO:0015631 GO:0016567 GO:0016874 GO:0019899 GO:0030521 GO:0030529 GO:0031398 GO:0031436 GO:0031572 GO:0031625 GO:0035066 GO:0035067 GO:0042127 GO:0042981 GO:0043009 GO:0043234 GO:0043627 GO:0044030 GO:0044212 GO:0045717 GO:0045739 GO:0045766 GO:0045892 GO:0045893 GO:0045893 GO:0045944 GO:0045944 GO:0046600 GO:0050681 GO:0051571 GO:0051572 GO:0051573 GO:0051574 GO:0051865 GO:0070512 GO:0070531 GO:0071158 GO:0071356 GO:0071681 GO:0085020 GO:1902042 GO:2000378 GO:2000617 GO:2000620

由于GO太多了,我简单讲几个

ubiquitin ligase complex

double-strand break repair via homologous recombination

double-strand break repair via homologous recombination

condensed nuclear chromosome

DNA binding

damaged DNA binding

transcription coactivator activity

RNA binding

ubiquitin-protein transferase activity

protein binding

都是描述这个基因的功能的。

到这里我们大致了解了这个基因的功能,但是还不够。

然后可以查到它有一下6个转录本,都有二十多个外显子。

NR_027676

NM_007300

NM_007299

NM_007298

NM_007297

NM_007294

在hg19这个参考基因组的起始终止坐标,还有各个外显子的起始终止坐标都能找到。

41196311,41199659,41201137,41203079,41209068,41215349,41215890,41219624,41222944,41226347,41228504,41234420,41242960,41243451,41247862,41249260,41251791,41256138,41256884,41258494,41267742,41276033,41277198

 

41197819,41199720,41201211,41203134,41209152,41215390,41215968,41219712,41223255,41226538,41228631,41234592,41243049,41246877,41247939,41249306,41251894,41256278,41256973,41258550,41267796,41276132,41277340

http://asia.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000012048;r=17:43044295-43125483

在ensembl里面关于这个基因的描述如下。

breast cancer 1, early onset [Source:HGNC Symbol;Acc:HGNC:1100]

BRCC1, FANCS, PPP1R53, RNF53

Chromosome 17: 43,044,295-43,125,483 reverse strand.

chromosome:GRCh38:CM000679.2:43044295:43125483:1

This gene has 29 transcripts (splice variants), 63 orthologues, is a member of 4 Ensembl protein families and is associated with 11 phenotypes.

RefSeq Gene ID 672

Uniprot identifiers: P38398

而且ensembl里面可以可视化这个基因的所有信息。

然后简单检索一下关于这个BRCA1基因的文献发表状况,居然多达2111篇文献,看来这个基因很火呀!!!

awk '{if ($1==9606 && $2==672) print }' gene2pubmed |wc

9606 672 1676470

9606 672 2001833

9606 672 2270482

9606 672 4506230

9606 672 7481765

9606 672 7545954

9606 672 7550349

9606 672 7795652

9606 672 7894491

9606 672 7894492

第三列1676470等编号是pubmed数据库的文献编号,可以直接找到关于这个基因的文献发表情况。

而直接在NCBI的pubmed数据库里面可以搜到多达11339篇文献。

esearch -db pubmed -query 'BRCA1'

Esearch这个程序是NCBI提供的,挺好用的,希望大家可以熟悉一下。

esearch -db pubmed -query 'BRCA1' | efetch -format docsum |   xtract -pattern DocumentSummary -present Author -and Title     -element Id -first "Author/Name" -element Title  >BRCA1.pubmed

一个基因的生信之旅4634

用这个代码,可以找到所有关于这个BRCA1基因的文献的作者及标题,这样可以统计在这个基因领域的研究者最出名的是谁。

至于这个基因的序列,及其转录本翻译的蛋白我就不列了,太长了,而且占位子

 

 

 

 

25

人为创造几个测序数据然后用soap组装成基因组

这里我选取酵母基因组来组装,以为它只有一条染色体,而且本身也不大!

人为创造几个测序数据然后用soap组装成基因组130

这个文件就4.5M,然后第一行就是序列名,第二列就是序列的碱基组成。共4641652个碱基。

我写一个perl程序来人为的创造一个测序文件

人为创造几个测序数据然后用soap组装成基因组58

这样我们的4.5M基因组就模拟出来了486M的单端100bp的测序数据,而且是无缝连接,按照道理应该很容易就拼接的。

/home/jmzeng/bio-soft/SOAPdenovo2-bin-LINUX-generic-r240/SOAPdenovo-127mer

all -s config_file -K 63 -R -o graph_prefix 1>ass.log 2>ass.err

人为创造几个测序数据然后用soap组装成基因组331

可以看到组装效果还不错哦,然后我模拟了一个测试数据,再进行组装一次,这次更好!

其实还可以模拟双端测序,应该就能达到百分百组装了。

但是由于我代码里面选取的是80在随机错开,所以我把kmer的长度设置成了81来试试看,希望这样可以把它完全组装成一条e-coli基因组。

/home/jmzeng/bio-soft/SOAPdenovo2-bin-LINUX-generic-r240/SOAPdenovo-127mer

all -s config_file -K 81 -R -o graph_prefix 1>ass.log 2>ass.err

但是也没有什么实质性的提高,虽然理论上是肯定可以组装到一起!

那我再模拟一个双端测序吧,中间间隔200bp的。

 

25

基因组组装软件SOAPdenovo安装使用

一.下载并安装这个软件

下载地址进下面,但是下载源码安装总是很困难,我直接下载bin文件可执行程序。

基因组组装软件SOAPdenovo安装使用104

解压进入目录

首先make

然后make install即可

基因组组装软件SOAPdenovo安装使用731

安装总是失败,我也不知道怎么回事,懒得解决了。

直接去我老师那里把这个程序拷贝进来了。

https://github.com/aquaskyline/SOAPdenovo2/archive/master.zip

http://sourceforge.net/projects/soapdenovo2/files/SOAPdenovo2/bin/r240/SOAPdenovo2-bin-LINUX-generic-r240.tgz/download

http://sourceforge.net/projects/soapdenovo2/files/latest/download?source=files

也可以直接下载bin程序

基因组组装软件SOAPdenovo安装使用1035

二.准备测试数据

基因组组装软件SOAPdenovo安装使用742

类似于这样的几个文库的左右两端测序数据。

我这里用一个小样本的单端数据做测试

基因组组装软件SOAPdenovo安装使用783

三,参考命令

You may run it like this:

参考:http://www.plob.org/2012/07/06/2537.html

https://github.com/aquaskyline/SOAPdenovo2

总共就四个步骤,介绍如下。

 

./pregraph_sparse [parameters]
./SOAPdenovo-63mer contig [parameters]
./SOAPdenovo-63mer map [parameters]
./SOAPdenovo-63mer scaff [parameters]

 

i) preparing the pregraph. This step is similar to velveth for velvet.
ii) Determining contigs. This step is similar to velvetg for velvet.
iii) Mapping back reads on to contigs.
iv) Assembling contigs into scaffolds.

 

 SOAPdenovo-63mer  sparse_pregraph  -s config_file -K 45 -p 28 -z 1100000000 -o outPG
 SOAPdenovo-63mer contig  -g outPG
 SOAPdenovo-63mer map  -s config_file -g outPG -p 28
 SOAPdenovo-63mer  scaff   -g outPG -p 28

基因组组装软件SOAPdenovo安装使用1629

官网给出的步骤如下

基因组组装软件SOAPdenovo安装使用1641

这个命令还需要一个配置文件

max_rd_len=99 设置最大reads长度,具体情况具体定义

[LIB] 第一个文库数据

avg_ins=225

reverse_seq=0

asm_flags=3

rank=1

q1=runPE_1.fq

q2=runPE_2.fq

[LIB] 第二个文库数据

avg_ins=2000

reverse_seq=1

asm_flags=2

rank=2

q1=runMP_1.fq

q2=runMP_2.fq

也可以全部一次性的搞一个命令

all -s config_file -K 63 -R -o graph_prefix 1>ass.log 2>ass.err

我简单修改了一下参考博客的代码跟官网的代码,然后运行了我自己的代码

/home/jmzeng/bio-soft/SOAPdenovo2-bin-LINUX-generic-r240/SOAPdenovo-127mer

all -s config_file -K 63 -R -o graph_prefix 1>ass.log 2>ass.err

反正我也不懂,就先跑跑看咯

我选取的是7个单端数据,所以我的配置文件是

max_rd_len=500

[LIB]

avg_ins=225

reverse_seq=0

asm_flags=3

rank=1

p=SRR072005.fa

p=SRR072010.fa

p=SRR072011.fa

p=SRR072012.fa

p=SRR072013.fa

p=SRR072014.fa

p=SRR072029.fa

四.输出数据解读

好像我的数据都比较小,就7个三百多兆的fasta序列,几个小时就跑完啦

四个步骤都有输出数据

基因组组装软件SOAPdenovo安装使用2446

好像组装效果惨不忍睹呀!共86万的contig,50多万的scaffold

scaffolds>100  505473 99.60%

scaffolds>500  113523 22.37%

scaffolds>1K   48283 9.51%

scaffolds>10K  0 0.00%

scaffolds>100K 0 0.00%

scaffolds>1M   0 0.00%

这其实都相当于没有组装了,因为我的测序判断本来就很多是大于500的!

可能是我的kmer值选取的不对

Kmer为63跑出来的效果不怎么好,86万的contig,50万的scaffold的

Kmer为35跑出来的效果更惨,203万的contig,近60万的scaffold。

我觉得问题可能不是这里了,可能是没有用到那个20k和3k的双端测序库,唉,其实我习惯了illumina的测序数据,不太喜欢这个454的

感觉组装好难呀,业余时间搞不定呀,希望有高手能一起交流,哈哈,我自己再慢慢来试试。

 

 

 

 

 

 

 

 

24

草莓基因组数据预处理

今天先 对7个单端数据做处理,是454数据,平均长度300bp左右,明天再处理3KB和20KB的配对reads。

首先跑fastqc

打开一个个看结果

草莓基因组数据预处理28

可以看到前面一些碱基的质量还是不错的, 因为这是454平台测序数据,序列片段长度差异很大,一般前四百个bp的碱基质量还是不错的,太长了的测序片段也不可靠

草莓基因组数据预处理39

重点在下面这个图片,可以看到,前面的4个碱基是adaptor,肯定是要去除的,不是我们的测序数据。是TCAG,需要去除掉。

草莓基因组数据预处理118

所以我们用了 solexaQA 这个套装软件对原始测序数据进行过滤

草莓基因组数据预处理214

可以看到过滤的非常明显!!!甚至有个样本基本全军覆没了!然后我查看了我的批处理脚本,发现可能是perl DynamicTrim.pl -454 $id这个参数有问题

for id in *fastq

do

echo $id

perl DynamicTrim.pl -454 $id

done

for id in *trimmed

do

echo $id

perl LengthSort.pl $id

done

 

可以看到末尾的质量差的碱基都被去掉了,但是头部的TCAG还是没有去掉。

草莓基因组数据预处理425

处理完毕后的数据如下:

草莓基因组数据预处理475

 

 

 

 

 

 

24

solexaQA 对测序数据进行简单过滤

一.下载该软件

http://solexaqa.sourceforge.net/index.htm

下载解压开

现在已经把它的三个功能整合到一起啦

之前是分开的程序,我主要用它的两个perl 程序,我比较喜欢之前的版本,所以下面的讲解也是基于这两个perl程序。

这两

solexaQA 对测序数据进行简单过滤295

个主要是对reads进行最大子串的截取

solexaQA 对测序数据进行简单过滤476

 

二.准备数据。

就是我们测序得到的原始数据。

第一个就是质量控制,一般是以20为标准,当然你也可以自己设定,该软件质控的原理如下:

使用默认的参数值(defaults to P = 0.05, or equivalently, Q = 13)

solexaQA 对测序数据进行简单过滤846

 

基本上就是取符合阈值的最大子串。

二:命令使用很简单一般使用DynamicTrim与LengthSort.pl就可以了

for id in *fastq

do

echo $id

perl DynamicTrim.pl -454 $id

done

for id in *trimmed

do

echo $id

perl LengthSort.pl $id

done

首先使用DynamicTrim.pl程序,非常耗时间

几个小时完毕之后

solexaQA 对测序数据进行简单过滤2335

查看,产出文件如下

solexaQA 对测序数据进行简单过滤2497

 

可以看到丢弃的不多,也就三五百M的

简单查看丢弃的,都是短的。

perl -lne '{print length if $.%4==2}' SRR1793918.fastq.trimmed.discard |head

solexaQA 对测序数据进行简单过滤2678

用这个脚本查看,可知好像都是短于25个碱基的被舍弃掉了,这个参数可以调整的。

接下来就可以用这些数据进行数据分析了

 

24

旧版本blast详解

其实我现在一般都用的是blast++了,也专门写了篇日志介绍它!

但是看到一些就的服务器上面只有blast,所以就搜了一些它的用法。

主要参考 http://www.bio.ku.dk/nuf/resources/BLAST_index.htm

很简单的两个步骤

首先建库formatdb -i Cad16_aa.fasta -p T -o F

就是把 Cad16_aa.fasta这个序列文件变成blast专用的库,-p选项中的T是代表蛋白库

然后就比对咯,比对程序有六个,需要用-p来选择

blastall -p blastx -d nr -i 19A.fa -o 19A.outm -v 1 -b 1 -m 8

上面这个命令就是选择了blastx这个比对程序,数据库是nr ,输入的查询序列是 19A.fa

然后我们输出格式的m8,这个格式很重要,我们还可以设置-a控制cpu数量,和-e控制阈值

BLAST programs

blastp Protein query > Protein database
blastn Nucleotide query > Nucleotide database
blastx Nucleotide query > Protein database (via translated query)
tblastn Protein query > Nucleotide database (via translated database)
tblastx Nucleotide query > Nucleotide database (via translated query and database) 

 

Formatting database for local BLAST

- Show a list of all arguments.
-i Input file(s) for formatting. Optional.
-p Type of file [T/F]. T = protein, F = nucleotide. Default = T.
-o Parse option [T/F]. T = Parse SeqId and create indexes, F = Do not parse or create indexes.
23

与hg19的突变相关的一些数据解释。

http://statgenpro.psychiatry.hku.hk/limx/kggseq/download/resources/

这个网站收集了大部分资料,我们就用它的,如果它倒闭了,大家再想办法去搜索吧。

与hg19的突变相关的一些数据解释210

其实这些文件都是基于NCBI以及UCSC和ensembl数据库的文件用一些脚本转换而来的,都是非常简单的脚本。

首先我们看看humandb/hg19_refGene.txt 这个文件,总共2.5万多个基因的共5万多个转录本。

与hg19的突变相关的一些数据解释325

19     可能是entrez ID,但是又不像。

NM_001291929    参考基因名

chr11   染色体

-

89057521

89223909

89059923

89223852

17      89057521,89069012,89070614,89073230,89075241,89088129,89106599,89133184,89133382,89135493,89155069,89165951,89173855,89177302,89182607,89184952,89223774,       89060044,89069113,89070683,89073339,89075361,89088211,89106660,89133247,89133547,89135710,89155150,89166024,89173883,89177400,89182692,89185063,89223909,

0

NOX4    基因的英文简称,通俗名

cmpl

cmpl

2,0,0,2,2,1,0,0,0,2,2,1,0,1,0,0,0,

然后我们看看hg19_snp141.txt这个文件

1       10229   A       -       .

1       10229   AACCCCTAACCCTAACCCTAAACCCTA     -       .

1       10231   C       A       .

1       10231   C       -       .

1       10234   C       T       .

1       10248   A       T       .

1       10250   A       C       .

1       10250   AC      -       .

1       10255   A       -       .

1       10257   A       C       .

1       10259   C       A       .

1       10291   C       T       .

1       10327   T       C       .

1       10329   ACCCCTAACCCTAACCCTAACCCT        -       .

1       10330   C       -       .

1       10390   C       -       .

1       10440   C       A       .

1       10440   C       -       .

1       10469   C       G       .

1       10492   C       T       .

1       10493   C       A       .

1       10519   G       C       .

1       10583   G       A       0.144169

1       10603   G       A       .

1       10611   C       G       0.0188246

1       10617   CGCCGTTGCAAAGGCGCGCCG   -

里面记录了以hg19为参考的所有的snp位点。

 

585

ENST00000518655 基因的ensembl ID号

chr1 + 11873 14409 14409 14409

4 基因有四个外显子

11873,12594,13402,13660, 12227,12721,13655,14409, 在基因的四个外显子的坐标

0

DDX11L1 基因的通俗英文名

none none -1,-1,-1,-1,

CTTGCCGTCAGCCTTTTCTTT·····gene的核苷酸序列

23

用annovar对snp进行注释

一、下载及安装软件

这个软件需要edu邮箱注册才能下载,可能是仅对科研高校开放吧。所以软件地址我就不列了。

用annovar对snp进行注释71

它其实是几个perl程序,比较重要的是这个人类的数据库,snp注释必须的。

用annovar对snp进行注释111

参考:http://annovar.readthedocs.org/en/latest/misc/accessory/

二,准备数据

既然是注释,那当然要有数据库啦!数据库倒是有下载地址

http://www.openbioinformatics.org/annovar/download/hg19_ALL.sites.2010_11.txt.gz

也可以用命令来下载

Perl ./annotate_variation.pl -downdb -buildver hg19 -webfrom annovar refGene humandb/

然后我们是对snp-calling流程跑出来的VCF文件进行注释,所以必须要有自己的VCF文件,VCF格式详解见本博客另一篇文章,或者搜索也行

http://vcftools.sourceforge.net/man_latest.html

三、运行的命令

首先把vcf格式文件,转换成空格分隔格式文件,自己写脚本也很好弄

perl convert2annovar.pl  -format vcf

/home/jmzeng/raw-reads/whole-exon/snp-calling/tmp1.vcf >annovar.input

用annovar对snp进行注释886

变成了空格分隔的文件

用annovar对snp进行注释900

然后把转换好的数据进行注释即可

./annotate_variation.pl -out ex1 -build hg19 example/ex1.avinput humandb/

 

四,输出文件解读

用annovar对snp进行注释1002

23

Snp-calling流程(BWA+SAMTOOLS+BCFTOOLS)

比对可以选择BWA或者bowtie,测序数据可以是单端也可以是双端,我这里简单讲一个,但是脚本都列出来了。而且我选择的是bowtie比对,然后单端数据。

首先进入hg19的目录,对它进行两个索引

samtools faidx hg19.fa

Bowtie2-build hg19.fa hg19

我这里随便从26G的测序数据里面选取了前1000行做了一个tmp.fa文件,进入tmp.fa这个文件的目录进行操作

Bowtie的使用方法详解见http://www.bio-info-trainee.com/?p=398

bowtie2 -x ../../../ref-database/hg19 -U  tmp1.fa -S tmp1.sam

samtools view -bS tmp1.sam > tmp1.bam

samtools sort tmp1.bam tmp1.sorted

samtools index tmp1.sorted.bam 

samtools mpileup -d 1000  -gSDf   ../../../ref-database/hg19.fa  tmp1.sorted.bam |bcftools view -cvNg -  >tmp1.vcf

然后就能看到我们产生的vcf变异格式文件啦!

 

当然,我们可能还需要对VCF文件进行再注释!

要看懂以上流程及命令,需要掌握BWA,bowtie,samtools,bcftools,

数据格式fasta,fastq,sam,vcf,pileup

 

如果是bwa把参考基因组索引化,然后aln得到后缀树,然后sampe对双端数据进行比对

首先bwa index 然后选择算法,进行索引。

然后aln脚本批量处理

==> bwa_aln.sh <==

while read id

do

echo $id

bwa aln hg19.fa $id >$id.sai

done <$1

然后sampe脚本批量处理

==> bwa_sampe.sh <==

while read id

do

echo $id

bwa sampe hg19.fa $id*sai $id*single >$id.sam

done <$1

然后是samtools的脚本

==> samtools.sh <==

while read id

do

echo $id

samtools view -bS $id.sam > $id.bam

samtools sort $id.bam $id.sorted

samtools index $id.sorted.bam

done <$1

然后是bcftools的脚本

==> bcftools.sh <==

while read id

do

echo $id

samtools mpileup -d 1000  -gSDf  ref.fa $id*sorted.bam |bcftools view -cvNg -  >$id.vcf

done <$1

 

 

==> mpileup.sh <==

while read id

do

echo $id

samtools mpileup -d 100000 -f hg19.fa $id*sorted.bam >$id.mpileup

done <$1

 

22

WordPress博客统计文章阅读次数及访客数并刷访问数

需要插件和自己修改主题下面的foot.php代码。

参考 http://jingyan.baidu.com/article/ae97a646ce37c2bbfd461d01.html

步骤如下:

1、登陆到wp后台,鼠标移动到左侧菜单的“插件”链接上,会弹出子菜单,点击子菜单的“安装插件”链接

2、WP-PostViews插件显示wordpress文章点击浏览量

在“安装插件”链接页面的搜索框中输入“WP-PostViews”,然后回车

3、WP-PostViews插件显示wordpress文章点击浏览量

在搜索结果页面点击“WP-PostViews”插件内容区域的“现在安装”按钮

4、WP-PostViews插件显示wordpress文章点击浏览量

程序自动下载插件到服务器并解压安装,一直等到安装成功信息出现,然后在安装成功提示页面点击“启动插件”链接。

5、WP-PostViews插件显示wordpress文章点击浏览量

页面会自动跳转到“已安装插件”页面,在已安装插件列表中我们可以看到“Form Manager”插件已经处于启用状态(插件名下是“停用”链接)。

Wordpress博客统计文章阅读次数及访客数并刷访问数599

有了这个插件之后,我们的整个网页环境里面就多了一个 the_views()函数,它统计着每个文章的点击数,这样我们之前的网页就能显示点击数了。

Wordpress博客统计文章阅读次数及访客数并刷访问数673

这个是我现在用的主题的php代码,把文章用span标记隔开了,而且显示着上面php代码里面的每一个内容包括日期,分类,标签,评论等等

Wordpress博客统计文章阅读次数及访客数并刷访问数742

其中thez-view()这个函数返回的不仅仅是一个访客数,但是我的文章的访客都太少了,所以我写了一个脚本帮我刷一刷流量。

[perl]

use List::MoreUtils qw(uniq);

$page='http://www.bio-info-trainee.com/?paged=';

foreach (1..5){   #我的文章比较少,就42个,所以只有5个页面

$url_page=$page.$_;

$tmp=`curl $url_page`;

#@p=$tmp=~/p=(\d+)/;

$tmp =~ s/(p=\d+)/push @p, $1/eg; #寻找p=数字这样的标签组合成新的网页地址

}

@p=uniq @p;

print "$_\n" foreach @p; #可以找到所有42个网页的地址

foreach (@p){

$new_url='http://www.bio-info-trainee.com/?'.$_;

`curl $new_url` foreach (1..100); #每个网页刷一百次

}

[/perl]

大家可以看到这个网页被刷的过程,从15到21到27直到100
Wordpress博客统计文章阅读次数及访客数并刷访问数1231

大家现在再去看我的网页,就每个文章都有一百的访问量啦!

http://www.bio-info-trainee.com/

 

21

Hg19基因组的分析

下载地址我就不贴了,随便谷歌一下即可!

Genome Reference Consortium Human  ---》  GRCh3

Feb. 2009 (hg19, GRCh37)这个是重点

Mar 2006 assembly = hg18 = NCBI36.

May 2004 assembly = hg17 = NCBI35.

July 2003 assembly = hg16 = NCBI34

以前的老版本就不用看啦,现在其实都已经有hg38出来啦,GRCh38 (NCBI) and hg38(UCSC)

参考:http://age.wang.blog.163.com/blog/static/119252448201092284725460/

http://www.ncbi.nlm.nih.gov/projects/genome/assembly/grc/human/

Hg19基因组的分析570

人的hg19基因组是3G的大小,因为一个英文字符是一个字节,所以也是30亿bp的碱基。

包括22条常染色体和X,Y性染色体及M线粒体染色体。

Hg19基因组的分析643

查看该文件可以看到,里面有很多的N,这是基因组里面未知的序列,用N占位,但是觉得部分都是A.T.C.G这样的字符,大小写都有,分别代表不同的意思。

然后我用linux的命令统计了一下里面这个文件的行数,

perl -lne 'END { print $. }'  hg19.fa

awk 'END { print NR }'  hg19.fa

wc -l hg19.fa

Hg19基因组的分析834

然后我写了一个脚本统计每条染色体的长度,42秒钟完成任务!

Hg19基因组的分析1125

看来这个服务器的性能还是蛮强大的,读取文件非常快!

[perl]

while(<>){

        chomp;

        if  (/>/){

if  (exists $hash_chr{$key} ){

$len = length $hash_chr{$key};

print "$key   =>   $len\n";

}

undef %hash_chr;

$key=$_;

}

else {

$hash_chr{$key}.=$_;

}

}

[/perl]

 

然后我用seed统计了一下hg19的词频(我不知道生物信息学里面的专业描述词语是什么)

Hg19基因组的分析1171

我的程序耗费了42分钟才跑完,感觉我写的程序应该是没有问题的,让我吃惊的是总共竟然只有105万条独特的10bp短序列。然后我算了一下4的10次方,(⊙o⊙)…悲剧,原来只有1048576,之所以出现这种情况,是因为里面有N这个字符串,不仅仅是A.T.C.G四个字符。我用grep -v N seed10.txt |wc -l命令再次统计了一下,发现居然就是1048576,也就是说,任意A.T.C.G四个字符组成的10bp字符串短序列在人的基因组里面都可以找到!!!

Hg19基因组的分析1407

然后我测试了一下,还是真是这样的,真是一个蛮有意思的现象。虽然我无法解释为什么,但是根据这个结果我们可以得知连续的A或者T在人类基因组里面高频出现,而连续的G或者C却很少!

如果我们储存这个10bp字符串的同时,也储存着它们在基因组的位置,那么就可以根据这个seed来进行比对,这就是blast的原理之一!

 

 

21

积累的一些perl代码分享

以前的一下perl代码分享

今天去参加了开源中国的一个源创会,感觉好隆重的样子,近五百人,BAT的工程师都过来演讲了,可都是数据库相关的, 我一个的都没有听懂,但是茶歇的披萨我倒是吃了不少。

说到开源中国,我想起来了我以前在上面分享的代码,上去看了看,竟然有那么多的访问量了,让我蛮意外的,那些代码完全是我学习perl的历程的真实写照。

http://www.oschina.net/code/list_by_user?id=1990747

Continue reading

21

Linux服务器基础知识

想了想,既然是菜鸟教程,那就索性再介绍点更基础的东西,基本上只要是大学毕业的都能看懂,不需要懂计算机了。首先讲讲linux服务器吧,因为生物信息也算是半个大数据分析,所以我们平常的办公电脑一般都是不能满足需求的,大部分实验室及公司都会自己配置好服务器给菜鸟们用,菜鸟们首先要拿到服务器的IP和高手给你的用户名和密码。

一般我们讲服务器,大多是linux系统,而我这里所讲的linux系统呢,特指ubuntu,其余的我懒得管了,大家也不要耗费无谓的时间纠结那些名词的不同!

登录到服务器有两种方法,一种是ssh,传输你的命令给服务器执行,另一种是ftp,和服务器交换文件。而ssh我们通常用putty,xshell等等。ftp呢,我们可以用winscp,xshell,所以我一直都用xshell,因为它两者都能搞定!

Xshell软件自行搜索下载,打开之后新建一个连接,然后登陆即可。

Linux服务器基础知识405

然后输入以下命令,可以查看服务器配置,包括cpu。内存,还有硬盘

cat /proc/cpuinfo |grep pro|wc -l

free -g

df -h

Linux服务器基础知识488

 

这个服务器配置好一点,有80个cpu,内存256G,硬盘有2个11T的,是比较成熟的配置。

Linux服务器基础知识536

 

这个是一个小型服务器。也就24个核,64G的内存,但是存储量有点小呀,其实可以随便花几百块钱买个1T的硬盘挂载上去的。

然后linux的其它命令大家就得自己去搜索一个个使用,然后熟悉,记牢,然后创新啦!

我随便敲几个我常用的吧: ls cd mkdir rm cp cat head tail more less diff grep awk sed grep perl 等等!

呀,突然间发现我才介绍了ssh的方法登陆服务器并且发送命令在服务器上面运行,下面贴图如何传输文件。一般xshell的菜单里面有绿的文件夹形式的标签就是打开ftp文件传输,这种可视化的软件,大家慢慢摸索吧!

Linux服务器基础知识830

 

 

 

 

20

自己动手写bowtie第一讲:BWT算法详解并建立索引

首先,什么是BWT,可以参考博客

http://www.cnblogs.com/xudong-bupt/p/3763814.html

他讲的非常好。

一个长度为n的串A1A2A3...An经过旋转可以得到

A1A2A3...An

A2A3...AnA1

A3...AnA1A2

...

AnA1A2A3...

n个串,每个字符串的长度都是n。

对这些字符串进行排序,这样它们之前的顺序就被打乱了,打乱的那个顺序就是index,需要输出。

首先我们测试一个简单的字符串acaacg$,总共六个字符,加上一个$符号,下次再讲$符号的意义。

BWT算法详解之一建立索引348

实现以上功能是比较简单的,代码如下

BWT算法详解之一建立索引563

但是这是对于6个字符串等小片段字符串,如果是是几千万个字符的字符串,这样转换就会输出千万的平方个字符串组成的正方形数组,是很恐怖的数据量。所以在转换的同时就不能把整个千万字符储存在内存里面。

在生物学领域,是这样的,这千万个 千万个碱基的方阵,我们取每个字符串的前20个字符串就足以对它们进行排序,当然这只是近视的,我后面会讲精确排序,而且绕过内存的方法。

Perl程序如下

[perl]

while (<>){

next if />/;

chomp;

$a.=$_;

}

$a.='$';

$len=length $a;

$i=0;

print "first we transform it !!!\n";

foreach (0..$len-1){

$up=substr($a,0,$_);

$down=substr($a,$_);

#print "$down$up\n";

#$hash{"$down$up"}=$i;

$key=substr("$down$up",0,20);

$key=$key.”\t”.substr("$down$up",$len-1);

$hash{$key}=$i;

$i++;

}

print "then we sort it\n";

foreach  (sort keys  %hash){

$first=substr($_,0,1);

$len=length;

$last=substr($_,$len-1,1);

#print "$first\t$last\t$hash{$_}\n";

print "$_\t$hash{$_}\n";

}

[/perl]

运行的结果如下

BWT算法详解之一建立索引1289

个人觉得这样排序是极好的,但是暂时还没想到如何解决不够精确的问题!!!

参考:

http://tieba.baidu.com/p/1504205984

http://www.cnblogs.com/xudong-bupt/p/3763814.html

 

20

bowtie简单使用

首先进入bowtie的主页,千万要谷歌!!!

http://bowtie-bio.sourceforge.net/bowtie2/index.shtml

image001

主页里面有下载链接,也有索引文件,当然索引文件只有人类等模式生物

下载之后是个压缩包,解压即可使用

image003

可以看到绿色的就是命令,可以添加到环境变量使用,也可以直接用全路径使用它!!!

然后example文件夹里面有所有的测试文件。

二、准备数据

我们就用软件自带的测试数据

image005

三、运行命令

分为两步,首先索引,然后比对!!!

  • 索引,bowtie2-build your-fastq-file.fa your-index-name

image008

然后你的目录就产生了六个索引文件,我给索引取名是tmp,你们可以随便取名字

  • 然后比对,分两种,一是单端测序数据,二是双端数据

重点参数的-x 和 –S ,单端是-U 双端是-1 -2

Bowtie –x tmp –U reads.fa –S hahahhha.sam

Bowtie –x tmp -1 reads1.fa -2 reads2.fa –S hahahha.sam

四:输出文件解读

就是输出了sam文件咯,这个就看我的sam文件格式讲解哈

 

 

20

perl实现二分法查找

perl实现二分法查找

在perl里面字符串跟数字是不区分的,所以写代码需要考虑到它们的区别!

首先是对数字查找来说,所有的操作符都是 < ,> ,==等等

@a=1..1000;
$b=56; #just a example
sub half_search{
 my($ref,$key,$low,$high)=@_;
 $high=@{$ref}-1 unless $high;
 if ($key < $ref->[$low] or $key > $ref->[$high]){
 print "not exists !!!\n";
 last;
 }
 if ($ref->[$low] > $ref->[$high]){
 print "not sort array !!!\n" ;
 last;
 }
 $mid=int (($low+$high)/2);
 if ($ref->[$mid] == $key) { 
 return $mid+1;
 }
 elsif ($ref->[$mid] < $key) {
 &half_search($ref,$key,$mid+1,$high);
 }
 else{
 &half_search($ref,$key,$low,$mid-1);
 }
}
print &half_search(\@a,$b);

对排序好的数字数组来说,是非常简单的,非常快速。

接下来是字符串数组的查找。

@a=qw(a b c d e f g h );
$b='d';
sub half_search{
 my($ref,$key,$low,$high)=@_;
 $high=@{$ref}-1 unless $high;
 if ($key lt $ref->[$low] or $key gt $ref->[$high]){
 print "not exists !!!\n";
 last;
 }
 if ($ref->[$low] gt $ref->[$high]){
 print "not sort array !!!\n" ;
 last;
 }
 $mid=int (($low+$high)/2);
 if ($ref->[$mid] eq $key) { 
 return $mid+1;
 }
 elsif ($ref->[$mid] gt $key) {
 &half_search($ref,$key,$mid+1,$high);
 }
 else{
 &half_search($ref,$key,$low,$mid-1);
 }
}
print &half_search(\@a,$b);

本人亲测可用,就不贴图啦!

19

个人网站的计划

转录组方向:

数据来源是NCBI里面的一个文献

网站的计划32

其中转录组方向的那些软件流程大多已经跑完了,大家可以见我的转录组总结。

trinity,tophat,cufflinks,RseQC,RNAseq,GOseq,MISO,RSEM,khmer,screed,trimmomatic,transDecoder,vast-tools,picard-tools,htseq,cuffdiff,edgeR,DEseq,funnet,davidgo,wego,kobas,KEGG,Amigo,go

 

基因组方向:

数据来源是strawberry草莓的文献

网站的计划282

velvet,SOAPdenovo2,repeatmasker,repeatscount,piler,

Chip-seq方向:

网站的计划348

这个群里有高手说要跟我合作,他来帮我写,希望是真的!

免疫组库方向:

网站的计划385

这个其实没有成熟软件,也就是一个igblastn, 然后是IMGT数据库,但是是我主打的产品,所以我会详细介绍一下。

全外显子组方向:

网站的计划455

这方面我不是很懂,。好像主要就是snp-calling

Snp-calling方向:

这个我准备自己写软件,不仅仅是用别人的软,它的数据本身也是前面几个方向的数据

bwa,bowtie,samtools,GATK,VarScan.jar,annovar

进化方向:

数据就是基因组数据

orthMCL,inparanoid, clustw,muscle,MAFFT,quickparanoid,blast2go,RAxML,phyML

 

 

19

Linux基础之shell脚本的批处理

脚本类似于下面的样子,大家可以读懂之后就仿写

for i in *sra

do

echo $i

/home/jmzeng/bio-soft/sratoolkit.2.3.5-2-ubuntu64/bin/fastq-dump --split-3 $i

Done

这个脚本是把当前目录下所有的NCBI下载的sra文件都加压开来成测序fastq格式文件

有这些数据,分布在不同的目录,如果是写命令一个个文件处理,很麻烦,如果有几百个那就更麻烦了,所以需要用shell脚本

Linux基础之shell脚本的批处理254

这样只需要bash这个脚本即可一次性处理所有的数据

Linux基础之shell脚本的批处理282

还有很多类似的脚本,非常简单的

for i in *fq

do

echo $i

bowtie2 -p 13 -x ../../RNA.fa -U $i -S  $i.sam

done

 

for i in */accepted_hits.bam

do

echo $i

out=`echo $i |cut -d'/' -f 1`_clout

samtools mpileup -guSDf  /home/immune/refer_genome/hg19/hg19.fa $i  | bcftools view -cvNg - >snp-vcf/$out.vcf

done

 

 

while read id

do

echo $id

out=`echo $id |cut -d'/' -f 2`

reads=`echo $id |cut -d'/' -f 3|sed 's/\r//g'`

tophat2 -p 13 -o $out /home/immune/refer_genome/hg19/hg19 $reads

done <$1

 

等等

19

转录组总结

网站成立也快一个月了,总算是完全搞定了生信领域的一个方向,当然,只是在菜鸟层面上的搞定,还有很多深层次的应用及挖掘,仅仅是我所讲解的这些软件也有多如羊毛的参数可以变幻,复杂的很。其实我最擅长的并不是转录组,但是因为一些特殊的原因,我恰好做了三个转录组项目,所以手头上关于它的资料比较多,就分享给大家啦!稍后我会列一个网站更新计划,就好谈到我所擅长的基因组及免疫组库。我这里简单对转录组做一个总结:

首先当然是我的转录组分类网站啦

转录组总结317

http://www.bio-info-trainee.com/?cat=18

 

 

同样的我用脚本总结一下给大家

转录组总结335

 

http://www.bio-info-trainee.com/?p=370阅读更多关于《转录组-GO和KEGG富集的R包clusterProfiler》

http://www.bio-info-trainee.com/?p=359阅读更多关于《转录组-GO通路富集-WEGO网站使用》

http://www.bio-info-trainee.com/?p=346阅读更多关于《转录组-TransDecoder-对trinity结果进行注释》

http://www.bio-info-trainee.com/?p=271阅读更多关于《转录组cummeRbund操作笔记》

http://www.bio-info-trainee.com/?p=255阅读更多关于《转录组edgeR分析差异基因》

http://www.bio-info-trainee.com/?p=244阅读更多关于《转录组HTseq对基因表达量进行计数》

http://www.bio-info-trainee.com/?p=166阅读更多关于《转录组cufflinks套装的使用》

http://www.bio-info-trainee.com/?p=156阅读更多关于《转录组比对软件tophat的使用》

http://www.bio-info-trainee.com/?p=125阅读更多关于《Trinity进行转录组组装的使用说明》

http://www.bio-info-trainee.com/?p=113阅读更多关于《RSeQC对 RNA-seq数据质控》

同时我也讲了如何下载数据

http://www.bio-info-trainee.com/?p=32

转录组总结1058

 

原始SRA数据首先用SRAtoolkit数据解压,然后进行过滤,评估质量,然后trinity组装,然后对组装好的进行注释,然后走另一条路进行差异基因,差异基因有tophat+cufflinks+cummeRbund,也有HTseq 和edgeR等等,然后是GO和KEGG通路注释,等等。

在我的群里面共享了所有的代码及帖子内容,欢迎加群201161227,生信菜鸟团!

http://www.bio-info-trainee.com/?p=1

线下交流-生物信息学
同时欢迎下载使用我的手机安卓APP

http://www.cutt.com/app/down/840375

19

转录组-GO和KEGG富集的R包clusterProfiler

PS: 请不要在问我关于这个包的任何问题,直接联系Y叔,我就两年前用过一次而已,再也没有用过。

Y叔的包更新太频繁了,这个教程已经作废,请不要再照抄了,可以去我们论坛看新的教程:http://www.biotrainee.com/thread-1084-1-1.html

一:下载安装该R包

clusterProfiler是业界很出名的YGC写的R包,非常通俗易懂,也很好用,可以直接根据cuffdiff等找差异的软件找出的差异基因entrez ID号直接做好富集的所有内容; Continue reading