【直播】我的基因组74:快速给测序reads比对到物种

其实这一讲只是把未比对到人类基因组的序列快速比对到细菌基因组,并得到各个种类的菌的占比。

在这之前我们讲的是对几亿条reads定位到指定参考基因组的具体某个坐标,那是因为我们预先知道那些reads来自于人类,就是我本人血液的测序结果。

直播】我的基因组(十五):提取未比对的测序数据

但是前面也说到了那8.9亿reads里面是有部分(850万)无法比对上的,如果我们需要探究它们到底是什么东西,会不会是其它物种的DNA物质掺和进来了呢?

我们不可能把所有物种的参考基因组都下载一遍,然后一个个的去比对。这时候就需要一些特殊的数据库和比对工具了。

Kraken软件:

https://ccb.jhu.edu/software/kraken/MANUAL.html

这个工具要求太高,一般人的电脑无法满足, Construction of Kraken's standard database will require at least 160 GB of disk space.内存要求也很大,The default database size is 75 GB (as of Feb. 2015), and so you will need at least that much RAM

taxoner比对软件:

如果仅仅是对微生物进行分类,硬件要求就低很多啦 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4117525/

GOTTCHA软件:

http://nar.oxfordjournals.org/content/early/2015/03/12/nar.gkv180.full

但是数据库的要求还是很高,可以就用NCBI的NR库,其实做blast就好了,但是很多加速的软件发表了。但是我看了一下NR库,太可怕了,压缩包就24.7G 了,里面包含了所有生物的已知基因。但是我只是想看看我的一些reads的物种分布而已, ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/

考虑到这个需求其实并不紧急,我也懒得下载那么大的数据库来做比对,而且blast这样的比对非常之慢。所以我就试用其中最简单的GOTTCHA这个软件吧!

github地址:https://github.com/LANL-Bioinformatics/GOTTCHA

说明书非常详细。

当然,也是需要下载数据库,但是稍微小一点,如下:

软件及数据库下载及安装方法是:

cd ~/biosoft
# https://github.com/LANL-Bioinformatics/GOTTCHA
mkdir GOTTCHA && cd GOTTCHA
git clone https://github.com/LANL-Bioinformatics/GOTTCHA
cd GOTTCHA
wget ftp://ftp.lanl.gov/public/genome/gottcha/GOTTCHA_database_v20150825/GOTTCHA_lookup.tar.gz
wget ftp://ftp.lanl.gov/public/genome/gottcha/GOTTCHA_database_v20150825/GOTTCHA_BACTERIA_c4937_k24_u30_xHUMAN3x.species.tar.gz
tar -zxvf GOTTCHA_lookup.tar.gz
tar -zxvf GOTTCHA_BACTERIA_c4937_k24_u30_xHUMAN3x.species.tar.gz

首先把未匹配的reads做成fastq文件,用来比对

cat P_jmzeng.unmapped.sam |perl -alne '{ ++$h{$F[0]};print q{@}."$F[0]_$h{$F[0]}\n$F[9]\n+\n$F[10]"}' >P_jmzeng.unmapped.fq

有了输入文件,前面软件和数据库也准备好了,就可以直接运行啦:

~/biosoft/GOTTCHA/GOTTCHA/bin/gottcha.pl \
--threads 5 \
--mode all \
--input ~/data/project/myGenome/bamFiles/P_jmzeng.unmapped.fq \
--database database/GOTTCHA_BACTERIA_c4937_k24_u30_xHUMAN3x.species

这一个命令里面包含着3个步骤:

(1) split-trimming the input data

(2) mapping reads to a GOTTCHA database using BWA

(3) profiling/filtering the results.

看起来这个软件使用很简单,但事实上,我被虐了两个多小时,它里面的代码下载YAML的代码似乎有问题

perl -MYAML::XS -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext YAML::XS;
perl -MYAML -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext YAML;

导致它自己的脚本总是找不到YAML.pm ,所以我把系统自带的YAML文件夹拷贝到了ext文件夹。 伤心~~~

解决了YAML模块问题,就可以运行成功啦,logo日志如下:

我打开了结果文件看了看,先卖个关子,下次再讲结果解读。

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-05-03

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏大内老A

The .NET of Tomorrow

Ed Charbeneau(http://developer.telerik.com/featured/the-net-of-tomorrow/) Exciti...

31710
来自专栏张善友的专栏

Silverlight + Model-View-ViewModel (MVVM)

     早在2005年,John Gossman写了一篇关于Model-View-ViewModel模式的博文,这种模式被他所在的微软的项目组用来创建Expr...

2968
来自专栏pangguoming

Spring Boot集成JasperReports生成PDF文档

由于工作需要,要实现后端根据模板动态填充数据生成PDF文档,通过技术选型,使用Ireport5.6来设计模板,结合JasperReports5.6工具库来调用渲...

1.2K7
来自专栏Ceph对象存储方案

Luminous版本PG 分布调优

Luminous版本开始新增的balancer模块在PG分布优化方面效果非常明显,操作也非常简便,强烈推荐各位在集群上线之前进行这一操作,能够极大的提升整个集群...

3145
来自专栏闻道于事

js登录滑动验证,不滑动无法登陆

js的判断这里是根据滑块的位置进行判断,应该是用一个flag判断 <%@ page language="java" contentType="text/html...

6818
来自专栏Golang语言社区

【Golang语言社区】GO1.9 map并发安全测试

var m sync.Map //全局 func maintest() { // 第一个 YongHuomap := make(map[st...

4708
来自专栏张善友的专栏

Mix 10 上的asp.net mvc 2的相关Session

Beyond File | New Company: From Cheesy Sample to Social Platform Scott Hansel...

2577
来自专栏落花落雨不落叶

canvas画简单电路图

61811
来自专栏杨龙飞前端

scrollto 到指定位置

2504
来自专栏我和未来有约会

Silverlight第三方控件专题

这里我收集整理了目前网上silverlight第三方控件的专题,若果有所遗漏请告知我一下。 名称 简介 截图 telerik 商 RadC...

4025

扫码关注云+社区