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

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

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

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

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

Kraken软件：

这个工具要求太高，一般人的电脑无法满足， 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比对软件：

如果仅仅是对微生物进行分类，硬件要求就低很多啦  

GOTTCHA软件：

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

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

github地址：

说明书非常详细。

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

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

1. cd ~/biosoft

2. # https://github.com/LANL-Bioinformatics/GOTTCHA

3. mkdir GOTTCHA && cd GOTTCHA

4. git clone https://github.com/LANL-Bioinformatics/GOTTCHA

5. cd GOTTCHA

6. wget ftp://ftp.lanl.gov/public/genome/gottcha/GOTTCHA_database_v20150825/GOTTCHA_lookup.tar.gz

7. wget ftp://ftp.lanl.gov/public/genome/gottcha/GOTTCHA_database_v20150825/GOTTCHA_BACTERIA_c4937_k24_u30_xHUMAN3x.species.tar.gz

8. tar -zxvf GOTTCHA_lookup.tar.gz

9. tar -zxvf GOTTCHA_BACTERIA_c4937_k24_u30_xHUMAN3x.species.tar.gz

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

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

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

1. ~/biosoft/GOTTCHA/GOTTCHA/bin/gottcha.pl \

2. --threads 5 \

3. --mode all \

4. --input ~/data/project/myGenome/bamFiles/P_jmzeng.unmapped.fq \

5. --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的代码似乎有问题

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

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

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

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

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