24

Genomemapper软件使用说明书

 我以前一直以为有了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

21

SAMStat软件使用说明书

这个软件是对我们的比对结果(通常是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

16

Annovar使用记录

至于如何安装该软件,请见上一个教程

一.首先把snp-calling步骤的VCF文件转为annovar软件要求的格式

convert2annovar.pl   -format vcf4   12.vcf >12.annovar

Annovar使用记录108

二.进行注释

命令行参数比较多,还是用脚本来运行

# define path

infolder=/home/jmzeng/hoston/diff

outfolder=$infolder

annovardb=/home/jmzeng/bio-soft/annovar/humandb

# start annotating

/home/jmzeng/bio-soft/annovar/annotate_variation.pl \

--buildver hg19 \

--geneanno \

--outfile ${outfolder}/12.anno \

${infolder}/12.annovar  \

${annovardb}

三.输出结果解读

2.6M Apr 14 22:32 12.anno.exonic_variant_function

1.9K Apr 14 22:32 12.anno.log

1.3M Apr 14 22:32 12.anno.variant_function

重点是后缀为exonic_variant_function,这个文件对每一个vcf的突变都进行了注释。

Annovar使用记录617

这个结果就可以用来解析了,可以根据实验设计来找到自己感兴趣的突变。

第5.6列是染色体及pos坐标

第4列信息非常复杂,是突变的注释

第12列是测序深度,一般要大于20

我这里是先把注释文件转换成以下格式

location:chr1:874467 SAMD11:NM_152486:exon6:c.G478A:p.D160N

location:chr1:888639 NOC2L:NM_015658:exon9:c.A918G:p.E306E

location:chr1:888659 NOC2L:NM_015658:exon9:c.A898G:p.I300V

location:chr1:916549 PERM1:NM_001291367:exon2:c.T58C:p.W20R

location:chr1:949608 ISG15:NM_005101:exon2:c.G248A:p.S83N

location:chr1:980552 AGRN:NM_198576:exon13:c.G2266A:p.A756T

location:chr1:1114699 TTLL10:NM_001130045:exon4:c.G104A:p.R35Q

location:chr1:1158631 SDF4:NM_016176:exon4:c.T570C:p.D190D

location:chr1:1158631 SDF4:NM_016547:exon4:c.T570C:p.D190D

location:chr1:1164073 SDF4:NM_016176:exon2:c.C101T:p.A34V

然后比较两个文件,取不同的突变来格式化输出。

29

NGS QC Toolkit 对测序reads进行简单过滤

这个软件其实我真心不需要讲些什么了,它的官网写的太好了,简直就是软件说明书的典范

http://www.nipgr.res.in/ngsqctoolkit.html

它列出了它的几个功能模块,还给出了下载地址,还给出了说明文档,下载压缩包,解压即可使用啦

更重要的是给出了测试数据和测试的结果,而且还专门测试了不同测序平台及不同的测序策略的使用说明

 

NGS QC Toolkit 对测序reads进行简单过滤264

里面就是一些perl测序,其实自己都可以写的,分成了四大类。

其中统计的那个平均测序质量,我在前面仿写fastqc就写过,至于那个统计N50,更是生信常用的脚本。

但是大家可以看看这个perl程序来学perl语言,蛮不错的这些程序,都写的很标准。

比如那个TrimmingReads.pl

NGS QC Toolkit 对测序reads进行简单过滤576

可以根据四个参数来选择性的对我们的原始reads进行过滤,当然很多其它的程序也有类似的功能,它的参数分别是铲掉5端的几个碱基或者3端的,或者根据测序质量来切除碱基,或者根据reads长度来取舍,都是挺实用的功能。但是我一般用LengthSort和DynamicTrim那两个程序,原因很简单,我老师是这样用的,所以我习惯了,哈哈

29

Samtools安装及使用

一、下载安装该软件。

网上可以搜索到下载地址,解压之后make即可

一般都会报错

In file included from bam_cat.c:41:0:

htslib-1.1/htslib/bgzf.h:34:18: fatal error: zlib.h: No such file or directory

 #include <zlib.h>

^

compilation terminated.

make: *** [bam_cat.o] Error 1

Samtools安装及使用265

然后,居然就通过了,晕。有时候我实在是搞不定linux系统一些具体的原理,但是反正就是能用!学会搜索,学会试错即可。

直到两年后我才理解(linux下 的软件安装需要指定路径,而且是自己有权限的路径,2016年11月23日10:12:11),比如安装下面的方式来安装软件:

mkdir -p ~/biosoft/myBin
echo 'export PATH=/home/jianmingzeng/biosoft/myBin/bin:$PATH' >>~/.bashrc
source ~/.bashrc
cd ~/biosoft
mkdir cmake && cd cmake
wget http://cmake.org/files/v3.3/cmake-3.3.2.tar.gz
tar xvfz cmake-3.3.2.tar.gz
cd cmake-3.3.2
./configure --prefix=/home/jianmingzeng/biosoft/myBin  ## 这里非常重要
make
make install

但是有些电脑会报另外一个错

#include <curses.h>

^

compilation terminated.

make: *** [bam_tview_curses.o] Error 1

我也顺便解决一下,因为以前我的服务器遇到过,也是很纠结的。

sudo apt-get install libncurses5-dev

二.准备数据及使用,见我的snp-caling流程

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

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

因为这个软件都是与bwa和bowtie等能产生sam文件的软件合作才能使用。

其中这个软件参数还是蛮多的,但是常用的就那么几个,网上也很容易找到教程

简单附上一点资料

 

 

samtools是一个用于操作sam和bam文件的工具合集。包含有许多命令。以下是常用命令的介绍

1. view

view命令的主要功能是:将sam文件转换成bam文件;然后对bam文件进行各种操作,比如数据的排序(不属于本命令的功能)和提取(这些操作是对bam文件进行的,因而当输入为sam文件的时候,不能进行该操作);最后将排序或提取得到的数据输出为bam或sam(默认的)格式。

bam文件优点:bam文件为二进制文件,占用的磁盘空间比sam文本文件小;利用bam二进制文件的运算速度快。

view命令中,对sam文件头部的输入(-t或-T)和输出(-h)是单独的一些参数来控制的。

Usage: samtools view [options] <in.bam>|<in.sam> [region1 [...]]默认情况下不加 region,则是输出所有的 region. Options:

-b       output BAM                  默认下输出是 SAM 格式文件,该参数设置输出 BAM 格式         -h       print header for the SAM output                  默认下输出的 sam 格式文件不带 header,该参数设定输出sam文件时带 header 信息         -H       print header only (no alignments)         -S       input is SAM                  默认下输入是 BAM 文件,若是输入是 SAM 文件,则最好加该参数,否则有时候会报错。

例子:

#将sam文件转换成bam文件$ samtools view -bS abc.sam > abc.bam$ samtools view -b -S abc.sam -o abc.bam

#提取比对到参考序列上的比对结果$ samtools view -bF 4 abc.bam > abc.F.bam #提取paired reads中两条reads都比对到参考序列上的比对结果,只需要把两个4+8的值12作为过滤参数即可$ samtools view -bF 12 abc.bam > abc.F12.bam #提取没有比对到参考序列上的比对结果$ samtools view -bf 4 abc.bam > abc.f.bam #提取bam文件中比对到caffold1上的比对结果,并保存到sam文件格式$ samtools view abc.bam scaffold1 > scaffold1.sam #提取scaffold1上能比对到30k到100k区域的比对结果$ samtools view abc.bam scaffold1:30000-100000 > scaffold1_30k-100k.sam #根据fasta文件,将 header 加入到 sam 或 bam 文件中$ samtools view -T genome.fasta -h scaffold1.sam > scaffold1.h.sam

2. sort

sort对bam文件进行排序。

Usage: samtools sort [-n] [-m <maxMem>] <in.bam> <out.prefix>  -m 参数默认下是 500,000,000 即500M(不支持K,M,G等缩写)。对于处理大数据时,如果内存够用,则设置大点的值,以节约时间。-n 设定排序方式按short reads的ID排序。默认下是按序列在fasta文件中的顺序(即header)和序列从左往右的位点排序。

例子:

$ samtools sort abc.bam abc.sort$ samtools view abc.sort.bam | less -S

3.merge

将2个或2个以上的已经sort了的bam文件融合成一个bam文件。融合后的文件不需要则是已经sort过了的。

Usage:   samtools merge [-nr] [-h inh.sam] <out.bam> <in1.bam> <in2.bam>[...] Options: -n       sort by read names         -r       attach RG tag (inferred from file names)         -u       uncompressed BAM output         -f       overwrite the output BAM if exist         -1       compress level 1         -R STR   merge file in the specified region STR [all]         -h FILE  copy the header in FILE to <out.bam> [in1.bam] Note: Samtools' merge does not reconstruct the @RG dictionary in the header. Users      must provide the correct header with -h, or uses Picard which properly maintains      the header dictionary in merging.

4.index

必须对bam文件进行默认情况下的排序后,才能进行index。否则会报错。

建立索引后将产生后缀为.bai的文件,用于快速的随机处理。很多情况下需要有bai文件的存在,特别是显示序列比对情况下。比如samtool的tview命令就需要;gbrowse2显示reads的比对图形的时候也需要。

Usage: samtools index <in.bam> [out.index]

例子:

#以下两种命令结果一样$ samtools index abc.sort.bam$ samtools index abc.sort.bam abc.sort.bam.bai

5. faidx

对fasta文件建立索引,生成的索引文件以.fai后缀结尾。该命令也能依据索引文件快速提取fasta文件中的某一条(子)序列

Usage: samtools faidx <in.bam> [ [...]] 对基因组文件建立索引$ samtools faidx genome.fasta#生成了索引文件genome.fasta.fai,是一个文本文件,分成了5列。第一列是子序列的名称;第二列是子序列的长度;个人认为“第三列是序列所在的位置”,因为该数字从上往下逐渐变大,最后的数字是genome.fasta文件的大小;第4和5列不知是啥意思。于是通过此文件,可以定位子序列在fasta文件在磁盘上的存放位置,直接快速调出子序列。 #由于有索引文件,可以使用以下命令很快从基因组中提取到fasta格式的子序列$ samtools faidx genome.fasta scffold_10 > scaffold_10.fasta

6. tview

tview能直观的显示出reads比对基因组的情况,和基因组浏览器有点类似。

Usage: samtools tview <aln.bam> [ref.fasta] 当给出参考基因组的时候,会在第一排显示参考基因组的序列,否则,第一排全用N表示。按下 g ,则提示输入要到达基因组的某一个位点。例子“scaffold_10:1000"表示到达第10号scaffold的第1000个碱基位点处。使用H(左)J(上)K(下)L(右)移动显示界面。大写字母移动快,小写字母移动慢。使用空格建向左快速移动(和 L 类似),使用Backspace键向左快速移动(和 H 类似)。Ctrl+H 向左移动1kb碱基距离; Ctrl+L 向右移动1kb碱基距离可以用颜色标注比对质量,碱基质量,核苷酸等。30~40的碱基质量或比对质量使用白色表示;20~30黄色;10~20绿色;0~10蓝色。使用点号'.'切换显示碱基和点号;使用r切换显示read name等还有很多其它的使用说明,具体按 ? 键来查看。

 

参考:samtools的说明文档:http://samtools.sourceforge.net/samtools.shtml

http://www.plob.org/2014/01/26/7112.html

 

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

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.
15

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

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

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

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

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

Continue reading