Metagenomic analysis reveals oropharyngeal microbiota alterations in patients with COVID-19

哈工大省医

宏基因组分析揭示了新冠病人口腔微生物的改变。

31 covid-19，29 flu，28 normal 三组对照。

宏基因组数据先区分出metagenomics，和COVID-19，reads，然后mapping to covid-19 reference.

组装出whole genome 看看snp，哪些位点是否有变异？

virus discovery？

看鉴别出是哪种virus

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7441049/

鉴定出的COVID-19 样本：定的标准是什么？

讨论：

meta对疾病的恶化，和预后的影响。

[]metaquast 提供reference

[]prokka参数 --kingdom 中Bacteria 和Viruses的区别

########################################################

用contig跟coronavirus reference做blast，

找到所有的coronavirus的contig，并画图：

cd /USB/metagenomics/analysis/my_project/blastn_coronavirus_result

python  gather_blastn_contig.py /USB/metagenomics/analysis/my_project/blastn_coronavirus_result /USB/metagenomics/analysis/my_project/megahit_result >all_sample.coronavirus.contigs.fa

下一步用mega或什么对齐软件画个树看看。

保留这个结果，并且验证跟kraken2的结果的一致性：

#######################################################

braken build database

/USB/tools/Bracken-2.6.2/bracken-build -d /USB/meta_genomics_pipeline/database/kraken2_database -t 32 -l 250

#######################################################

troubleshooting for kraken2 results

Question：kraken2的结果中有48% unclassified，并且没有找到virus相关的任何比例。

Loading database information... done.

275881 sequences (55.16 Mbp) processed in 10.696s (1547.6 Kseq/m, 309.44 Mbp/m).

  142595 sequences classified (51.69%)

  133286 sequences unclassified (48.31%)

解决：在kraken2构建完成的database中seqid2taxid.map文件里面没有找到virus的序列，故此推断在构建时没有把virus包含进去，因为重新使用kraken2-build 构建数据库。Note：重新构建数据库前必须要删掉原来的seqid2taxid.map，taxo.k2d，hash.k2d，opts.k2d四个文件，

第二次：after rebuild the database，重新跑kraken2:

$ /USB/meta_genomics_pipeline/bin/kraken2 --db /USB/meta_genomics_pipeline/database/kraken2_database  --threads 50  --report /USB/metagenomics/analysis/my_project/kraken2_result/houmaohu.kraken2.report --use-names --output /USB/metagenomics/analysis/my_project/kraken2_result/houmaohu.kraken2.output /USB/metagenomics/analysis/my_project/kneaddata_result/houmaohu.kneaddata.fastq

Loading database information... done.

275881 sequences (55.16 Mbp) processed in 3.922s (4221.0 Kseq/m, 843.98 Mbp/m).

  174591 sequences classified (63.28%)

  101290 sequences unclassified (36.72%)

finally we found Covid-19 in kraken2 result. yeah. but still 36.72% of sequences are unclassified. we doubt that they are human related sequencing .

从第二次的结果中能够找到virus结果了，但是仍然有36%是unclassified，验证是否是human genome，发现不是。

验证方法如下：

使用bowtie2 将用于kraken2的fastq跟hg19 reference比对，发现比对率为0.说明该fastq已经不含有human genome了。那么怀疑是其它fungi，plasmid，等数据库，为了验证，

需要下载其他数据库，并重新使用kraken2-build 构建index。

第三次更新数据库之后

$ /USB/meta_genomics_pipeline/bin/kraken2 --db /USB/meta_genomics_pipeline/database/kraken2_database  --threads 32  --report /USB/metagenomics/analysis/my_project/kraken2_result/houmaohu.kraken2.report --use-names --output /USB/metagenomics/analysis/my_project/kraken2_result/houmaohu.kraken2.output /USB/metagenomics/analysis/my_project/kneaddata_result/houmaohu.kneaddata.fastq 2>/USB/metagenomics/analysis/my_project/kraken2_result/houmaohu.kraken2.log

Loading database information... done.

275881 sequences (55.16 Mbp) processed in 9.124s (1814.2 Kseq/m, 362.74 Mbp/m).

  174783 sequences classified (63.35%)

  101098 sequences unclassified (36.65%)

still 36.65% unclassified

那么这些序列是什么呢？输出来看看吧，重新用kraken2 输出unclassified reads。在ncbi上

发现还是virus和一些bacteria，那下载nt库试试看吧。

kraken2-build --download-library nt --db ./

######################################################################

kraken2 report explaination: here

######################################################################

'a' : all taxonomic levels

'k' : kingdoms

'p' : phyla only

'c' : classes only

'o' : orders only

'f' : families only

'g' : genera only

's' : species only

######################################################################

记录cornonavirus的相关信息：

单独的cov-reference放在：/USB/meta_genomics_pipeline/database/SARS-CoV-2-reference

NCBI名称：

>NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome

长度：29903bp

######################################################################

添加contig注释功能：目的：鉴定组装出来的contig都是什么物种。目标：能够找到target物种。用于后续的分析。

思路：基于kraken2使用的bacteria，virus，fungi，archaea，等6个数据库的fasta文件，重新构建blastn 所需index，然后使用blastn 将sample组装好的contigalign到该数据库上。使用in-house脚本统计并整理每个contig的物种信息，并汇总。

1）构建数据库：

cd /USB/meta_genomics_pipeline/database/blastn_dbs

