三 25

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

Posted on 2015年3月25日 by ulwvfje

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

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

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

这样我们的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

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

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

但是由于我代码里面选取的是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安装使用

Posted on 2015年3月25日 by ulwvfje

一．下载并安装这个软件

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

解压进入目录

首先make

然后make install即可

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

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

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程序

二．准备测试数据

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

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

三，参考命令

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

官网给出的步骤如下

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

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序列，几个小时就跑完啦

四个步骤都有输出数据

好像组装效果惨不忍睹呀！共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

草莓基因组数据预处理

Posted on 2015年3月24日 by ulwvfje

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

首先跑fastqc

打开一个个看结果

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

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

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

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

for id in *fastq

echo $id

perl DynamicTrim.pl -454 $id

done

for id in *trimmed

echo $id

perl LengthSort.pl $id

done

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

处理完毕后的数据如下：

三 24

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

Posted on 2015年3月24日 by ulwvfje

一．下载该软件

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

下载解压开

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

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

这两

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

二．准备数据。

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

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

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

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

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

for id in *fastq

echo $id

perl DynamicTrim.pl -454 $id

done

for id in *trimmed

echo $id

perl LengthSort.pl $id

done

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

几个小时完毕之后

查看，产出文件如下

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

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

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

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

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

三 24

旧版本blast详解

Posted on 2015年3月24日 by ulwvfje

其实我现在一般都用的是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的突变相关的一些数据解释。

Posted on 2015年3月23日 by ulwvfje

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

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

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

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

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,

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

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, 在基因的四个外显子的坐标

DDX11L1 基因的通俗英文名

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

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

三 23

用annovar对snp进行注释

Posted on 2015年3月23日 by ulwvfje

一、下载及安装软件

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

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

参考：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

变成了空格分隔的文件

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

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

四，输出文件解读

三 23

Snp-calling流程（BWA+SAMTOOLS+BCFTOOLS）

Posted on 2015年3月23日 by ulwvfje

比对可以选择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

echo $id

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

done <$1

然后sampe脚本批量处理

==> bwa_sampe.sh <==

while read id

echo $id

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

done <$1

然后是samtools的脚本

==> samtools.sh <==

while read id

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

echo $id

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

done <$1

==> mpileup.sh <==

while read id

echo $id

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

done <$1

三 22

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

Posted on 2015年3月22日 by ulwvfje

需要插件和自己修改主题下面的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”插件已经处于启用状态（插件名下是“停用”链接）。

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

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

其中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

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

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

三 21

Hg19基因组的分析

Posted on 2015年3月21日 by ulwvfje

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

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基因组是3G的大小，因为一个英文字符是一个字节，所以也是30亿bp的碱基。

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

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

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

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

awk 'END { print NR }' hg19.fa

wc -l hg19.fa

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

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

[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的词频（我不知道生物信息学里面的专业描述词语是什么）

我的程序耗费了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字符串短序列在人的基因组里面都可以找到！！！

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

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

三 21

积累的一些perl代码分享

Posted on 2015年3月21日 by ulwvfje

以前的一下perl代码分享

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

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

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

Continue reading →

三 21

Linux服务器基础知识

Posted on 2015年3月21日 by ulwvfje

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

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

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

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

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

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

free -g

df -h

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

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

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

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

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

三 20

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

Posted on 2015年3月20日 by ulwvfje

首先，什么是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算法详解之一建立索引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简单使用

Posted on 2015年3月20日 by ulwvfje

首先进入bowtie的主页，千万要谷歌！！！

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

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

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

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

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

二、准备数据

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

三、运行命令

分为两步，首先索引，然后比对！！！

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

然后你的目录就产生了六个索引文件，我给索引取名是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实现二分法查找

Posted on 2015年3月20日 by ulwvfje

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

个人网站的计划

Posted on 2015年3月19日 by ulwvfje

转录组方向：

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

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

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草莓的文献

velvet，SOAPdenovo2，repeatmasker,repeatscount,piler，

Chip-seq方向：

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

免疫组库方向：

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

全外显子组方向：

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

Snp-calling方向：

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

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

进化方向：

数据就是基因组数据

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

三 19

Linux基础之shell脚本的批处理

Posted on 2015年3月19日 by ulwvfje

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

for i in *sra

echo $i

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

Done

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

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

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

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

for i in *fq

echo $i

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

done

for i in */accepted_hits.bam

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

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

转录组总结

Posted on 2015年3月19日 by ulwvfje

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

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

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

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

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

原始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

Posted on 2015年3月19日 by ulwvfje

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

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

一：下载安装该R包

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

三 19

转录组-GO通路富集-WEGO网站使用

Posted on 2015年3月19日 by ulwvfje

一，所谓的网站，其实就是一个网页版的可视化软件接口而已

看看网站主页，看看它需要什么数据

http://wego.genomics.org.cn/cgi-bin/wego/index.pl

二，所需要的数据

1，human.all.go.entrez，需要自己制作，每个基因名entrez ID号，对应着一堆GO通路，人有两万多个基因，所以应该有两万多行的文件。

2，差异基因的GO通路，需要用cuffdiff得到差异基因名，然后用然后用脚本做成下面的样子。记住，上面的那个人类的背景GO文件也是一样的格式，基因名是entrez ID号，与GO通路用制表符隔开，然后每个基因所对应的GO直接用空格隔开。格式要求很准确才行。

三，上传数据，出图

点击plot画图即可，就可以出来了一个GO通路富集图

顺便贴上wego上传数据制作的几个脚本，脚本这种东西都很难看，随便意思一下啦，用一下脚本处理就可以得到wego需要上传的数据了

1，得到差异基因名，并且转换为entrez ID号
grep yes gene_exp.diff |cut -f 3 |sort -u >diff.gene.name
cat diff.gene.name ../Homo_sapiens.gene_info |perl -alne '{$hash{$_}=1;print $F[1] if exists $hash{$F[2]}}' |sort -u >diff.gene.entrez
2，根据找到的差异基因的entrez ID号来找到它的GO信号，输出文件给wego网站
cat diff.gene.entrez ../gene2go |perl -alne '{$hash{$_}=1;print "$F[1]\t$F[2]" if exists $hash{$F[1]}}' |perl -alne '{$hash{$F[0]}.="$F[1] "}END{print "$_\t$hash{$_}" foreach keys %hash}' >diff.gene.entrez.go
3，得到entrez ID号跟ensembl ID号的转换hash表
perl -alne '{if (/Ensembl:(ENSG\d+)/) {print "$1=>$F[1]"} }' Homo_sapiens.gene_info >entrez.ensembl
4，得到人类entrez ID的go背景
grep '^9606' gene2go |perl -alne '{$hash{$F[1]}.="$F[2] "}END{print "$_\t$hash{$_}" foreach sort keys %hash}' >human.all.go.entrez
5，把人类entrez ID的go背景转换成ensembl的go背景
cat entrez.ensembl human.all.go.entrez |perl -F"=>" -alne '{$hash{$F[1]}=$F[0];print "$hash{$F[0]}\t$F[1]" if exists $hash{$F[0]}}' >human.all.go.ensembl

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

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

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

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

Page 107 of 109« First ‹ Previous 104 105 106107108 109 Next ›

2024年5月
一	二	三	四	五	六	日
« 四
		1	2	3	4	5
6	7	8	9	10	11	12
13	14	15	16	17	18	19
20	21	22	23	24	25	26
27	28	29	30	31