01

重复序列屏蔽第一讲RepeatMasker的一些参数调试

这是很久以前的一篇文章,我先贴出来给大家看看,然后讲一个实例

一:RepeatMasker的一些参数运行结果比较

从ncbi随便下载的zebrafish的一条sequence.fasta

不加上任何参数跑出来结果是 RepeatMasker   sequence.fasta

repeat-masker参数摸索138

加上物种的参数之后跑出来是: RepeatMasker -species Danio  sequence.fasta

效果里面出来了,之前得到的重复序列不到10%,这次可以达到70%以上,所以必须得选好对应的物种,这样才不会错过那么多要找的重复序列

repeat-masker参数摸索267

repeat-masker参数摸索269

再加上-low这个参数是 RepeatMasker -species Danio -low  sequence.fasta

感觉没有改变多少,就少了几个

repeat-masker参数摸索349 repeat-masker参数摸索351

比较-div参数:RepeatMasker -species Danio  sequence.fasta

RepeatMasker -species Danio -div 10  sequence.fasta

repeat-masker参数摸索459

而加上-div 10之后

repeat-masker参数摸索475

第二列小于10%的全部被剔除掉了

输出参数,本来应该是用N把重复区域屏蔽掉的

repeat-masker参数摸索518

但是如果加上参数-x,原来输出是N的地方就都变成了X,感觉这个参数没啥子意义。

repeat-masker参数摸索560

还有一些类似的参数,意义也不大,加上-xsmall,就是把重复区域用小写字母,不再需要N来掩盖了

repeat-masker参数摸索613

如果加上-a这个参数,就多了一个文件

repeat-masker参数摸索637

查看可知其内容是

repeat-masker参数摸索648

The  alignments are in the cross_match/SWAT format, in which mismatches rather than matches are indicated: transitions

with an i and  transversions with a v. Note it exists some differences between the  alignment file and the map fi le.

The map fi le is produced by  ProcessRepeats that the main task is to defragment the original  map file and the alignment fi le is created from the original map fi le:  the difference between them comes from the defragmented hits.
如果加上-poly,也会多出一个文件

repeat-masker参数摸索1139

查看,可知其单独列出了微卫星的表格

repeat-masker参数摸索1159

The ‘-xm’, ‘-ace,’ and ‘-gff ’ options create an additional out put file in cross match, ACeDB, and Gene Feature Finding format  respectively.这几个参数都是为了生成适合其它处理的文件。

另外针对大文件的操作,可能需要-pa来设置运行速度,或者-s,-q,-qq

 

二:生成的文件的解释

会输出这些文件

repeat-masker参数摸索1387

1,。Out类文件

repeat-masker参数摸索1399

SW score 根据Smith-Waterman算法比对的分值 2555
Div% 比上区间与共有序列相比的替代率 5.7
Del% 在查询序列中碱基缺失的百分率(删除碱基) 0.0
Ins% 在repeat库序列中碱基缺失的百分率(插入碱基)  0.0
Query sequence 输入的待屏蔽重复的序列 gi|211853417|emb|CU633477.14|
Position begin 373
Position end  690
Query left 在查询序列中超出比上区域的碱基数

+= 比上了库中重复序列的正义链,如果是互补连用“c”表示

(50140)
Matching repeat 比上的重复序列的名称 C DNA13TA1a_DR
Repeat family(class) 比上的重复序列的类型   DNA/TcMar-Tc1
Position begin
Position end
Query left 比对区域距重复序列左端的碱基数
比对的顺序ID

3.cat文件基本类似于。Out文件
3。。Tbl类文件

repeat-masker参数摸索1917 repeat-masker参数摸索1919
4.masked文件,就是找到的重复序列被N给代替了,或者用参数改变代替形式

polyout文件。就是单独列出了微卫星表格

Align文件,其实就是把之前的。Out文件的每一行记录单独拿出来再进行表格化解释

repeat-masker参数摸索2027

把373到690的核苷酸序列列出来,说明这个DNA13TA1a_DR 重复具体的意义

但是没看懂这个i,v是什么意思

 

结果比较

从ncbi随便下载的zebrafish的一条sequence.fasta

不加上任何参数跑出来结果是 RepeatMasker   sequence.fasta

 

加上物种的参数之后跑出来是: RepeatMasker -species Danio  sequence.fasta

