宏转录组学习笔记--另一个教程

用户1075469

发布于 2020-03-31 15:10:38

2.6K0

发布于 2020-03-31 15:10:38

这项工作已获得Creative Commons Attribution-ShareAlike 4.0 International协议的许可。这意味着您可以复制，共享和修改作品，只要结果以相同的许可证分发即可。本教程由Mobolaji Adeolu（adeolum@mcmaster.ca），John Parkinson（john.parkinson@utoronto.ca）和Xuejian Xiong（xuejian@sickkids.ca）制作。

「注意，这个教程的软件运行环境为linux，没有相关环境需要使用docker或者虚拟机，而且，经过测试，python版本要求为2.7， biopython=1.67，在不停报错的教训中得到的结论。」

总览

本教程将带您完成处理元转录组数据的流程。实验室开发的流程包括以下各个步骤：

去除在文库制备和测序步骤中添加的接头序列，并修剪低质量的碱基和测序reads。
删除重复的reads以减少以下步骤的处理时间。
除去载体污染（来自克隆载体，刺突和引物的reads）。
除去宿主reads（如果要研究其中存在宿主的微生物组）。
尽管使用了rRNA去除试剂盒，但仍要删除通常主导转录组数据集的大量rRNA序列。
将重复的reads（在步骤2中删除）添加回数据集，以提高组装的质量。
将reads分类到已知的物种分类组，并可视化数据集的分类组成。
将reads组装成重叠群中以提高注释质量。
注释reads已知基因。
将已鉴定的基因映射到swiss-prot数据库中以鉴定酶功能
生成与每个基因相关的标准化表达值。
使用KEGG代谢途径作为Cytoscape的重叠群，可视化结果。

整个宏转录组学流程包括现有的生物信息学工具和一系列处理文件格式转换和输出解析的Python脚本。我们将通过以下步骤来说明流程的复杂性以及基础工具和脚本。

新的，更快和/或更准确的工具一直在开发，值得牢记的是，随着这些流程被社区采纳为标准，任何流程都需要灵活地整合这些工具。例如，在过去的两年中，我们的实验室已经从cross_match过渡到Trimmomatic，从BLAST过渡到DIAMOND。

「译者注：本研讨会旨在与DIAMOND v0.8.26一起使用。较新版本的DIAMOND将与我们在此练习中制作的预编译数据库文件不兼容」。

为了说明该过程，我们将使用从小鼠结肠内容产生的序列reads。这些是150 bp单端reads。也可以使用成对末端的reads，并且通常是首选的，因为当reads对中有足够的重叠以提高有效平均reads长度时，它们可以提高注释质量。使用成对末端数据需要一个额外的数据处理步骤（合并重叠的reads），从而在数据处理过程中生成更多文件（用于合并/单reads，正向reads和反向reads的文件），但是成对末端reads的结构数据类似于此处描述的reads，并且可以轻松进行调整。

本教程将带您逐步处理100000个reads的一部分，而不是使用整个2500万个reads的整个过程（后者在桌面上可能要花费几天的时间）。

开场

工作目录

创建一个新目录，该目录将存储在本实验中创建的所有文件。

mkdir -p ~/metatranscriptomics
cd ~/metatranscriptomics

Python脚本

我们已经编写了许多脚本来从您将要使用的工具中提取和分析数据。下载用于宏转录组学研讨会的软件包，并解压缩我们的python脚本。

#原文采用wget，这里为了加速，采用多线程的axel
axel https://github.com/ParkinsonLab/2017-Microbiome-Workshop/releases/download/Extra/precomputed_files.tar.gz
wget https://jiawen.zd200572.com/hla/precomputed_files.tar.gz --no-check-certificate
wget https://github.com/ParkinsonLab/Metatranscriptome-Workshop/archive/EC.zip --no-check-certificate
unzip EC.zip
tar -zxvf precomputed_files.tar.gz *.py  --wildcards

输入文件

我们的数据集包含从小鼠结肠内容物产生的150 bp单端Illumina reads。检查其内容：

tar -xvf precomputed_files.tar.gz mouse1.fastq
less mouse1.fastq

注意事项：

输入q以退出less。

使用FastQC检查reads质量

source ~/Miniconda/bin/activate
conda create -n metatranscripts python=2.7
conda activate metatranscripts
conda install -y bwa SPAdes fastqc biopython=1.67 diamond=0.826 blat blast samtools bwa
fastqc mouse1.fastq -t 4 #只有一个文件，所以4线程是可选的

FastQC报告在HTML文件中生成mouse1_fastqc.html。您还将找到一个zip文件，其中包含用于生成报告的数据文件。

要打开HTML报告文件，请使用浏览器浏览mouse1_fastqc.html并查找以下信息：

基本统计信息：小鼠RNA序列数据的基本信息，例如reads总数，reads长度，GC含量。
每碱基序列质量：每个位置上所有碱基的质量值范围的概述。
每碱基序列含量：显示跨序列长度的核苷酸偏差的图。
适配器内容：提供有关序列样品中适配器污染程度的信息。

