10

一步一步运行软件系列合集

这些是很久以前写的一些教程,是关于进化树构建和全基因组关联分析的!

gwas-plink分析教程.pdf
plink的统计基础.ppt
一步一步构建系统进化树.pdf
一步一步运行blast.pdf
一步一步运行inparanoid蛋白聚类.pdf
一步一步运行PLINK-part1.pdf
一步一步运行plink-part2.pdf
用PhyML构建系统发育树.pptx
进化树的构建分子原理.pdf

都在云盘(http://pan.baidu.com/s/1jIvwRD8 )里面,群空间(201161227)里面也有!

暂时应该不会写这些教程了,因为没有项目,实在没有动力去做那么多事情

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.
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的原理之一!