效果里面出来了,之前得到的重复序列不到10%,这次可以达到70%以上,所以必须得选好对应的物种,这样才不会错过那么多要找的重复序列

01

Perl及R及python模块碎碎念

老实说,模块其实是一个很讨厌的东西,但是它也实实在在的节省了我们很多时间,也符合我的理念:避免重复造轮子!此教程可能过期了,请直接看最新版(perl模块安装大全)

1,perl的那些模块

如果有root权限,用root权限

进入cpan然后install ExtUtils::Installed模块

这样就可以执行instmodsh这个脚本了,可以查看当前环境下所有的模块 Continue reading

01

R的包(package)

关于R语言包的一些操作,挺重要的!!!

R的包(package)通常有两种:
1 binary package:这种包属于即得即用型(ready-to-use),但是依赖与平台,即Win和Linux平台下不同。
2 Source package: 此类包可以跨平台使用,但用之前需要处理或者编译(compiled)。

以下一些常用的包相关的函数:
.libPaths():查看包的安装目录

ls('package:ggplot2')可以查看该包里面所有的函数
library():查看已经安装的包目录
library(mypackage):载入mypackage包

getOption("defaultPackages"):查看启动R时自动载入的包。
help(package = 'mypackage'):查看‘mypackage’的帮助
args(function):查看函数的参数
example(function):自动运行该函数帮助文档中的例子,很赞!
demo("package"):展示一些包中demostration,需要再看下??
vignette('mypackage'):有的包,特别是bioconductor的包有vignette,用函数查看
openVignette('mypackage'):这个函数也可以查看vignette,更好用一些
RSiteSearch("helpinfor"):搜索R网站上的“helpinfor”相关信息
help.start():查看已经安装包的详细HTML文档,这个命令非常爽
更新:
search():查看当前载入的包

sessionInfo():查看R中载入的包
methods():查看某个S3泛型函数中所有的方法或者一个类中所有的方法(S3:S version 3)

showMethods(class = "myClass"):查看S4类的方法

findMethods("myMethods"):查看method的代码

class(myObject):查看某个对象的类
getClass(“class/package”):查看某个class或者包的具体内容

getSlots("class"):查看某个class的slot

slotNames(MyObject):查看某个对象的slot。

可以使用Myobject@slotNames访问对象的slot值,这个@设计实在是太爽了,可以连续用。
查询包内信息:1. ?function/method:查看某个“函数”或者“方法”的详细内容
2. class?graph::graph:查看“组”的详细内容的一个例子。这个例子的来源是查询graph包时候,查看其中class的信息,输入??graph后出现一个graph::graph-class
ls("package:mypackage"):查看"mypackage"中的所有对象。

安装source package方法

1 在终端输入 # R CMD INSTALL /.../mypackage.tar.gz
使用此方法,需要解决包依赖问题,即安装此包所依赖的包,安装过程有提示

2 也可以使用R的install.packages()函数安装
回答:可以使用install.packages()函数安装,而且比较简便,即联网即可装,装了就可用。
# R
> install.packages('mypackage')

回答2:可以使用install.packages()安装本地下载的包,尤其适用于在服务器上安装包

$ R

> install.packages( c("XML_0.99-5.tar.gz", "http://www.cnblogs.com/Interfaces/Perl/RSPerl_0.8-0.tar.gz"), repos = NULL, configure.args = c(XML = '--with-xml-config=xml-config', RSPerl = "--with-modules='IO Fcntl'"))
3 Bioconductor的安装方法
> source("http://bioconductor.org/biocLite.R")
> biocLite("mypackage")

 

4 卸载package

remove.packages("mypackage")

 

5 查看R及其package的version

R version: version 或者 R.version

R package version:

 

6 更新包

update.packages( )  可以定期执行以下

 

7 使用别人安装的包

修改.bashrc文件,添加环境变量R的lib路径

export R_LIBS=/home/.../R/lib64/R/library

R中用.libPaths()函数查看lib路径,如果有多个lib,install.packages()默认是安装在第一个目录下

 

01

ubuntu服务器解决方案第十讲–虚拟机屏幕及联网设置

发现自己搞服务器遇到的困难还是蛮多的,所以记录了一下,给菜鸟们指个路。
ubuntu对生信菜鸟来说是最好用的linux服务器,没有之一,因为它有apt-get。