处理reads

步骤1.移除接头序列并修剪低质量序列。

Trimmomatic可以快速识别和修剪接头序列，以及识别和删除低质量序列数据

ln -s ~/Miniconda/envs/tera/share/trimmomatic-*/adapters/TruSeq3-SE.fa Adapters
#ln -s /usr/local/prg/Trimmomatic-0.36/adapters/TruSeq3-SE.fa Adapters
Trimmomatic SE mouse1.fastq mouse1_trim.fastq ILLUMINACLIP:Adapters:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:50
#运行过程的输出
Automatically using 4 threads
Using Long Clipping Sequence: 'AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA'
Using Long Clipping Sequence: 'AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC'
ILLUMINACLIP: Using 0 prefix pairs, 2 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences
Quality encoding detected as phred33
Input Reads: 100000 Surviving: 94415 (94.42%) Dropped: 5585 (5.58%)
TrimmomaticSE: Completed successfully

ln -s /usr/local/prg/Trimmomatic-0.36/adapters/TruSeq3-SE.fa Adapters用于创建指向Trimmomatic提供的单端接头序列文件的符号链接，该文件适用于HiSeq和MiSeq机器生成的序列。但是，如果可能，应使用您自己的测序项目中的已知适配器文件替换此文件。

「命令行参数是：」

SE：输入数据是单端reads。
ILLUMINACLIP:Adapters:2:30:10：卸下适配器。
LEADING:3：如果reads的质量得分低于3，则在reads开始时将其作为基础。
TRAILING:3：如果它们的质量得分低于3，则在reads结束时修剪基数。
SLIDINGWINDOW:4:15：使用大小为4的窗口进行扫描，以reads局部质量低于15的reads，如果发现则进行修剪。
MINLEN:50：删除长度小于50的序列。

「问题1：删除了多少个低质量序列？」(5585)

使用FastQC检查reads质量：fastqc mouse1_trim.fastq -t 4与上一份报告进行比较，以查看以下各节中的更改：

基本统计
每碱基序列质量

可选：配对读合并

如果您使用的是配对末端数据集，我们可以识别重叠的序列读对，因此可以合并为单个序列。为此，我们使用工具VSEARCH，可以在以下网站上找到该工具：

Exmaple only, do not run!
vsearch --fastq_mergepairs mouse1_trim.fastq --reverse mouse2_trim.fastq --fastqout mouse_merged_trim.fastq --fastqout_notmerged_fwd mouse1_merged_trim.fastq --fastqout_notmerged_rev mouse2_merged_trim.fastq

命令行参数是：

--fastq_mergepairs 指示VSEARCH使用reads合并算法来合并重叠的成对末端reads
--reverse 指示具有3'至5'（反向）配对末端reads的文件名
--fastqout 指示输出文件包含重叠的配对末端reads
--fastqout_notmerged_fwd和--fastqout_notmerged_rev指示输出文件包含不重叠的成对末端reads
如果要查看合并的reads长度的分布，可以使用fastqc检查reads属性：

fastqc mouse_merged_trim.fastq
mouse_merged_trim_fastqc.html

reads质量过滤

Trimmomatic用于删除reads中的接头和修剪低质量碱基，它使用滑动窗口方法删除reads中低质量碱基的连续区域。但是，值得加一个总体reads质量阈值，以确保在我们的分析中使用的所有reads均具有足够的正确率。为此，我们使用VSEARCH （在处理配对末端数据时，此步骤应在reads合并步骤之后执行）：

vsearch --fastq_filter mouse1_trim.fastq --fastq_maxee 2.0 --fastqout mouse1_qual.fastq

「命令行参数是：」

--fastq_filter 指示VSEARCH使用质量过滤算法删除低质量reads
--fastq_maxee 2.0预期的错误阈值。设置为1。任何质量得分表明平均预期错误数量大于1的reads都将被过滤。
--fastqout 指示输出文件包含高质量的过滤reads

使用FastQC检查reads质量：

fastqc mouse1_qual.fastqmouse1_qual_fastqc.html与以前的报告进行比较，以查看以下各节中的更改：

基本统计、每碱基序列质量、每序列质量

**问题2：每次reads序列质量曲线如何变化？**变好了一些

步骤2.删除重复的reads

为了大大减少识别和过滤rRNAreads所需的计算时间，我们使用CD-HIT执行去重复步骤，以删除重复的reads。

tar zxvf cd-hit-v4.8.1-2019-0228.tar.gz
cd cd-hit-v4.8.1-2019-0228/
make openmp=no #我这里黑苹果不支持这个命令，一般应该可以不用后边的参数
cd cd-hit-auxtools
make
cd-hit-v4.8.1-2019-0228/cd-hit-auxtools/cd-hit-dup -i mouse1_qual.fastq -o mouse1_unique.fastq

「命令行参数是：」

