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

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

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

 

等等

18

查某个基因家族在某物种的具体信息

查某个基因家族在某物种的具体信息

我很伤心,不知道是不是我写的教程还是不够人性化,一个朋友在群里面问如何知道NAC基因家族在拟南芥里面的105个基因信息,我随便给他示范了一下在人类里面如何找,希望他能触类旁通,结果他不会linux,啥生信基础都没有,我只会诱导他简单学习一下,希望他至少明白什么的taxid。所以我给了他我之前写的教程,只希望他告诉我拟南芥的taxid我就帮他把那105个基因找出来。 Continue reading

18

生信常用论坛seq-answer里面所有帖子爬取

生信常用论坛seq-answer里面所有帖子爬取

这个是爬虫专题第二集,主要讲如何分析seq-answer这个网站并爬去所有的帖子列表,及标签列表等等,前提是读者必须掌握perl,然后学习perl的LWP模块,可以考虑打印那本书读读,挺有用的!

其实爬虫是个人兴趣啦,跟这个网站没多少关系,本来一个个下载,傻瓜式的重复也能达到目的。我只是觉得这样很有技术范,哈哈,如何大家不想做傻瓜式的操作可以自己学习学习,如果不懂也可以问问我!

http://seqanswers.com/这个是主页

http://seqanswers.com/forums/forumdisplay.php?f=18 这个共570个页面需要爬取

其中f=18 代表我们要爬去的bioinformatics板块里面的内容

http://seqanswers.com/forums/forumdisplay.php?f=18&order=desc&page=1

http://seqanswers.com/forums/forumdisplay.php?f=18&order=desc&page=570

<tbody id="threadbits_forum_18">这个里面包围这很多<tr>对,

前五个<tr>对可以跳过,里面的内容不需要

Continue reading

18

生信常用论坛bio-star里面所有帖子爬取

生信常用论坛bio-star里面所有帖子爬取

这个是爬虫专题第一集,主要讲如何分析bio-star这个网站并爬去所有的帖子列表,及标签列表等等,前提是读者必须掌握perl,然后学习perl的LWP模块,可以考虑打印那本书读读,挺有用的!

http://seqanswers.com/ 这个是首页

http://seqanswers.com/forums/forumdisplay.php?f=18 这个共570个页面需要爬取

http://seqanswers.com/forums/forumdisplay.php?f=18&order=desc&page=1

http://seqanswers.com/forums/forumdisplay.php?f=18&order=desc&page=570

<tbody id="threadbits_forum_18">这个里面包围这很多<tr>对,

前五个<tr>对可以跳过,里面的内容不需要

生信常用论坛bio_star462

Continue reading

17

转录组cummeRbund操作笔记

转录组cummeRbund操作笔记

这是跟tophat和cufflinks套装紧密搭配使用的一个R包,能出大部分文章要求的标准化图片。

一:安装并加装该R包

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

转录组cummeRbund操作笔记220

然后把cuffdiff输出的文件目录拷贝到R的工作目录,或者自己设置工作目录

 

二:读取FN目录下面的所有文件。

转录组cummeRbund操作笔记239

可以看到把cuffdiff下面的文件夹所有的文件都读取到了,里面有如下文件,包括genes,isoforms,cds,tss这四种差异情况都读取了。

转录组cummeRbund操作笔记316

 

三:表达水平分布图

转录组cummeRbund操作笔记328

转录组cummeRbund操作笔记330
四、表达水平箱线图

csBoxplot(genes(cuff_data))

转录组cummeRbund操作笔记371
五、画基因表达差异热图

转录组cummeRbund操作笔记386

画出热图如下

转录组cummeRbund操作笔记396

 

六、得到差异的genes,isoforms,TSS,CDS等等

 

  • 得到上调下调基因列表

diffData <- diffData(myGenes )

转录组cummeRbund操作笔记430

转录组cummeRbund操作笔记474

 

只有一百个有表达差异的基因

转录组cummeRbund操作笔记490

 

 

最后贴出一个综合性的代码,算了,太浪费空间了,把整个空间搞得不好看,就不贴了。

这个代码可以自动运行出图;

转录组cummeRbund操作笔记3781

16

转录组edgeR分析差异基因

转录组edgeR分析差异基因

edgeR是一个研究重复计数数据差异表达的Bioconductor软件包。一个过度离散的泊松模型被用于说明生物学可变性和技术可变性。经验贝叶斯方法被用于减轻跨转录本的过度离散程度,改进了推断的可靠性。该方法甚至能够用最小重复水平使用,只要至少一个表型或实验条件是重复的。该软件可能具有测序数据之外的其他应用,例如蛋白质组多肽计数数据。可用性:程序包在遵循LGPL许可证下可以从Bioconductor网站。

一:下载安装该软件

下载安装edgeR这个R包,因为这是一次讲R包的下载,我就啰嗦一点,这种生物信息学的包不同于普通的R包,是需要用biocLite来安装的,命令如下

转录组edgeR分析差异基因304

 

Continue reading

15

仿写fastqc软件的一些功能-R代码

仿写fastqc软件的一些功能(下)

文件来自于上面perl代码的输出文件,好像算法有点问题,26G的文件居然处理近一个小时才出数据!

仿写fastqc软件的一些功能-下-R代码263

R语言本身自带的画图工具都很丑,懒得说了,可以用ggplot2来重新画一个,不是项目要求没有报酬我就懒得画了,大家面前看看画图原理即可。

Continue reading

14

仿写fastqc软件的部分功能-perl代码

  仿写fastqc软件的部分功能(上)

前面我们介绍了fastqc这个软件的使用方法 http://www.bio-info-trainee.com/?p=95 ,这是一个java软件,但是有些人服务器没有配置好这个java环境,导致无法使用,这里我贴出几个perl代码,也能实现fastqc的部分功能

统一测试文件是illumina的phred33格式的fastq文件,共100000/4=25000条reads,读长都是101个碱基

程序名-fastq2quality.pl

使用命令:perl fastq2quality.pl SRR504517_1.fastq >quality.txt

功能: 把fastq格式的每条原始reads的第四行ascii码质量值,转换为Q值并输出一个矩阵,有多少条reads就有多少行,每条reads的碱基数就是列数。

Continue reading