很多人可能并没有自己的服务器,那么就只能通过虚拟机来试试ubuntu啦

我想起来我以前玩虚拟机的时候遇到过一些困难,记录了一些,分享给大家, Continue reading

01

ubuntu服务器解决方案第九讲-mysql和apache的安装

发现自己搞服务器遇到的困难还是蛮多的,所以记录了一下,给菜鸟们指个路。
ubuntu对生信菜鸟来说是最好用的linux服务器,没有之一,因为它有apt-get。

很多时候大家的服务器可能并不是想联网,只是想玩一下,或者只是因为生信的某些软件要求数据库,所以大家可能会单独安装mysql,或者想学习perl的CGI模块,需要apache。

ubuntu上安装mysql

非常简单只需要几条命令就可以完成。

1. sudo apt-get install mysql-server

2. sudo apt-get install mysql-client

3.  sudo apt-get install libmysqlclient-dev

安装过程中会提示设置密码什么的,注意设置了不要忘了,安装完成之后可以使用如下命令来检查是否安装成功:

sudo netstat -tap | grep mysql

通过上述命令检查之后,如果看到有mysql 的socket处于 listen 状态则表示安装成功。

登陆mysql数据库可以通过如下命令:

mysql -u root -p

-u 表示选择登陆的用户名, -p 表示登陆的用户密码,上面命令输入之后会提示输入密码,此时输入密码就可以登录到mysql。

Ubuntu上安装Apache

Ubuntu上安装Apache,有两种方式:1 使用开发包的打包服务,例如使用apt-get命令;2 从源码构建Apache。本文章将详细描述这两种不同的安装方式。

方法一:使用开发包的打包服务——apt-get

安装apache,在命令行终端中输入一下命令:

$ sudo apt-get install apache2

安装完成后,重启apache服务,在命令行终端中输入一下命令:

$ sudo /etc/init.d/apache2 restart

可能会出现的问题1: NameVirtualHost *:80 has no VirtualHosts,

出现上述问题的原因:定义了多个NameVirtualHost,故将/etc/apache2/ports.conf中的NameVirtualHost *:80注释掉即可。

可能会出现的问题2: Could not reliably determine the server's fully qualified domain name, using 127.0.1.1 for ServerName

原因:

根据提示,无法可靠的确定服务器的有效域名,使用127.0.1.1作为服务器域名。应此,在下面的测试中,应该使用127.0.1.1,而不是127.0.0.1!

解决:

$ vim /etc/apache2/httpd.conf,在文件中添加:

ServerName localhost:80,再次重启apache2,就可以使用127.0.0.1来访问web服务器啦!

测试:

在浏览器里输入http://localhost或者是http://127.0.0.1,如果看到了It works!,那就说明Apache就成功的安装了,Apache的默认安装,会在/var下建立一个名为www的目录,这个就是Web目录了,所有要能 过浏览器访问的Web文件都要放到这个目录里。

01

ubuntu服务器解决方案第八讲–网络服务器配置lamp

发现自己搞服务器遇到的困难还是蛮多的,所以记录了一下,给菜鸟们指个路。
ubuntu对生信菜鸟来说是最好用的linux服务器,没有之一,因为它有apt-get。

主流的网络服务器配置就是linux+apache+mysql+php咯,简称LAMP

在ubuntu系统里面安装这个是非常easy的

sudo apt-get install apache2 mysql-server mysql-client php5 php5-gd php5-mysql Continue reading

01

ubuntu服务器解决方案第七讲-perl安装模块

此教程可能过期了,请直接看最新版(perl模块安装大全)

发现自己搞服务器遇到的困难还是蛮多的,所以记录了一下,给菜鸟们指个路。
ubuntu对生信菜鸟来说是最好用的linux服务器,没有之一,因为它有apt-get。

前面我简单写了一个perl的cpan安装模块,但是前些天突然发现有些perl模块在cpan里面找不到,所以又总结了一下不同的perl模块安装方法。 Continue reading

01

ubuntu服务器解决方案第六讲-添加环境变量

发现自己搞服务器遇到的困难还是蛮多的,所以记录了一下,给菜鸟们指个路。
ubuntu对生信菜鸟来说是最好用的linux服务器,没有之一,因为它有apt-get。