-i：输入的fasta或fastq文件。
-o：包含去重复序列的输出文件，其中唯一的代表序列用于表示具有多个重复的每组序列。
mouse1_unique.fastq.clstr创建第二个输出文件，该文件确切显示由去复制的文件中的每个唯一序列表示的复制序列，mouse1_unique.fastq2.clstr还创建了第三个空的输出文件，该文件仅用于配对末端reads。

「问题3：您能找到多少个独特的reads内容吗？」

cat mouse1_unique.fastq | wc -l # 341616

尽管在这个小型数据集中复制的reads次数相对较少，但对于较大的数据集，此步骤可以将文件大小减少多达50-80％

步骤3.去除载体污染

为了识别和过滤来自载体，接头，接头和引物污染源的reads，我们使用了Burrows Wheeler序列比对器（BWA）和BLAST类似的比对工具（BLAT）来搜索这些序列数据库。作为用于识别污染性载体和接头序列的参考数据库，我们依赖于UniVec_Core数据集，该数据集是从NCBI Univec数据库中已知载体以及的常见测序接头，接头和PCR引物的fasta文件。请首先将其下载到您的工作目录中。

axel ftp://ftp.ncbi.nih.gov/pub/UniVec/UniVec_Core现在，我们必须使用以下命令为BWA和BLAT的这些序列生成索引：

bwa index -a bwtsw UniVec_Core
samtools faidx UniVec_Core
makeblastdb -in UniVec_Core -dbtype nucl

接下来，我们可以使用BWA对reads进行比对，并使用以下命令使用Samtools筛选出与载体数据库比对的所有reads：

接下来，我们可以使用BWA对reads进行比对，并使用以下命令使用Samtools筛选出与数据库比对的所有reads：

bwa mem -t 4 UniVec_Core mouse1_unique.fastq > mouse1_univec_bwa.sam
samtools view -bS mouse1_univec_bwa.sam > mouse1_univec_bwa.bam
samtools fastq -n -F 4 -0 mouse1_univec_bwa_contaminats.fastq mouse1_univec_bwa.bam
samtools fastq -n -f 4 -0 mouse1_univec_bwa.fastq mouse1_univec_bwa.bam

用于执行以下任务的命令：

bwa mem：生成与载体污染物数据库的reads比对
samtools view：将bwa的.sam输出转换为.bam，以进行以下步骤
samtools fastq：生成所有的fastq输出reads映射到污染物数据库（-F 4）和所有reads没有映射到向量污染物数据库（-f 4）

「问题4：您能否找到映射到载体数据库的BWAreads数目？」

cat mouse1_univec_bwa_contaminats.fastq | wc -l #304

现在，我们想对使用BLAT的reads执行其他比对，以滤除与载体污染数据库比对的所有剩余reads。但是，BLAT仅接受fasta文件，因此我们必须将reads内容从fastq转换为fasta。可以使用VSEARCH完成。

vsearch --fastq_filter mouse1_univec_bwa.fastq --fastaout mouse1_univec_bwa.fasta注意事项：

所使用的VSEARCH命令--fastq_filter与步骤1中用于过滤低质量reads的命令相同。但是，这里我们没有提供过滤条件，因此所有输入reads都传递到输出fasta文件。现在，我们可以使用BLAT对载体污染数据库进行额外的比对。

blat -noHead -minIdentity=90 -minScore=65 UniVec_Core mouse1_univec_bwa.fasta -fine -q=rna -t=dna -out=blast8 mouse1_univec.blatout注意事项：

命令行参数是：

-noHead：禁止.psl标头（因此它只是一个制表符分隔的文件）。
-minIdentity：设置最小序列同一性为90％。
-minScore：设置最低分数为65。这是匹配减去不匹配减去某种空位罚分。
-fine：用于高质量的mRNA。
-q：查询类型为RNA序列。
-t：数据库类型为DNA序列。

最后，我们可以运行一个小的python脚本来过滤BLAT不能可靠地与我们的载体污染数据库中的任何序列比对的reads。

./1_BLAT_Filter.py mouse1_univec_bwa.fastq mouse1_univec.blatout mouse1_univec_blat.fastq mouse1_univec_blat_contaminats.fastq

注意事项：

该脚本的参数结构为：1_BLAT_Filter.py <Input_Reads.fq> <BLAT_Output_File> <Unmapped_Reads_Output> <Mapped_Reads_Output>

在这里，BLAT不会识别与载体污染物数据库比对的任何其他序列。但是，我们发现BLAT通常能够找到BWA无法识别的比对，特别是在搜索由全基因组组成的数据库时。

在数百万个大型reads数据集中对BWA遗漏的污染进行了一些比对。

步骤4.删除宿主reads

为了识别和过滤宿主reads（这里是小鼠来源的reads），我们使用小鼠DNA序列数据库重复上述步骤。为了我们的目的，我们使用从Ensembl下载的小鼠基因组数据库。

