【直播】我的基因组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 条评论
登录 后参与评论

相关文章

使用Bluemix,NoSQL DB和Watson创建云应用程序

大家好,因为近几年工作很忙,我已经很久没有写过文章了。我现在是IBM的Bluemix平台的云架构师。我曾经使用Tomcat服务器上的Web应用程序编写了一个在B...

1846
来自专栏老九学堂

有个程序猿很忧桑:一个命令rm -rf/ ,他把整个公司删没了...

▼ 话说 最近有个程序猿很忧桑 ....... ? 因为弄错了一行代码 这哥们不小心把他整个公司 删没了 没了 了 ... ? ? 好吧.. 事情是这样的.....

2615
来自专栏王大锤

iOS各种调试技巧豪华套餐

3519
来自专栏数据小魔方

是时候展现真正的技术了——让你的图表舞动起来~

最近更新的频率越来越少了,因为我个人原因,各种事情在忙,实在对不起大家,今天写一篇补偿一下~ 既然写的少了,那么以后每写一篇,都会让技术含量高一些,这样才能对得...

2937
来自专栏CSDN技术头条

分享11款主流的开源编程工具

导读:有了开源编程工具,在基于开源许可证的情况下您可以轻松学习、修改、提高代码的质量,本文收集了11款最主流的且有价值的开源编程工具。或许会给您带来一丝惊喜。一...

1796
来自专栏程序员的碎碎念

Hello World,Python!

早上在床上挣扎了一个多小时,才起床的我,思来想去,还是觉得学一些python可能会有用些。 于是乎,刷牙洗脸,背上书包,吃个早饭,跑到图书馆四...

3044
来自专栏Keegan小钢

App项目实战之路(四):UI篇

上一篇文章[原型篇]发布之后,就开始设计UI了,包括Icon和界面UI,周一到周五晚上一般花两到三到小时,周六日的时候则有五六个小时,最终用了一个星期多才设计完...

1003
来自专栏SeanCheney的专栏

Python工程的文档结构

Python工程的文档结构,可以参考https://stackoverflow.com/questions/193161/what-is-the-best-pr...

862
来自专栏C/C++基础

程序是什么

程序(Program)是计算机系统的必备元素,因为计算机系统由硬件、操作系统以及软件构成,而程序又是软件的组成部分。操作系统是管理和控制计算机硬件与软件资源的计...

772
来自专栏前端新视界

浅谈 Angular 项目实战

我不是 Angular 的布道者,但如今自称 Angular 派,使用 Angular 做项目让我有一种兴奋感。目前的三大主流前端框架都研究过,博客中也有三者的...

660

扫描关注云+社区