在我的第一讲里面,对JAVA的安装,其实就需要添加环境变量,大家可以回头看看!

添加PATH环境变量,第1种方法:

[root@lx_web_s1 ~]# export PATH=/usr/local/webserver/mysql/bin:$PATH

再次查看:

[root@lx_web_s1 ~]# echo $PATH

/usr/local/webserver/mysql/bin:/usr/local/webserver/mysql/bin/:/usr/kerberos/sbin:/usr/kerberos/bin:/usr/local/sbin:/usr/local/bin:/sbin:/bin:/usr/sbin:/usr/bin:/root/bin

说明添加PATH成功。

上述方法的PATH 在终端关闭 后就会消失。所以还是建议通过编辑/etc/profile来改PATH,也可以修改家目录下的.bashrc(即:~/.bashrc)。

第2种方法:需要管理员权限。

# vim /etc/profile

在最后,添加:

export PATH="/usr/local/webserver/mysql/bin:$PATH"

保存,退出,然后运行:

#source /etc/profile,不报错则成功。

01

ubuntu服务器解决方案第五讲-配置ssh供远程登录

发现自己搞服务器遇到的困难还是蛮多的,所以记录了一下,给菜鸟们指个路。
ubuntu对生信菜鸟来说是最好用的linux服务器,没有之一,因为它有apt-get。

同样,这个ssh也非常简单

sudo apt-get install openssh-server

SSH分客户端openssh-client和openssh-server

如果你只是想登陆别的机器的SSH只需要安装openssh-client(ubuntu有默认安装,如果没有则sudo

apt-get install openssh-client),如果要使本机开放SSH服务就需要安装openssh-server

sudo apt-get install openssh-server

然后确认sshserver是否启动了:

ps -e |grep ssh

如果看到sshd那说明ssh-server已经启动了。

如果没有则可以这样启动:sudo /etc/init.d/ssh start 或者 service ssh start

ssh-server配置文件位于/etc/ssh/sshd_config,在这里可以定义SSH的服务端口,默认端口是22,你可以自己定义成其他端口号,如222。

然后重启SSH服务:

sudo

/etc/init.d/ssh stop

sudo /etc/init.d/ssh start

然后使用以下方式登陆SSH:

ssh username@192.168.1.112 username为192.168.1.112 机器上的用户,需要输入密码。

我给七八个虚拟机都配置成功了,但是呢,偏偏别人的一个我始终不能解决,感觉linux里面的学问还是蛮多的

01

ubuntu服务器解决方案第四讲-输入法-中文

发现自己搞服务器遇到的困难还是蛮多的,所以记录了一下,给菜鸟们指个路。
ubuntu对生信菜鸟来说是最好用的linux服务器,没有之一,因为它有apt-get。

这个主要是针对有界面的服务器来说的,不是我们通常意义的ssh登陆,一般ssh登陆的可以把中文复制张贴进去即可。

Ubuntu上的输入法主要有小小输入平台(支持拼音/二笔/五笔等),Fcitx,Ibus,Scim等。其中Scim和Ibus是输入法框架。

在Ubuntu的中文系统中自带了中文输入法,通过Ctrl+Space可切换中英文输入法。这里我们主要说下Ubuntu英文系统中,中文输入法的安装。

安装输入法的第一步,是安装语言包。我们选择System Settings-->Language Support-->Install/Remove Languages,这里面可以选择简体中文

输入密码后,系统会安装简体中文语言包。

第二步,安装完毕后切换到终端,安装IBus框架,在终端输入以下命令:

sudo apt-get install ibus ibus-clutter ibus-gtk ibus-gtk3 ibus-qt4

启动IBus框架,在终端输入:

im-switch -s ibus

安装完IBus框架后注销系统,保证更改立即生效。

第三步:安装拼音引擎

有下面几种常用选择:

IBus拼音:sudo apt-get install ibus-pinyin

IBUS五笔:sudo apt-get install ibus-table-wubi

谷歌拼音输入法:sudo apt-get install ibus-googlepinyin

Sun拼音输入法:sudo apt-get install ibus-sunpinyin

第四步:设置IBus框架

终端输入ibus-setup 此时,IBus Preference设置被打开。我们在Input Method选项卡中,选择自己喜欢的输入方式,并配置自己喜欢的快捷键即可。

