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

相关文章

来自专栏机器学习从入门到成神

校招面试知识点复习之计算机网络

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/sinat_35512245/articl...

611
来自专栏iOS 开发杂谈

iOS开发之iOS10适配

公司的项目,临上线之前做了一下iOS10的适配,发现一大堆的坑,瞬间觉得苹果不友好了。 一、证书问题 打开xcode8.0时编译运行时出现下面问题:

732
来自专栏葡萄城控件技术团队

ComponentOne使用技巧——从Winform穿越到WPF

WPF 和 Winform 是两个单独的平台,但二者又都是基于 .NET 4.0 以上版本开发的,所以很多.NET开发人员就开始研究如何在WPF中使用Winfo...

602
来自专栏FreeBuf

CVE-2017-8759完美复现(另附加hta+powershell弹框闪烁解决方案)

CVE-2017-8759 是前几天出的 0 DAY ,搜了下,国内几乎没有人复现,这个洞总体来说,危害很大,而且比CVE-2017-0199 更难防御。 漏...

24610
来自专栏YouMeek

1.5 Elasticsearch DSL 聚合语法介绍

课程环境 CentOS 7.3 x64 JDK 版本:1.8(最低要求),主推:JDK 1.8.0_121 Elasticsearch 版本:5.2.0 相关软...

3097
来自专栏JackeyGao的博客

大话西游答题器 command line

571
来自专栏hotqin888的专栏

电子规范管理系统(2)

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/hotqin888/article/det...

681
来自专栏菩提树下的杨过

javascript读写本机文本文件

近日在工作中遇到了一个情况:一张纯html的网页,用它一条一条输入数据,然后由JS运算出结果(这些数据多半都是临时的,所以也没考虑保存到数据库),每次用完后换台...

1847
来自专栏Python中文社区

进击的爬虫:用Python搭建匿名代理池

專 欄 ❈ 苍冥,Python中文社区专栏作者,澳洲华裔,目前在墨尔本某国际咨询公司任职Splunk Developer,擅长网络安全及攻防,热爱Python...

2695
来自专栏北京马哥教育

腾讯反病毒实验室为你揭秘Xshell软件后门技术!

背景: 最近,XShell远程终端工具发现被加入了恶意代码,目前官方就此事也做出了回应,要求使用者尽快下载最新版本。腾讯安全反病毒实验室就此跟进分析,对此次带...

2756

扫码关注云+社区