wget ftp://ftp.ensembl.org/pub/current_fasta/mus_musculus/cds/Mus_musculus.GRCm38.cds.all.fa.gz
gzip -d Mus_musculus.GRCm38.cds.all.fa.gz
mv Mus_musculus.GRCm38.cds.all.fa mouse_cds.fa

然后，我们重复上述步骤，为BWA和BLAT的这些序列生成索引：

bwa index -a bwtsw mouse_cds.fa
samtools faidx mouse_cds.fa
makeblastdb -in mouse_cds.fa -dbtype nucl

现在，我们使用BWA和Samtools比对并过滤出与我们的宿主序列数据库比对的所有reads：

bwa mem -t 4 mouse_cds.fa mouse1_univec_blat.fastq > mouse1_mouse_bwa.sam
samtools view -bS mouse1_mouse_bwa.sam > mouse1_mouse_bwa.bam
samtools fastq -n -F 4 -0 mouse1_mouse_bwa_contaminats.fastq mouse1_mouse_bwa.bam
samtools fastq -n -f 4 -0 mouse1_mouse_bwa.fastq mouse1_mouse_bwa.bam

最后，我们使用BLAT对我们的宿主序列数据库进行额外的比对。

vsearch --fastq_filter mouse1_mouse_bwa.fastq --fastaout mouse1_mouse_bwa.fasta
blat -noHead -minIdentity=90 -minScore=65  mouse_cds.fa mouse1_mouse_bwa.fasta -fine -q=rna -t=dna -out=blast8 mouse1_mouse.blatout
./1_BLAT_Filter.py mouse1_mouse_bwa.fastq mouse1_mouse.blatout mouse1_mouse_blat.fastq mouse1_mouse_blat_contaminats.fastq

「问题5：BWA和BLAT与小鼠宿主序列数据库比对了多少次？」

cat mouse1_mouse_blat_contaminats.fastq | wc -l #24

可选：在您自己的未来分析中，您可以选择使用载体污染数据库和宿主序列数据库的组合来同时完成步骤3和4 cat UniVec_Core mouse_cds.fa > contaminants.fa。但是，一起执行这些步骤使您很难分辨您的reads中有多少是专门来自宿主生物的。

步骤5.删除大量rRNA序列

rRNA基因倾向于在所有样品中高度表达，因此必须进行筛选，以避免组装和注释步骤的下游下游处理时间过长。您可以在此步骤中使用序列相似性工具（例如BWA或BLAST），但是我们发现[Infernal]（http://infernal.janelia.org/）速度较慢，但它更敏感，因为它依赖于协方差模型数据库（ CMs）描述基于Rfam数据库的rRNA序列图。由于对CMs的依赖，Infernal在单个内核上进行约100,000次reads最多可能需要4个小时。因此，我们将跳过此步骤，并使用mouse1_rRNA.infernalouttar文件中的预计算文件precomputed_files.tar.gz。

tar -xzf precomputed_files.tar.gz mouse1_rRNA.infernalout注意事项：

下面给出了用于生成此输出的命令：

vsearch --fastq_filter mouse1_mouse_blat.fastq --fastaout mouse1_mouse_blat.fasta
cmsearch -o mouse1_rRNA.log --tblout mouse1_rRNA.infernalout --anytrunc --rfam -E 0.001 Rfam.cm mouse1_mouse_blat.fasta

命令行参数是：

-o：输出日志文件。
--tblout：简单的表格输出文件。
--noali：从主输出中省略比对部分。这样可以大大减少输出量。
--anytrunc：放宽截断比对的阈值
--rfam：使用针对大型数据库设计的严格过滤策略。这将加快搜索速度，但可能会降低灵敏度。
-E：报告E值为0.001的靶序列。从此输出文件中，我们需要使用脚本来过滤rRNAreads：

从此输出文件中，我们需要使用脚本来过滤rRNAreads：

/Users/zd200572/Miniconda/envs/py2/bin/python ./2_Infernal_Filter.py mouse1_mouse_blat.fastq mouse1_rRNA.infernalout mouse1_unique_mRNA.fastq mouse1_unique_rRNA.fastq注意事项：

该脚本的参数结构为：2_Infernal_Filter.py <Input_Reads.fq> <Infernal_Output_File> <mRNA_Reads_Output> <rRNA_Reads_Output>

在这里，我们只删除了数千个reads，而不是映射到rRNA，但是在某些数据集中，rRNA最多可以代表80％的测序reads。

「问题6：鉴定了多少rRNA序列？现在还剩下多少？」

2253 reads were aligned to the rRNA database
81700 reads were not aligned to the rRNA database

步骤6.复制

去除污染物，宿主序列和rRNA后，我们需要将以前去除的重复reads替换回我们的数据集中。

./3_Reduplicate.py mouse1_qual.fastq mouse1_unique_mRNA.fastq mouse1_unique.fastq.clstr mouse1_mRNA.fastq注意事项：