第五步:通常情况下,IBus图标(一个小键盘)会出现在桌面右上角的任务栏中。有时候这个图标会自行消失,可使用以下命令,找回消失的IBus图标:

ibus-daemon –drx

01

ubuntu服务器解决方案第三讲-perl最新版的安装

发现自己搞服务器遇到的困难还是蛮多的,所以记录了一下,给菜鸟们指个路。
ubuntu对生信菜鸟来说是最好用的linux服务器,没有之一,因为它有apt-get。

理论上perl是不需要更新,但是我就不巧碰到了这个情况,所以也记录一下

linux下升级系统默认安装的perl版本,不建议先rm

先下载tar.gz ...然後手动安装..default 安装到/usr/local/目录下..

然後修改/usr/bin/perl的symbolic link到/usr/local/bin/perl

下载方式不用说了吧,各显神通,笔者习惯用wget.

所以wget[url]http://www.cpan.org/src/perl-5.10.0.tar.gz[/url] .现在最新是5.20

下载完以后解压安装

#tar zxvf perl-5.10.0.tar.gz

#cd perl-5.10.0

#./Configure -des -Dprefix=/usr/local/perl

参数-Dprefix指定安装目录为/usr/local/perl

#make

#make test

#make install

如果这个过程没有错误的话,那么恭喜你安装完成了.是不是很简单?

接下来替换系统原有的perl,有最新的了咱就用嘛.

#mv /usr/bin/perl/ usr/bin/perl.bak

#ln -s /usr/local/perl/bin/perl/ usr/bin/perl

#perl –v

然后就可以了用它来安装一些其它你需要的perl模块了

#perl -MCPAN-e shell

第一次执行的话,会提示安装cpan并要求连接网络下载最新的模块列表.然后就可以安装东西了

cpan[1]> install DBI

01

ubuntu服务器解决方案第二讲-R程序包最新版的安装

