Skip to content

前期准备

本次分析流程准备这样做，使用fastp进行质控，再使用bwa进行比对到人类基因组后去除人源reads。使用centrifuge或kraken2对进行微生物检出，最后使用pavian进行结果可视化。

数据准备

本次使用的数据来源于文章PMID32020836，样本类型是支气管肺泡灌洗液。提取样本的总RNA，然后使用试剂盒建库后使用Miseq进行了双端150bp测序。事实上，目前市面的mNGS（病原微生物方向）以75bp读长的单端检测为主，一般检测数据量为20m的reads，如果需要更高深度，一般会调整数据量到100m。使用mNGS的方式检测出了新冠病毒。

1
2
3
4
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR109/001/SRR10903401/SRR10903401_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR109/001/SRR10903401/SRR10903401_2.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR109/002/SRR10903402/SRR10903402_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR109/002/SRR10903402/SRR10903402_2.fastq.gz

上述文章的分析流程是在质控后，去除人源，再比对到病毒的参考数据库以及ncbi的nr库，使用blast来比对到病毒，然后使用Metaphlan2来识别病原细菌。再使用Megahit对上诉blast到病毒结果的reads进行组装。

还有以下文章PMID32615925，样本类型有尿液和血液，检测分析得到Enterococcus faecalis。所有数据可在SRA中获取。

参考数据准备

对于使用的centrifuge与kraken2两个软件，均需要参考数据库。

人类参考基因组

使用ucsc的hg38就好了。

1
2
wget https://hgdownload.cse.ucsc.edu/goldenpath/hg38/bigZips/hg38.fa.gz
gunzip hg38.fa.gz

centrifuge

由于centrifuge提供的数据库最新的版本是2016年的，提供的nr库也是2018年的，比较老旧，因此我们需要自行下载新的数据库进行更新。

建立archaea、bacteria、viral参考库 

1
2
3
4
5
6
7
centrifuge-download -o taxonomy taxonomy
centrifuge-download -o library -m -d "archaea,bacteria,viral" refseq > seqid2taxid.map
cat library/*/*.fna > input-sequences.fna
centrifuge-build -p 8 --conversion-table seqid2taxid.map \
    --taxonomy-tree taxonomy/nodes.dmp \
    --name-table taxonomy/names.dmp \
    input-sequences.fna abv

但是，因为网络原因，使用以上流程下载会经常性的失败（其实我感觉是因为NCBI把数据库从原来的ftp协议转移到了https协议，但是centrifuge里没有修改，我也不知道在哪改）。幸好，我们也可以使用genexa提供的数据库，最新版本是2020年3月29日的，足够我们使用。

1
2
wget https://zenodo.org/record/3732127/files/h+p+v+c.tar.gz?download=1
tar -zxvf h+p+v+c.tar.gz

还是建议挂一个代理来下载，快很多。

建立nt库（并不一定需要）

1
2
3
4
5
6
7
8
9
wget https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nt.gz
gunzip nt.gz && mv -v nt nt.fa
wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/gi_taxid_nucl.dmp.gz
gunzip -c gi_taxid_nucl.dmp.gz | sed 's/^/gi|/' > gi_taxid_nucl.map
centrifuge-build -p 8 --bmax 1342177280 \
    --conversion-table gi_taxid_nucl.map \
    --taxonomy-tree taxonomy/nodes.dmp \
    --name-table taxonomy/names.dmp \
    nt.fa nt

像nt库就没办法了，因为genexa已经停止了这部分的分享，毕竟一直占着别人的带宽也不好，还是得去ncbi下载。ncbi的nt库目前已经有100G大了，对于在自家进行练习并不合适，因此先不下载了，后续也会跳过此步。另外一个原因是gi_taxid_nucl.dmp.gz这个文件我还没有找到新的获得链接。

kraken2

下载kraken2使用的数据库，由于我的电脑内存只有16G，因此使用minikraken的数据库。

1
wget ftp://ftp.ccb.jhu.edu/pub/data/kraken2_dbs/minikraken_8GB_202003.tgz

如果需要自行建立数据库，可下载archaea、bacteria、plasmid、viral、human、fungi、plant、protozoa、nr、nt、env_nr、env_nt、UniVec等。 

1
2
3
4
kraken2-build --download-taxonomy --threads 24 --db databases/
kraken2-build --download-library viral --threads 24 \
    --db databases/
kraken2-build --build --threads 24 --db databases/

比对索引准备

流程上使用bwa或bowtie2进行比对。一般来说，常用的是bowtie2，因此这里对人类基因组建立bowtie2的索引。

1
bowtie2-build hg38.fa hg38

数据质控

使用fastp进行数据质控。

1
2
3
4
5
6
fastp -i rawdata/SRR10903401_1.fastq.gz \
    -I rawdata/SRR10903401_2.fastq.gz \
    -o cleandata/SRR10903401_1.clean.fastq.gz \
    -O cleandata/SRR10903401_2.clean.fastq.gz \
    -w 8 -l 36 -p \
    -j QC/SRR10903401.json -h QC/SRR10903401.html

参考

Metagenomics

杭州市疾控中心测序实验室