cat ../kraken2_database/library/*/library.fna > my_kraken_cat6dbs_20210726.fna

makeblastdb  -parse_seqids  -dbtype  nucl  -in  my_kraken_cat6dbs_20210726.fna

2）mapping 添加到main pipeline中，做一个测试：

blastn -db /USB/meta_genomics_pipeline/database/blastn_dbs/my_kraken_cat6dbs_20210726.fna -query /USB/metagenomics/analysis/my_project/megahit_result/houmaohu/houmaohu.megahit.final.contigs.fa -out ./sample01.corona.out -outfmt 7 -num_threads 30 -subject_besthit

3）写脚本统计每条contig都是什么物种（besthit）

输入文件1：contig文件，

输入文件2：blast out 文件，cat 所有blast out文件到一个文件中，

输入文件3：blast 时用到的reference文件中就有id，叫做seqid2taxid.6abfpuv.map

输出文件1：contig fasta

思路二：

无需cat 6个database，直接分别比对，这样节省时间，提高运算速度，只是一个组装文件会得到6个比对结果，再用脚本整理到一起

1）那么对6个database分别构建blastn index

2）分别比对6次

3）统计结果

ls /USB/meta_genomics_pipeline/database/kraken2_database/library/*/library.fna |while read line;do name=`ls $line|awk -F '/' '{print$(NF-1)}'`;echo "makeblastdb  -parse_seqids  -dbtype  nucl  -in $line -out ./$name ";done |sh

有个数据库出错了，检查blastn的报错，并重新构建blast index，

######################################################################

遗留问题：

[]下一步用mega或什么对齐软件画个树看看。

[x]rebuild kraken2 database (running)

[x]rerun kraken2 to check ratio after 2 was done

[x]rerun blastn when finished build reference index.

[x]kraken2 结果画图，不画了，krona的图就够了。

[x]处理结果统计脚本，完成，将fasta文件用于galaxy展示。

[x]检查blast index 的错误，重跑blast

######################################################################

使用krona画物种分类图：

注意：kraken2.output中第三列为taxa_id,so do not using --use-names when using kraken2.

le ./houmaohu.kraken2.output |cut -f 2,3 >houmaohu.kraken.krona

ktImportTaxonomy houmaohu.kraken.krona -o houmaohu.taxonomy.krona.html

[ WARNING ] Too many query IDs to store in chart; storing supplemental files in 'taxonomy.krona.html.files'.

会报一个warning，但不影响html上的图

done

######################################################################

验证kneaddata之后的reads是否还有hg19.genome.

$ bowtie2  -p 20 -x /USB/meta_genomics_pipeline/database/kneaddata_database/hg37dec_v0.1 -U ../kneaddata_result/houmaohu.kneaddata.fastq -S houmaohu_r1_bowtie.fw.sam --nofw --local

275881 reads; of these:

  275881 (100.00%) were unpaired; of these:

    275870 (100.00%) aligned 0 times

    1 (0.00%) aligned exactly 1 time

    10 (0.00%) aligned >1 times

0.00% overall alignment rate

(prokka_env) ceid 15:08:09 /USB/metagenomics/analysis/my_project/test_kneaddata_2hg19

$ bowtie2  -p 20 -x /USB/meta_genomics_pipeline/database/kneaddata_database/hg37dec_v0.1 -U /USB/metagenomics/analysis/my_project/fastqc_result/houmaohu.raw.fastq.gz -S houmaohu_raw_bowtie.fw.sam --nofw --local

488767 reads; of these:

  488767 (100.00%) were unpaired; of these:

    474782 (97.14%) aligned 0 times

    3171 (0.65%) aligned exactly 1 time

    10814 (2.21%) aligned >1 times

2.86% overall alignment rate

1）原始fastq中就只有2.8%的human reads，

2）kneaddata后没有human reads，证明kneaddata会去除humanreads

3）那么问题来了，kraken2记过中37% no hit是什么？是什么？是什么？

############################################################################

看VIP工具文章，写思路，出结果。

############################################################################

data from paper:

https://www.nature.com/articles/srep23774/tables/2

补充数据分析，新增PE数据，适配流程。

seqkit fx2tab -lg UniVec_Core_library.fixed.fna |cut -f 1,4,5|head -1

python /USB/meta_genomics_pipeline/bin/blastout2contig_annotation.py -c /USB/metagenomics_analysis/analysis/my_project/megahit_result/houmaohu/houmaohu.megahit.final.contigs.fa -b /USB/metagenomics_analysis/analysis/my_project/blastn_result/houmaohu/houmaohu.cat.all.blast.out -m /USB/meta_genomics_pipeline/database/blastn_dbs/seqid2taxid.6abfpuv.map -o /USB/metagenomics_analysis/analysis/my_project/blastn_result/houmaohu.annotated.final.contigs.fa

contig_name  annotation  ref_name    ref_gc    contig_len    contig_gc Genome_coverage

k141_259        Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome    29903  37.97  7433    39.63  24.857%

个性化分析：

提取所有样本中的target 物种的组装结果序列，并检查N50，N90等，看看全长多少，能不能画进化树之类的，或者后续snp分析。

################################################

method for phylogeny：

conda create -n phylo blast mafft raxml iqtree bedtools

conda activate phylo

mafft --auto all.coronavirus.with.ref.fa > all.coronavirus.with.ref.out

raxmlHPC -s all.coronavirus.with.ref.out -m GTRGAMMA -n coronavirus_tree -p 12345

# or

iqtree -s all.coronavirus.with.ref.out

cat all.coronavirus.with.ref.out.treefile

paste treefile to https://itol.embl.de/upload.cgi

图不直观，应该用gene,

用contig 注释出gene，在用gene作图。

补充：

1）Genome Coverage (%)  ，contig长度与reference长度的比例，在已有的脚本上修改。

2）Total reads  （过滤kneaddata后的fastq的reads，）

3）Number of reads (%) 统计比对上每个contig后的reads，咋统计？

评估每种virus的结果。

参数修正：

megahit --min-len-contig 1000

blastn -evalue 10

reference:

https://genomics.sschmeier.com/ngs-orthology/index.html