发现自己搞服务器遇到的困难还是蛮多的,所以记录了一下,给菜鸟们指个路。
ubuntu对生信菜鸟来说是最好用的linux服务器,没有之一,因为它有apt-get。
首先,如果你的服务器里面有旧的R,那么删除Linux Ubuntu系统中原有的R软件包,代码如下:
sudo apt-get autoremove r-base-core # 删除系统中原有的R软件包
接下来,随便找到一个Ubuntu的软件源镜像
(http://mirror.bjtu.edu.cn/cran/bin/linux/ubuntu/ )
因为看我教程的大部分在国内,所以我拿北京大学举例子。
 当然,厦门大学也很赞,http://mirrors.xmu.edu.cn/CRAN/bin/linux/ubuntu/xenial/  不过,版本代号一定药搞清楚哦。
PS(2017/12/20, 没想到随便找到来一个,就挂掉了,唉:

厦门大学不再提供R语言镜像

)

Linux Ubuntu 12.04对应的名字是 precise,
ubuntu14.04,那么就应该是 trusty
ubuntu15.04 ,其代号为 vivid
Ubuntu 16.04 LTS,代号为Xenial Xerus(非洲的一种地松鼠),于UTC时间2016年4月21日正式发布。
比如我的Ubuntu 12.04就需要进入到 precise/目录,找到r-base-core相关的文件,发现有多个R的版本。
把这个软件源,增加到apt的sources.list文件中,代码如下:
deb http://mirror.bjtu.edu.cn/cran/bin/linux/ubuntu precise/
其余的ubuntu版本类似即可:
在/etc/apt/sources.list文件最下面,新加一行
~ sudo apt-get update # 更新源
~ sudo apt-get install r-base-core # 再次安装R语言软件包
~ R –version # 检查R的版本
这时我们就安装了最新的R语言版本—3.0.3版。
require(ggplot2)
Loading required package: ggplot2
Failed with error: ‘package ‘ggplot2’ was built before R 3.0.0: please re-install it’
这个失败原因是怎么回事?R 3.0.0 的问题吗?怎么解决?
R 2.x 升级 3.x 需要重新(编译)安装所有包:
update.packages(checkBuilt = TRUE, ask = FALSE)
当然如果你不是在R里面用install.package来安装包的话
你还需要
sudo apt-get install r-base-dev
这样你才能从源代码编译R的包
但是如果你导入的R源被你的服务器拒绝,你就惨了
The following signatures couldn't be verified because the public key is not
以下签名不能因为公钥未验证
01

ubuntu服务器解决方案第一讲-java安装

ubuntu对生信菜鸟来说是最好用的linux服务器,没有之一,因为它有apt-get。

1、JDK官网上http://www.oracle.com/technetwork/java/javase/downloads/index.html选择:

但是,如果你的服务器是64位的,请不要选择i586,选择你自己的机器对应的!

2、将打开终端,建立目录:

Sudo mkdir /usr/lib/java

3、将下载的 jdk-7u3-linux-i586.tar.gz移到这个文件夹下面并进行解压,改名字:

sudo mv jdk-7u3-linux-i586.tar.gz /usr/lib/java

sudo tar –xvf jdk-7u3-linux-i586.tar.gz

mv jdk1.7.0_03java-7-sun

4、修改环境变量:

在终端输入:vim /etc/profile

然后添加以下代码:

export JAVA_HOME=/usr/lib/java/jdk1.8.0_25

export JRE_HOME=${JAVA_HOME}/jre

export CLASSPATH=.:${JAVA_HOME}/lib:${JRE_HOME}/lib

export PATH=${JAVA_HOME}/bin:$PATH

保存之后,再运行下面命令更新电脑的配置文件

source /etc/profile  这个千万要记得!!!!

5、在终端中输入 java –version,显示:

jeydragon@jeydragon-VirtualBox:~$ java -version

java version "1.7.0_03"

Java(TM) SE Runtime Environment (build 1.7.0_03-b04)

Java HotSpot(TM) Client VM (build 22.1-b02, mixed mode)

表示安装成功

31

Ensembl数据库在线网页工具biomart简单教程

这个工具主要是针对不会bioperl不会API调取数据的生信纯菜鸟准备的,主要是方便大家批量研究某些感兴趣的基因,需要准备的数据就是基因名或者基因的ID号,能从该网站获取的资料非常多,可以是关于你的输入的基因名的各种数据库有的信息。

http://www.ensembl.org/biomart/

第一步:选取数据库,我一般选取人的

Ensembl数据库在线网页工具biomart简单教程243

第二步,选择上传数据的格式

Ensembl数据库在线网页工具biomart简单教程259

这个下拉框里面可以选取很多种格式,你随便张贴进去哪一种格式的基因ID都可以,也可以把做好的ID文件上传进去,批量获取基因信息。

Ensembl数据库在线网页工具biomart简单教程325

我这里输入的是几个免疫基因。

第三步,选择下载数据的格式

首先可以选择你上传的gene的可以转换的各种ID

Ensembl数据库在线网页工具biomart简单教程356

然后可以选择你上传的gene的各种序列

Ensembl数据库在线网页工具biomart简单教程358

可以选择的信息非常多,基本上可以想到的转换在这里都能做!!!

但是,始终没有脚本方便,只适合不太懂编程的菜鸟使用!

然后点击result即可,看到结果还可以导出成txt文档,点击右上角的GO即可

Ensembl数据库在线网页工具biomart简单教程458

 

30

云服务器试用报告

刚买了一个空的云服务器,所以就试用体验了一下!

一.硬盘容量,系统盘很小,就30个G,但是有一个1T的空盘,自己格式化安装并挂载即可。

云服务器试用报告53

 

二.内存状况,11G

云服务器试用报告67

三。Cup数量,12核

 

云服务器试用报告79

四.磁盘文件状况

云服务器试用报告90

 

五.开通其它账号

adduser

六.服务器其它信息

云服务器试用报告119

 

七.测试下载速度

云服务器试用报告130

 

八.磁盘分区管理

 

首先看到的是32.2GB的系统盘,是标示符是xvda盘,被分成了三个区

云服务器试用报告141

然后是一个1T的硬盘,需要分区然后挂载

云服务器试用报告199

我分区后,分割成两个盘,一个给各个用户,还有一个放公共数据

apt-get install nfs-common

mkfs -t ext4 /dev/xvdd

mkfs.ext4 /dev/xvdd1

mount -t ext4 /dev/xvdd1 /home/

mount -t ext4 /dev/xvdd2  /data

云服务器试用报告366

 

九.脚本环境

This is perl 5, version 18, subversion 2 (v5.18.2)

/usr/lib/python2.7/site.pyc matches /usr/lib/python2.7/site.py

十.库文件状况

apt-get install unzip

apt-get install make

apt-get install gcc

十一.软件安装状况

云服务器试用报告570

 

常用生物信息学软件都可以自己安装,并且可以使用。

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

 

28

自己动手写bowtie第三讲:序列查询。

查询需要根据前面建立的索引来做。

BWT算法第二讲532

这是一个比较复杂的过程,我也是看了bowtie的作者的ppt才慢慢弄懂的,感觉自己也不可能三言两语就说清楚,一般都是辅助图片,动画,再经过多方交流才能慢慢理解。

所以大家呢,就自己去看ppt,看懂那个查询算法。(ppt及代码在我的群里面有共享,欢迎大家加群交流)

这里我简单讲讲我的程序

首先读取索引文件,统计好A,C,G,T的总数

然后把查询序列从最后一个字符往前面回溯。

我创建了一个子函数,专门来处理回溯的问题

每次接受四个参数(左右两端的碱基,上下的阈值),并返回两个参数(新的上下两个阈值)

自己动手写bowtie之三序列查询261

大家要看懂阈值是如何更新迭代,这样动态的一个个回溯字符串,一个个迭代阈值。

直到四种临界情况的出现。

自己动手写bowtie之三序列查询477

第一是上下阈值已经相等了,但是我们还没有回溯完全,那就说明字符串只能查找后几个字符,前面还有字符是无法匹配的

第二种情况是上下阈值已经相等了,正巧我们也回溯到了最后一个字符串,那么我们就找到了精确匹配。

第三种情况是已经进行到了最后一个字符串,但是上下阈值还有差值,那么就找到了多个精确匹配点。

最后一种情况是各种非法字符。

然后我简单的测序了一下在病毒的5K基因组里面的精确匹配情况,好像效果还挺好的

自己动手写bowtie之三序列查询518

但是在酵母里面还有一个问题没有解决,就是取前二十个字符串排序的问题,不够精确,需要重新审视排序结果进行局部优化,可能是需要用堆排序发,具体我还得考虑一个星期,只能等下周上课再看看了,平时太忙了,基本没时间码代码。

这里贴上我的代码给大家看看,

[perl]

$a='CGCTATGTACTGGATGCGCTGGCAAACGAGCCTGCCGTAAG';

while(<>){

chomp;

@F=split;

$hash_count_atcg{$F[0]}++;

$hash{$.}=$_;

}

$all_a=$hash_count_atcg{'A'};

$all_c=$hash_count_atcg{'C'};

$all_g=$hash_count_atcg{'G'};

$all_t=$hash_count_atcg{'T'};

#print "$all_a\t$all_c\t$all_g\t$all_t\n";

$len_a=length $a;

$end_a=$len_a-1;

print "your query is $a\n";

print "and the length of your query is $len_a \n";

foreach (reverse (0..$end_a)){

$after=substr($a,$_,1);

$before=substr($a,$_-1,1);

#对第一个字符进行找阈值的时候,我们需要人为的定义起始点!

if($_ == $end_a){

if ($after eq 'A') {

$start=1;

$end=$all_a;

}

elsif ($after eq 'C') {

$start=$all_a+1;

$end=$all_a+$all_c;

}

elsif ($after eq 'G') {

$start=$all_a+$all_c+1;

$end=$all_a+$all_c+$all_g;

}

elsif ($after eq 'T'){

$start=$all_a+$all_c+$all_g+1;

$end=$all_a+$all_c+$all_g+$all_t;

}

else {print "error !!! we just need A T C G !!!\n";exit;}

}

#如果阈值已经无法继续分割,但是字符串还未查询完

if ($_  > 0 && $start == $end) {

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

$find_len=length $find_char;

#这里需要修改,但是不影响完全匹配了

print "we can just find the last $find_len char ,and it is $find_char \n";

exit;

}

#如果进行到了最后一个字符

if ($_ == 0) {

if ($start == $end) {

print "It is just one perfect match ! \n";

my @F_start=split/\s+/,$hash{$start};

print "The index is $F_start[1]\n";

exit;

}

else {

print "we find more than one perfect match!!!\n";

#print "$start\t$end\n";

foreach  ($start..$end) {

my @F_start=split/\s+/,$hash{$_};

print "One of the index is $F_start[1]\n";

}

exit;

}

}

($start,$end)=&find_level($after,$before,$start,$end);

}

sub find_level{

my($after,$before,$start,$end)=@_;

my @F_start=split/\s+/,$hash{$start};

my @F_end=split/\s+/,$hash{$end};

if ($before eq 'A') {

return ($F_start[2],$F_end[2]);

}

elsif ($before eq 'C') {

return ($all_a+$F_start[3],$all_a+$F_end[3]);

}

elsif ($before eq 'G') {

return ($all_a+$all_c+$F_start[4],$all_a+$all_c+$F_end[4]);

}

elsif ($before eq 'T') {

return ($all_a+$all_c+$all_g+$F_start[5],$all_a+$all_c+$all_g+$F_end[5]);

}

else {print "sorry , I can't find the right match!!!\n";}

}

#perl -alne '{next if />/;$all.=$_;}END{print substr($all,308,10)}'   lambda_virus.fa 

[/perl]

28

自己动手写bowtie第二讲:优化索引

BWT算法第二讲192

其中第一讲我提到了一个简单的索引产生方式,因为是课堂就半个小时想的,很多细节没有考虑到,对病毒那种几K大小的基因组来说是很简单的,速度也非常快,但是我测试了一下酵母,却发现好几个小时都没有结果,我只好kill掉重新改写算法,我发现之前的测序最大的问题在于没有立即substr函数的实现方式,把一个5M的字符串不停的截取首尾字符串好像是一个非常慢的方式。

BWT算法第二讲241

所以我优化了那个字符串的函数,虽然代码量变多了,实现过程也繁琐了一点,但是速度提升了几千倍。

 

time perl bwt_new_index.pl e-coli.fa >e-coli.index

测试了一下我的脚本,对酵母这样的5M的基因组,索引耗费时间是43秒

real 0m43.071s

user 0m41.277s

sys 0m1.779s

输出的index矩阵如下,我简单的截取头尾各10行给大家看,一点要看懂这个index。

BWT算法第二讲532

首先第一列就是我们的BWT向量,也就是BWT变换后的尾字符

第二列是之前的顺序被BWT变换后的首字符排序后的打乱的顺序。

第三,四,五,六列分别是A,C,G,T的计数,就是在当行之前累积出现的A,C,G,T的数量,是对第一列的统计。

 

这个索引文件将会用于下一步的查询,这里贴上我新的索引代码,查询见下一篇文章

[perl]

while (<>){

        next if />/;

        chomp;

        $a.=$_;

}

$len=length $a;

open FH_F,">tmp_forward.txt";

open FH_R,">tmp_reverse.txt";

for(my $i=0;$i<=$len-1;$i+=20){

        print FH_F  substr($a,$i,20);

        print FH_F "\n";

}

$rev_a=reverse $a;

for(my $i=0;$i<=$len-1;$i+=20){

        print FH_R  substr($rev_a,$i,20);

        print FH_R "\n";

}

close FH_F;

close FH_R;

$a='';

open FH_F,"tmp_forward.txt";

open FH_R,"tmp_reverse.txt";

#把前一行的所有20bp碱基当做后一行的头部信息

$residue_F=<FH_F>;

$residue_R=<FH_R>;

$i=0;

while  ($F_reads=<FH_F>){

        $R_reads=<FH_R>;

        $F_merge=$residue_F.$F_reads;

        $R_merge=$residue_R.$R_reads;

#这样每次就需要处理20个碱基

foreach  (0..19) {

$up  =substr($F_merge,$_,20);

                        $down=substr($R_merge,$_,1);

                        $hash{"$up\t$down"}=$i;

                        $i++;

}

#处理完毕之后再保存当行的20bp碱基做下一行的头部信息

$residue_F=$F_reads;

$residue_R=$R_reads;

}

#print "then we sort it\n";

$count_a=0;

$count_c=0;

$count_g=0;

$count_t=0;

foreach  (sort keys  %hash){

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

$len=length;

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

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

$count_a++ if $last eq 'A';

$count_c++ if $last eq 'C';

$count_g++ if $last eq 'G';

$count_t++ if $last eq 'T';

print "$last\t$hash{$_}\t$count_a\t$count_c\t$count_g\t$count_t\n";

}

unlink("tmp_forward.txt");

unlink("tmp_reverse.txt");

[/perl]