Skip to content

前期准备

本次分析流程准备这样做,使用fastp进行质控,再使用bwa进行比对到人类基因组后去除人源reads。使用centrifugekraken2对进行微生物检出,最后使用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

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