该脚本的参数结构为：3_Reduplicate.py <Duplicated_Reference_File> <Deduplicated_File> <CDHIT_Cluster_File> <Reduplicated_Output>

问题7：确定了多少假定的mRNA序列？有多少个独特的mRNA序列？

81700 dereplicated reads produce 82623 rereplicated reads.

既然我们已经过滤了载体，接头，接头，引物，宿主序列和rRNA，请使用FastQC检查reads质量：

fastqc mouse1_mRNA.fastq -t 4mouse1_mRNA_fastqc.html「问题8：滤出了多少总污染物，宿主和rRNAreads？」

100000-82623=17377

步骤7.分类

现在我们有了推定的mRNA转录本，我们可以开始推断我们的mRNAreads的来源了。首先，我们将尝试使用基于引用的短阅读分类器来推断我们阅读的分类起点。在这里，我们将使用[Kaiju]（https://github.com/bioinformatics-centre/kaiju）基于参考数据库为我们的reads生成分类学分类。Kaiju可以在内存少于16GB（约13GB）的系统上使用proGenomes数据库，以每分钟数百万次reads的速度对原核生物reads进行分类。使用整个NCBI nr数据库作为参考大约需要43GB。同样，快速分类工具需要大于100GB的RAM才能对大型数据库的reads进行分类。但是，Kaiju对于我们的系统仍然占用了过多的内存，因此我们已经预先编译了分类，mouse1_classification.tsv，在tar文件中precomputed_files.tar.gz。

tar --wildcards -xzf precomputed_files.tar.gz kaiju*
chmod +x kaiju*
tar -xzf precomputed_files.tar.gz mouse1_classification.tsv nodes.dmp names.dmp

注意事项：

您使用的kaiju命令如下(ram硬件限制，不需要运行，kaiju_db.fmi文件也没提供)：./kaiju -t nodes.dmp -f kaiju_db.fmi -i mouse1_mRNA.fastq -z 4 -o mouse1_classification.tsv命令行参数是：

-t：分类ID的层次表示
-f：kaiju的预先计算的索引 -i：输入内容为
-z：系统上支持的线程数
-o：kaiju分类标准的输出文件

