宏基因组：一日一工具之-kneaddata -完成质控加去宿主
一日一条命令-kneaddata
写在前面
kneaddata是一分结合质控和去除宿主的过程，集合了两个软件Trimmomatic和 Bowtie2，Trimmomatic作为质控软件应用的非常多，但是去除宿主的还有BMTagger软件。。这部分的核心是对序列和宿主的序列进行比对，并去除匹配上的序列，所以需要宿主基因组序列，不同的软件可能要构建不同的所以索引。值得注意的是这是Huttenhower lab开发专门质控宏基因组的工具。
KneadData需要FASTQ文件。它支持单端和双端序列。
安装要求
Trimmomatic （版本== 0.33）（自动安装）
Bowtie2 （版本> = 2.2）（自动安装）
Python （版本> = 2.7）
Java运行时环境
TRF （可选）
Fastqc （可选）
SAMTools （仅当输入文件为BAM格式时才需要）
内存（如果使用Bowtie2，> = 4 Gb，如果使用BMTagger，> = 8 Gb）
操作系统（Linux或Mac）
今日份—kneaddata
这是一个综合软件，用于指控的去除宿主污染，所以这个软件封装了 两个子软件。
第一个就是：Trimmomatic
Trimmomatic软件，2014年首次发表在Bioinformatics期刊上，是一款专门对Illumina平台测序产生的reads进行修剪和过滤的软件。自发表以来，Trimmomatic软件凭借其简单的安装方法、较快的运行速度（支持多线程）、强大的去接头能力（simple和palindrome两种模式）、多元化的低质量数据处理方式、人性化的输出格式（clean reads为一一对应的pair-end形式，无需再次处理）等特点，深受数据处理者喜爱！
安装过程
直接使用conda即可安装成功
# 测试 Ubuntu 20
conda install kneaddata
Bowtie2
Bowtie2 只能索引基因组文件。是将测序reads与长参考序列比对工具，注意着个工具虽然类似我们RNAseq所使用的hisat2，并且用法也类似，但是Bowtie2 只能索引基因组文件。
这个软件已经在安装kneaddata的时候作为依赖安装进去了。
但是却需要基因组索引文件，这里我们需要手动下载，或者使用命令
软件提供的这个地方只有人类和小鼠的基因组索引文件，和几个测试文件，基本上是不够我们做生态的使用的。
下载索引到当前文件夹中
kneaddata_database --download human_genome bowtie2 ./
可选的数据框有以下几个：可以通过kneaddata_database —available 命令查看。
human_genome : bmtagger = http://huttenhower.sph.harvard.edu/kneadData_databases/Homo_sapiens_BMTagger_v0.1.tar.gz
human_genome : bowtie2 = http://huttenhower.sph.harvard.edu/kneadData_databases/Homo_sapiens_Bowtie2_v0.1.tar.gz
mouse_C57BL : bowtie2 = http://huttenhower.sph.harvard.edu/kneadData_databases/mouse_C57BL_6NJ_Bowtie2_v0.1.tar.gz
human_transcriptome : bowtie2 = http://huttenhower.sph.harvard.edu/kneadData_databases/Homo_sapiens_hg38_transcriptome_Bowtie2_v0.1.tar.gz
ribosomal_RNA : bowtie2 = http://huttenhower.sph.harvard.edu/kneadData_databases/SILVA_128_LSUParc_SSUParc_ribosomal_RNA_v0.1.tar.gz
但是如果有其他物种，这里可以单独使用bowtie2构建索引，就可以用使用在kneaddata 中了。kneaddata也没有自己构建很多的基因组供大家下载，还是成本太高啦。
使用Bowtie2构建索引
这里我构建了玉米的索引，同样的方式适用于植物或者动物：
Ensembl数据库下载参考基因组（植物的基因组）：
http://plants.ensembl.org/index.html
动物的基因组：
http://asia.ensembl.org/index.html?redirect=no
注意下载toplevel的全基因组。
bowtie2-build ./Zea_mays.B73_RefGen_v4.dna.toplevel.fa index
参考学习：http://www.chenlianfu.com/?p=178
Trimmomatic
illumina的几个接头和trimmomatic
一般illumina就那么几个接头，如果你想知道接头的序列可以去illumina下载，或者再trimmomatic的安装文件里都有的。
点击下载
应该选择什么接头呢？
我们一般选择这个接头：trimmomatic-0.39-1//adapters/TruSeq3-PE.fa
为什么选择这个呢？因为这和TruSeq protocol 协议有关系，一般来讲V2 使用在一些比较老的测序机器上，例如：GAII，比较新的及其，例如：HiSeq 。
SE 或者 PE分别代表是单端测序还是双端。
所以 TruSeq-3-PE.fa 就是最近几年的及其测的双端序列。
为什么我们还需要了解其他格式的接头呢？因为在NABI等数据库中的宏基因组数据可能会用到。
reference
https://www.biostars.org/p/323087/
https://www.genepattern.org/modules/docs/Trimmomatic/
kneaddata命令
是否去除重复序列
注意参数 —run-trim-repetitive ，为去除重复序列，这对于大量的重复序列非常重要，所以我们在运行的过程中添加上，例如：扩增子序列中就有大量的重复序列，添加该参数有助于减少运行时间。
如果不知道是否添加这个参数，可以使用fastqc来检测，添加参数：—fastqc FastQC；
kneaddata --input demo.fastq -db demo_db -o kneaddata_output --run-trim-repetitive --fastqc FastQC
如何指定bowtie2和Trimmomatic工具自己的参数
使用—bowtie2-options=后面用双引号接参数；Trimmomatic工具也是一样的道理；—trimmomatic-options=后面用双引号接参数；如果同一个工具添加多个参数，那就多用几个这种格式: 例如：
kneaddata --input demo.fastq --output kneaddata_output --reference-db database_folder --bowtie2-options="--very-fast" --bowtie2-options="-p 2"
kneaddata --input demo.fastq --output kneaddata_output --reference-db database_folder --trimmomatic-options="LEADING:3" --trimmomatic-options="TRAILING:3"
取宿主可可以去除多个宿主
这个软件专门设计用于宏基因组，对于处理宿主数据就显得很得心应手，这里可以添加多个宿主，用于复杂环境去除宿主数据，例如：
kneaddata --input demo.fastq --output kneaddata_output --reference-db database_folder --reference-db database_folder2/demo
根据示例文件运行数据
time kneaddata -i seq/p144C_1.fq.gz -i seq/p144C_2.fq.gz -o temp/qc -v -t 3 --remove-intermediate-output --trimmomatic /home/wentao/sra/Trimmomatic-0.38/ --trimmomatic-options 'SLIDINGWINDOW:4:20 MINLEN:50' --bowtie2-options '--very-sensitive --dovetail' -db ~/db/kneaddata/human_genome/Homo_sapiens
进行完成后的序列统计
kneaddata_read_count_table --input temp/qc --output result/01kneaddata_sum.txt
cat result/01kneaddata_sum.txt
结果文件
$OUTPUT_DIR/$INPUTNAME_kneaddata.fastq ：质控后并过滤完成的最终文件
$OUTPUT_DIR/$INPUTNAME_kneaddata_$DATABASE_$SOFTWARE_contam.fastq 宿主比对上的序列
$OUTPUT_DIR/$INPUTNAME_kneaddata.log 运行的日志文件，如果报错可以查看这里。
$OUTPUT_DIR/$INPUTNAME_kneaddata.trimmed.fastq 这仅仅是Trimmomatic过滤后的文件
reference : https://huttenhower.sph.harvard.edu/kneaddata/
reference
https://www.sohu.com/a/252244330_100269283
https://github.com/biobakery/kneaddata
常见问题- 可以不做阅读
pip安装
此安装方式会同时安装Trimmomatic和Bowtie2，以为好多人都已经安装过这两个依赖了，所以可以在安装过程中添加参数
pip install kneaddata
pip install kneaddata –install-option=’–bypass-dependencies-install
手动安装过程 Trimmomatic
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.38.zip
unzip Trimmomatic-0.38.zip
java -jar ~/sra/Trimmomatic-0.38/trimmomatic-0.38.jar –h
创建BMTagger数据库
默认情况下，这将生成参考数据库，其名称以开头reference.fasta。
kneaddata_build_database reference.fasta
# bowtie2 数据库或者索引构建
bowtie2-build Homo_sapiens.fasta -o Homo_sapiens_db
Homo_sapiens_db让多个索引文件具有统一的扩展名。
关于adapters的问题
按照刘老师说的，目前公司都去除了接头了，那么如果没有去除接头，我们该如何选择adapters序列文件呢？
进入trimmomatic路径我们发现有很多的adapters。
cd ./metagenome_env/share/trimmomatic-0.39-1/adapters/
这里有六个，那么如何决定选择哪个adapters来做trimmomatic，去除接头呢呢？
NexteraPE-PE.fa TruSeq2-SE.fa TruSeq3-PE.fa
TruSeq2-PE.fa TruSeq3-PE-2.fa TruSeq3-SE.fa
java错误: -d64
报错信息：
Error message returned from Trimmomatic :
Unrecognized option: -d64
Error: Could not create the Java Virtual Machine.
Error: A fatal exception has occurred. Program will exit.
我不知道这是什么问题：首先开始修改java的版本 查找ubuntu下面到底安装了多少个java的版本：
sudo update-alternatives --config java
发现我安装了好几个java。
0 /usr/lib/jvm/java-11-openjdk-amd64/bin/java 1111 auto mode
1 /usr/lib/jvm/java-11-openjdk-amd64/bin/java 1111 manual mode
* 2 /usr/lib/jvm/java-8-openjdk-amd64/jre/bin/java 1081 manual mode
3 /usr/local/java/jdk1.7.0_80/bin/java
这种方法选择数字没有更改了默认版本，所以我在这个文件中修改：
sudo gedit /etc/profile
在文件中修改JAVA_HOME即可（输入sudo update-alternatives —config java即可看到这些路径）
export JAVA_HOME=/usr/local/java/jdk1.7.0_80
export JRE_HOME=${JAVA_HOME}/jre
export CLASSPATH=.:$JAVA_HOME/lib:$JRE_HOME/lib:$CLASSPATH
export PATH=$JAVA_HOME/bin:$JRE_HOME/bin:$PATH
完成之后source一下
source /etc/profile