现在手里有一批NGS测序，经初步比对，95%的序列是细菌，比对到人的特别低，因为建库的时候没有做去rRNA处理，导致绝大部分序列是rRNA序列，几乎占据了90%。所以现在的分析思路是：
1# 对rRNA部分进行taxonomy分析，阿尔法、beta分析
2# 对非rRNA部分（暂时就叫mRNA）,做功能注释、DEGs、富集分析, kegg，ko
3# 对每个样本病原体分析
4# 对某病毒的覆盖度、深度统计
5# 对病毒覆盖度>80%, 深度>100X的，做intra-host variant分析
首先，做数据过滤，因为是外面测序公司给的，他们已经做了基本过滤了，这一步我就省了；
接着，用sormeRNA比对到自带的几个ref，拆分出rRNA和other的序列：
sortmerna -ref xxx1 -ref xxx2 -ref xxx3 ...-reads R1.fq.gz -reads R2.fq.gz --fastx --aligned sortmenra.aligned --other --paired_in --num_alignments 1 -m 1000000 -v
注意：-m 参数我也不是很清楚怎么设，反正这一步比较慢
得到other.fq和sortmeRNA.aligned.fq
因为我的是PE reads，输出的在一个文件里，于是写个脚本把reads拆分为两个文件：
perl split_fq.pl other.fq
perl split_fq.pl sortmeRNA.aligned.fq
因为我们还没有去人源，理论上人源的在other.fq部分，于是把这一部分比对到human（bwa-mem2 mem）, 然后筛选出非比对到人的reads：
grep -vE "@SQ|@RG|@PG" all_other_map2huaman.sam|awk '$3~/\*/'|awk '{print $1}' >all_other_nonhuman.list
写个脚本根据reads ID提取fastq序列：
perl extract_reads_based_on_ID.pl all_other_nonhuman.list all_other.R1.fq.gz all_other.R2.fq.gz
接着对去人源后的lnRNA部分和rRNA部分分别做kraken2注释：这就需要安装kraken2和braken2了。
Manual · DerrickWood/kraken2 Wiki (github.com)
首先下载kraken2的数据库，看文献，需要包含archaea、bacteria、protozoa、fungi、human、virus（包括SARS-Cov-2）。可以用kraken2直接下载标准库：
kraken2-build --standard --db $DBNAME （为你要下载到的路径）
也可以参考Index zone by BenLangmead 看这里的数据库，直接下载：
wget -c https://genome-idx.s3.amazonaws.com/kraken/k2_standard_20200919.tar.gz
注意：网速太慢了，耗时很长
kraken2-build --download-taxonomy --db $DBNAME （为你要下载到的路径）它其实是在$DBNAME下面建立taxonomy目录，里面放nucl_gb.accession2taxid.gz（这个暂时存疑，我还没下完）文件，这里展示的是下载中样子：
然后需要用braken（一个字母的差别，看manual，像同一家写的）建库：
braken-build -d $DBNAME -t 24 -k 35 -l 94 -x /WhereYourCommanderInstalled/kraken2
注意：-k 35是我看到文献里一般都这么设置的，-l 94是因为我的测序读长是PE94，默认是100，根据自己实际的测序读长设置。它也会要求$DBNAME下面有library目录，要不然会报错，library下面是taxonomy，里面放着seqid2taxid.map。会生成三个文件：
database94mers.kmer_distrib
database94mers.kraken
database.kraken
于是seqid2taxid.map这个文件哪里来的？在kranken2的官网搜一下这个文件，然后直接下载：
wget -c http://ccb.jhu.edu/software/bracken/dl/seqid2taxid.map
因为服务器上有kraken2的数据库（商家自带的，他们自己整理的5w+条目，从GTDB下载的；不是标准库），暂时先用这个测试下：
/biostack/tools/microbiome/kraken2-2.0.9b/kraken2 --db /biostack/database/kraken2/kraken-db --threads 12 --report all_sortmeRNA-kraken2.report -output all_sortmeRNA-kraken2.output --paired all_sortmeRNA.R1.fq.gz all_sortmeRNA.R2.fq.gz --gzip-compressed --report-zero-counts
注意哦：这里加了--report-zero-counts参数
对kraken2的结果用braken跑个report（注意：这里有-l参数，为读长，实际测序为PE94，但没有对应长度的databaseXXXmers.kmer_distrib和databaseXXXmers.kraken文件），这里还是用-r 100试一下：
/biostack/tools/microbiome/Bracken-2.5/bracken -d /biostack/database/kraken2/kraken-db -i mRNA-kraken2.report -o mRNA.bracken -r 100 -l S -t 12
默认是100，商家也有150、200的
再对结果用krona画个图瞅瞅：
提取kraken2结果的ID和taxonomyID两列：cat mRNA-kraken2.output|cut -f 2,3 >mRNA-kraken2.output.krona
再作为krona的输入文件：
/project/baowenjuan/APP/Krona-2.7.1/KronaTools/bin/ktImportTaxonomy mRNA-kraken2.output.krona -n Bacteria
提示以上ID没有找到
我也不知道为什么。。。。在debug。。。。
2020年12月08日：
今天测试用metaphlan：
安装：conda install -c bioconda python=3.6 metaphlan
测试跑：metaphlan --bowtie2db /usr/local/lib/python3.6/site-packages/MetaPhlAn-3.0.4-py3.6.egg/metaphlan/metaphlan_databases --bowtie2out rRNA.bowtie2.out --input_type fastq all_sortmeRNA.R1.fq.gz,all_sortmeRNA.R2.fq.gz -o rRNA.metaphlan.o
提示错误：
Use of uninitialized value $bt2_args[2] in join or string
可以忽略，这里说的：[Warnings] MetaPhlAn version 3.0 · Issue #101 · biobakery/MetaPhlAn (github.com)
但是比对率特别低:3.09% overall alignment rate。可能因为数据是采用的marker gene吧。anyway，把文献下载了现看看。
晚上看完文献，发现这个不适用于我的数据（参见我写的另外一篇 文章正在审核中... - 简书），嘻嘻嘻。白干~明天继续折腾！
网友评论