然后，我们可以进行分类阅读并进行补充分析。首先，我们将分类的特异性限制在属级分类单元上，这限制了虚假分类的数量。``/Users/zd200572/Miniconda/envs/py2/bin/python ./4_Constrain_Classification.py genus mouse1_classification.tsv nodes.dmp names.dmp mouse1_classification_genus.tsv` 注意事项：

该脚本的参数结构为：4_Constrain_Classification.py <Minimum_Taxonomic_Rank> <kaiju_Classification> <nodes_file> <names_file> <Output_Classifications>

然后，我们使用Kaiju生成人类可读的分类摘要。

./kaijuReport -t nodes.dmp -n names.dmp -i mouse1_classification_genus.tsv -o mouse1_classification_Summary.txt -r genus注意事项：

命令行参数是：

-t：分类ID的层次表示
-n：与每个分类ID对应的分类名称
-i：海归类分类
-o：摘要报告输出文件
-r：将为其生成摘要的分类等级

问题9：kaiju分类了多少reads？

20141到属,好少

0.062936	52	Viruses
33.335352	27543	cannot be assigned to a genus
42.225018	34888	unclassified

最后，我们将使用Krona生成数据集分类组成的分层多层饼图摘要。

./kaiju2krona -t nodes.dmp -n names.dmp -i mouse1_classification_genus.tsv -o mouse1_classification_Krona.txt
tar -xzf precomputed_files.tar.gz KronaTools
sudo KronaTools/install.pl
KronaTools/scripts/ImportText.pl -o mouse1_classification.html mouse1_classification_Krona.txt

然后，我们可以使用网络浏览器查看此数据集的饼图表示形式：mouse1_classification.html看这张交互式的图还是很漂亮的，展开后相当惊艳。

「问题10：在我们的数据集中，最丰富的科是什么？什么是最丰富的门？」提示：尝试减小Max depth屏幕左上方的值，和/或双击特定的分类单元。

Lachnospiraceae科和厚壁门

步骤8.组装reads

以前的研究表明，将reads组装成较大的重叠群，会大大提高我们通过序列相似性搜索将它们注释到已知基因上的能力。在这里，我们将SPAdes基因组组装者的转录本组装算法应用于我们推定的mRNAreads集。

spades.py --rna -s mouse1_mRNA.fastq -o mouse1_spades
mv mouse1_spades/transcripts.fasta mouse1_contigs.fasta

注意事项：

命令行参数是：

--rna：使用mRNA转录组装算法
-s：单端输入reads
-o：输出目录
SPAdes将reads组合成重叠群，这些重叠群被放置在名为的文件中 mouse1_spades/transcripts.fasta

「问题11：SPAdes生产了多少条组装序列？」(1082) 提示：尝试使用命令tail mouse1_contigs.fasta为了提取未组装的reads，我们需要通过BWA将所有推定的mRNAreads映射到我们的组装重叠群中。

首先，我们需要建立索引以允许BWA针对我们的重叠群进行搜索：

bwa index -a bwtsw mouse1_contigs.fasta接下来，我们尝试将整个假定的mRNA读图映射到此contig数据库：

bwa mem -t 4 mouse1_contigs.fasta mouse1_mRNA.fastq > mouse1_contigs.sam然后，我们将未映射的reads提取到fastq格式文件中以进行后续处理，并生成一个映射表，其中每个重叠群都与用于组装该重叠群的reads数目相关联。该表对于确定映射到重叠群的reads数目很有用，并用于确定相对表达（请参阅步骤6和8）。

/Users/zd200572/Miniconda/envs/py2/bin/python ./5_Contig_Map.py mouse1_mRNA.fastq mouse1_contigs.sam mouse1_unassembled.fastq mouse1_contigs_map.tsv注意事项：

该脚本的参数结构为：5_Contig_Map.py <Reads_Used_In_Alignment> <Output_SAM_From_BWA> <Output_File_For_Unassembed_Reads> <Output_File_For_Contig_Map>

「问题12：重叠群装配中没有使用多少个reads？重叠群装配中使用了多少个reads？」

22814 reads were aligned to the assemblies
59809 reads were not aligned to the assemblies

步骤9.注释reads已知的基因/蛋白质

在这里，我们将尝试推断推定的mRNAreads所源自的特定基因。在我们的reads中，我们依赖于精度递减的分层序列相似性搜索集-BWA和DIAMOND。虽然BWA提供高严格性，但在核苷酸水平上发生的序列多样性导致在这些过程中观察到的匹配很少。尽管如此，它还是很快的。为了避免在核苷酸水平上发生多样性的问题，尤其是在没有参考微生物基因组的情况下，我们使用DIAMOND搜索来提供更敏感的基于肽的搜索，这种搜索不太容易出现菌株之间的序列变化。

由于BWA利用核苷酸搜索，因此我们依赖于从微生物基因组数据库获得的数据NCBI包含5231个ffn文件。然后，我们将所有5231个ffn文件合并为一个fasta文件，microbial_all_cds.fasta并为此数据库建立索引以允许通过BWA搜索。对于DIAMOND搜索，我们使用NCBI的非冗余（NR）蛋白质数据库。

「注意事项：」

研讨会中使用的系统没有足够的内存来处理索引或搜索大型数据库，例如microbial_all_cds.fasta（9GB）和nr（> 60GB）。本节中的描述仅供参考。请使用我们预先计算的基因，蛋白质，并从tar文件中reads映射文件tar -xzf precomputed_files.tar.gz mouse1_genes_map.tsv mouse1_genes.fasta mouse1_proteins.fasta虽然我们在这里仅使用BWA，但可以使用BWA跟BLAT进行更彻底的搜索，microbial_all_cds.fasta如步骤3和4中所述。这限制了需要针对更大的NCBI nr数据库完成的搜索次数钻石/高炉。用于建立索引数据库的命令如下（您无需执行这些操作！）

bwa index -a bwtsw microbial_all_cds.fasta
samtools faidx microbial_all_cds.fasta
diamond makedb -p 8 --in nr -d nr

BWA搜索微生物基因组数据库

如果您要自己运行BWA，则可以使用以下命令搜索microbial_all_cds.fasta数据库：

bwa mem -t 4 microbial_all_cds.fasta mouse1_contigs.fasta > mouse1_contigs_annotation_bwa.sam
bwa mem -t 4 microbial_all_cds.fasta mouse1_unassembled.fasta > mouse1_unassembled_annotation_bwa.sam

然后，您将运行以下python脚本以提取与microbial_all_cds.fasta数据库的高可信度比对并生成reads基因映射表。在这里，我们每个重叠群只取一个基因，但是重叠群可能有多个基因（例如，共转录的基因）。

./6_BWA_Gene_Map.py microbial_all_cds.fasta mouse1_contigs_map.tsv mouse1_genes_map.tsv mouse1_genes.fasta mouse1_contigs.fasta mouse1_contigs_annotation_bwa.sam mouse1_contigs_unmapped.fasta mouse1_unassembled.fastq mouse1_unassembled_annotation_bwa.sam mouse1_unassembled_unmapped.fasta

该脚本的参数结构为：

6_BWA_Gene_Map.py <Gene_database> <Contig_Map> <Output_File_For_Gene_Map> <Output_File_For_Gene_sequences> <Contigs_File> <Contig_BWA_SAM> <Unmapped_Contigs> <Unassembled_Reads_File> <Unassembled_Reads_BWA_SAM> <Unmapped_Unassembled_Reads>

DIAMOND针对非冗余（NR）蛋白数据库

DIAMOND是一种BLAST样的本地aligner，用于将翻译的DNA查询序列映射到蛋白质参考数据库（BLASTX比对模式）。根据数据和设置，短读时BLAST的加速高达20,000，相对于BLAST的典型灵敏度为90-99％。但是，nr数据库的搜索时间仍然很长（主要是通过参考数据库的大小来确定少量reads的时间）。

如果您自己运行DIAMOND，则可以使用以下命令：

mkdir -p dmnd_tmp
diamond blastx -p 4 -d nr -q mouse1_contigs_unmapped.fasta -o mouse1_contigs.dmdout -f 6 -t dmnd_tmp -k 10 --id 85 --query-cover 65 --min-score 60
diamond blastx -p 4 -d nr -q mouse1_unassembled_unmapped.fasta -o mouse1_unassembled.diamondout -f 6 -t dmnd_tmp -k 10 --id 85 --query-cover 65 --min-score 60

命令行参数是：

-p：搜索中使用的线程数为4。
-q：输入文件名。
-d：数据库名称。
-e：保存匹配的期望值（E）阈值。
-k：要保留的最大比对序列数为10。
t：临时文件夹。-o：输出文件名。
-f：输出文件为表格格式。

从这些搜索的输出中，您需要使用以下脚本提取最匹配的蛋白质。在这里，如果85％的序列同一性超过reads长度的65％，则我们考虑匹配-这可能会导致非常差的e值（E = 3！），但是匹配仍然是合理的。

./7_Diamond_Protein_Map.py nr mouse1_contigs_map.tsv mouse1_genes_map.tsv mouse1_proteins.fasta mouse1_contigs_unmapped.fasta mouse1_contigs.dmdout mouse1_contigs_unannotated.fasta mouse1_unassembled_unmapped.fasta mouse1_unassembled.dmdout mouse1_unassembled_unannotated.fasta

该脚本的参数结构为：

7_Diamond_Protein_Map.py <Protein_database> <Contig_Map> <Gene_Map> <Output_File_For_Protein_sequences> <Unmapped_Contigs_File> <Contig_Diamond_Output> <Output_For_Unannotated_Contigs> <Unmapped_ Unassembled_Reads_File> <Unassembled_Reads_Diamond_Output> <Output_For_Unannotated_Unassembled_Reads>由于非冗余蛋白质数据库包含许多物种（包括真核生物）的条目，因此我们经常发现序列reads可以匹配具有相同分数的多种蛋白质。当前，从这些多个匹配中，我们选择第一个（即“热门”）。正如在宏基因组学讲座中提到的那样，可以应用更复杂的算法，但是我们目前的理念是，共享相同序列匹配的蛋白质在任何情况下都可能具有相似的功能；分类法是一个单独的问题！

「问题13：每个步骤映射了多少次reads？reads映射到多少个基因？基因定位到多少蛋白质？」

BWA = 3356次reads的映射reads总数定位基因总数（BWA）= 1234 DIAMOND = 51936个reads的映射reads总数映射的蛋白质总数（DIAMOND）= 21699 因此，在〜83000个假定的微生物mRNA来源的reads中，我们可以将〜55000个注释中的〜23000个基因注释！

请记住，要提取此步骤的预先计算的输出文件：

tar -xzf precomputed_files.tar.gz mouse1_genes_map.tsv mouse1_genes.fasta mouse1_proteins.fasta

步骤10.酶功能注释

为了从功能的角度帮助解释我们的转录组学数据集，我们依赖于将数据映射到功能性网络，例如代谢途径和蛋白质复合物图。在这里，我们将使用KEGG碳水化合物的代谢途径。

首先，我们需要首先将注释的基因与KEGG途径中的酶进行匹配。为此，我们将使用Diamond来从SWISS-PROT数据库中识别已分配酶功能的基因/蛋白质的同源物。菱形是通过同源性注释酶功能的相对粗略而直接的方法。我们选择在这里使用它是为了避免引入其他工具。但是，文献中存在更强大的酶功能注释方法，例如我们自己的基于概率密度的酶功能注释工具Detect。

mkdir -p dmnd_tmptar -xzf precomputed_files.tar.gz swiss_db.dmnd swiss_map.tsv通过我们的BWA搜索确定的微生物基因：

diamond blastx -p 4 -d swiss_db -q mouse1_genes.fasta -o mouse1_genes.diamondout -f 6 -t dmnd_tmp -e 10 -k 1对于通过我们的DIAMOND搜索确定的蛋白质：

diamond blastp -p 4 -d swiss_db -q mouse1_proteins.fasta -o mouse1_proteins.diamondout -f 6 -t dmnd_tmp -e 10 -k 1然后，我们需要生成一个映射文件，其中列出了我们的基因/蛋白质和酶(EC)号，描述了与之相对应的酶功能：

./8_Gene_EC_Map.py swiss_map.tsv mouse1_genes.diamondout mouse1_proteins.diamondout mouse1_EC_map.tsv该脚本的参数结构为：

8_Gene_EC_Map.py <SWISS-PROT_EC_Mappings> <Diamond_Output_For_Genes> <Diamond_Output_For_Proteins> <Output_EC_Mapping_File>

「问题14：我们的数据集中识别出多少种独特的酶功能？」

cat mouse1_EC_map.tsv|wc -l #1016

步骤11.生成与每个基因相关的标准化表达值

我们已尽我们所能删除了低质量的碱基/reads，载体，接头，接头，引物，宿主序列和rRNA序列，并注释了reads-现在让我们总结一下发现。我们通过观察我们每个基因在微生物组中的相对表达来做到这一点。

./9_RPKM.py nodes.dmp mouse1_classification.tsv mouse1_genes_map.tsv mouse1_EC_map.tsv mouse1_RPKM.txt mouse1_cytoscope.txt注意事项：

该脚本的参数结构为：

9_RPKM.py <nodes_file> <kaiju_Classification> <Gene_Mapping_File> <EC_Mapping_File> <RPKM_Output> <Cytoscope_Attributes_Output>输出文件的结构mouse1_RPKM.txt为：

[geneID/proteinID, length, #reads, EC#, Total RPKM, RPKM per phylum]gi|110832861|ref|NC_008260.1|:414014-415204 1191 1 3.9.3.5 106.98 0 0 45.89 6.86 20.77 7.35 2.3 0 4.63「问题15：看一下mouse1_RPKM.txt文件。最高表达的基因是什么？哪个门似乎最活跃？」

WP_016228325.1和Bacteroidetes

步骤10.使用KEGG Pathway作为Cytoscape的支架可视化结果。

为了在碳水化合物代谢途径的背景下可视化我们处理过的微生物组数据集，我们使用了网络可视化工具-Cytoscape以及EnhancedGraphics和KEGGscape插件。下面提供了一些用于加载网络，节点属性和更改视觉属性的有用命令（在线提供了许多cytoscape教程）。

下载代谢途径

首先，使用以下命令从KEGG下载碳水化合物代谢途径：

wget https://github.com/ParkinsonLab/2017-Microbiome-Workshop/releases/download/EC/ec00010.xmlwget https://github.com/ParkinsonLab/2017-Microbiome-Workshop/releases/download/EC/ec00500.xml您可以找到其他[KEGG上的路径]（http://www.genome.jp/kegg-bin/get_htext?htext=br08901.keg），也可以通过选择Download KGML页面顶部的选项将其导入Cytoscape 。每个途径。

安装Cytoscape插件

选择Apps->App Manager 搜索 enhancedGraphics 在中间列中选择，然后单击Install右下角搜索 KEGGScape KEGGScape在中间列中选择，然后单击Install右下角将XML从KEGG导入Cytoscape(我表示这个网络和速度实在太差了，用了好久才安装完全)

选择File-> Import-> Network->File... 选择XML文件，ec00010.xml或者ec00500.xml单击Open 检查Import pathway details from KEGG Database框，然后选择OK 加载节点属性文本文件（.txt）-这会将属性映射到网络中的节点，您随后可以将它们可视化

选择File-> Import-> Table->File... 选择mouse1_cytoscape.txt文件，然后单击Open Key Column for network从更改shared name为KEGG_NODE_LABEL 点击确定可视化您的节点属性

在左侧Control Panel选择Style标签勾选Lock node width and height方块单击Size面板最左侧的框，然后将默认节点大小更改为20.0。单击您单击的框右侧的空白框，以更改默认大小，将Column字段更改为RPKM，并将Mapping Type字段更改为Continuous Mapping 单击Image/Chart 1面板最左侧的框，切换至Charts选项卡，单击“甜甜圈”图标，然后>>在两列字段之间按“添加全部”按钮，然后单击“应用”（确保从添加的字段中删除整体RPKM到甜甜圈圈）如果看不到Image/Chart 1面板，请从控制面板的左上角选择Properties-> Paint-> Custom Paint 1->Image/Chart 1 为了改善可视化效果，您可以在Image/Chart 1-> Charts-> 下修改颜色属性Options，或修改其他属性，例如标签字体大小，标签位置，填充颜色，节点位置和边缘属性注意事项：

为了方便起见，提供了预先计算了节点属性的cytoscape文件tar -xzf precomputed_files.tar.gz Example.cys，可以随时打开它并以不同的可视化效果和不同的布局进行播放-例如，将圆形布局与spring嵌入式布局进行比较。如果要返回原始布局，则必须重新加载文件。Cytoscape可能很气质。如果看不到节点的饼图，它们将显示为空白圆圈，您可以手动显示它们。在左侧的“属性”面板下，有一个标记为“自定义图形1”的条目。双击左侧的空白框（这是默认行为）-这将弹出一个新窗口，其中包含“图像”，“图表”和“渐变”选项-选择“图表”，选择所需的图表类型（饼图或甜甜圈），然后将不同的细菌分类从“可用列”移到“选定列”。最后，点击窗口右下方的“应用”（在您移动窗口之前可能看不到）。