Nanopore 纳米孔 测序数据处理 微生物 16S全长 Centrifuge的安装和使用

linux Centrifuge 纳米孔测序 微生物 16S全长

EPI2ME在线网站处理

差不多2021.6这之前不用注册，后面要注册了。而且只对购买仪器的开发-555-

操作流程官方学习网站github

操作流程官方学习网站Centrifuge指南

官网也描述到用Centrifuge分类物种

软件权威文章推荐

本地服务器处理如下流程：

① 拿到公司提供给你fastq格式文件，进行barcode和质量过滤，裁剪，这部分参考宏基因组处理吧。我自己用了NanoFilt处理，FastQC质量评估

gunzip -c 00data/Control105A.fastq.gz | NanoFilt -q 7 -l 500 --maxlength 1780 --headcrop 75 --tailcrop  50 | gzip > 01nanofilt/clean_Control105A.fastq.gz

② 下载安装 github下载

网络不行就直接下载，解压后，拉进去服务器，体积小，很快。
切换到拉进去的文件夹 设置安装在的目录 /home/XXX/

cd centrifuge
make
make install prefix=/home/xxx

设置PATH环境变量，也是启动软件的命令

export PATH=$PATH:/home/xxx/bin

③ 根据需求下载参考的序列的数据库 indexes下载地址

下面画红色的地方，找到你对应要的，16S直接右边的bacteria
== 建议用迅雷破解版的等等下载，15M每秒，比较快==

解压

tar -xf  /home/xxx/nanopore/p_compressed_2018_4_15.tar.gz

④ 运行 代码解释看上面提到官方学习网站

centrifuge -p 8 --time --ignore-quals -S control105_result.tsv --report-file control105_report.tsv -x /home/xxx/nanopore/NCBI/p_compressed -U /home/xxx/dataset/nanopore/Nanopore/Control105A.fastq.gz

结果文件

⑤ 自己找的参考数据库建立 indexes ，我用口腔微生物数据库{HOMD数据库}，但是它无法提供要求的那三个文件，可以用centrifuge1.0.4\example\reference的文件代替，所以报告只有result，而report.tsv的结果是空的，需要自己提取。

centrifuge-build --conversion-table /home/zhuangzc/dataset/nanopore/taxonomy/merged.dmp --taxonomy-tree /home/zhuangzc/dataset/nanopore/taxonomy/nodes.dmp --name-table /home/zhuangzc/dataset/nanopore/taxonomy/names.dmp /home/zhuangzc/dataset/nanopore/HOMD_16S_rRNA_RefSeq_V15.22.aligned.fasta homd

cd 05result_homd
centrifuge -p 12 --time --ignore-quals --min-hitlen 22 -S Control105A_result.tsv --report-file Control105A_report.tsv -x /home/zhuangzc/dataset/nanopore/db3/homd -U /home/zhuangzc/dataset/nanopore/01nanofilt/clean_Control105A.fastq.gz

Github讨论里的一些理解 特别是如果处理result.tsv的

关于score

abundance的计算

⑥ 自己HOMD数据库用Centrifuge的report提取出单个样本的tax和num

lf <-list.files(pattern = "result.tsv$") #以report.tsv 结尾的
lf
files <- gsub("\\.tsv", "", lf)   #切掉后缀.tsv，获得这些名称，为循环准备
files
for (i in seq_along(files))
  assign(files[i], read.table(lf[i], sep = '\t', header = TRUE
                              ,fill = T,quote=""))
Control105A_result <-filter(Control105A_result,hitLength>=500) 
Control105A_result <- Control105A_result[,c('seqID','numMatches')]

HOMD <- read.delim("~/HOMD.taxonomy.txt", header=FALSE)
Control105A_result <- left_join(Control105A_result,HOMD,by=c('seqID'='V1')) 
Control105A_result <- tapply(Control105A_result$numMatches,Control105A_result$V2,sum)
Control105A_result <- as.data.frame(Control105A_result)
names(Control105A_result)[1] <- "Control105A"
Control105A_result <- rownames_to_column(Control105A_result,var='tax')
write_csv(Control105A_result,'zhujie.csv')

不能直接去掉那些hit小于500的，这些应该定义为unassigned。

##【保留<500计数的】1
#library(tidyverse)
#Control105A_result <- mutate(Control105A_result,
#                             seqID=ifelse(hitLength<500,"lowhit",seqID))
#Control105A_result <- Control105A_result[,c('seqID','numMatches')]
#Control105A_result <- left_join(Control105A_result,HOMD,by=c('seqID'='V1')) 
#Control105A_result[is.na(Control105A_result)] <-"unssigned"
#Control105A_result <- tapply(Control105A_result$numMatches,Control105A_result$V2,sum)
#Control105A_result <- as.data.frame(Control105A_result)
#names(Control105A_result)[1] <- "Control105A"
#Control105A_result <- rownames_to_column(Control105A_result,var='tax')

bug提示，需要变量需要character

#【保留<500计数的】2
library(tidyverse)
Control105A_result$seqID <- as.character(Control105A_result$seqID)
Control105A_result <- mutate(Control105A_result,
                             seqID=ifelse(hitLength<500,"lowhit",seqID))
Control105A_result <- Control105A_result[,c('seqID','numMatches')]
Control105A_result <- left_join(Control105A_result,HOMD,by=c('seqID'='V1')) 
Control105A_result$V2 <- as.character(Control105A_result$V2)
Control105A_result <- mutate(Control105A_result,V2=
                               ifelse(is.na(Control105A_result$V2)==TRUE,
                                                          "unssigned",V2))
Control105A_result <- tapply(Control105A_result$numMatches,Control105A_result$V2,sum)
Control105A_result <- as.data.frame(Control105A_result)
names(Control105A_result)[1] <- "Control105A"
Control105A_result <- rownames_to_column(Control105A_result,var='tax')

⑦ result.tsv的理解：readID重复的次数等于numMatches的个数（所以我们提取result的报告需要将numMatches全部转化成1）；seqID对应准确的taxonomy（但readID的hit是不同的，通过过滤hit，将相同的taxonomy合起来），而相同的readID对应不同的seqID（注意他们的hit是相同的，这里说明序列可能错误比对上多个物种的可能，而numMatces说明这个序列seqID比对上的多个，{下图3}，说明分配到种水平是不准确的，只能处理到属水平）

所以需要将numMatches全部转化成1。然后在seqID匹配上taxonomy。

继续⑦, 去掉重复匹配到多个种水平的序列后，==测序量的read序列数。24690

所需要将重复匹配到种水平的限制在genus水平。

⑧ result.tsv最终处理代码 最终理论一个readID对应多个seqID，多个seqID对应一个物种，所以直接readID比对上taxonomy

data <- read.csv('data(留底编辑) - 副本.csv')
HOMD <- read.delim("~/HOMD.taxonomy.txt", header=FALSE)
data <- read.table('Control105A_result.tsv', sep = '\t', 
                   header = TRUE,fill = T,quote="")
library(tidyverse)

#重复不一致的处理
#【一个readID对应多个seqID，多个seqID对应一个物种,就是一个物种有多个seqID】
#【所以直接readID比对上taxonomy】
data <- left_join(data,HOMD,by=c('seqID'='V1')) #taxonomy的变量名为V2
table(is.na(data$V2)) #23个缺失值
#
data2 <- table(data$readID,data$V2)
data2[data2>0] <- 1
data3 <- as.matrix(data2)
data4 <- data3[which(rowSums(data3) >1 ), ] 
a=row.names(data4)
data5 <- data[which(data$readID%in% a),]
write.csv(data5,'查验.csv')
#导出csv分列tax检查【不同的readIDID对应相同的taxonomy】
Ndata5<- data5[!duplicated(data5$readID), ]
dim(Ndata5)
#
data5$readID=droplevels(data5$readID)
d <- as.data.frame(table(data5$readID))
d  #重复readID重与不一致readIDID种水平数据的数据,后续直接合并进去，
#其numMatches还是为1，物种采用genus

#重复的readID与对应【一致的taxonomy】

data6 <- table(data$readID,data$V2)
data6[data6>0] <- 1
data7 <- as.matrix(data6)
data8 <- data7[which(rowSums(data7) ==1 ), ] 
a=row.names(data8)
data9 <- data[which(data$readID%in% a),]
data9_2 <- data9[!duplicated(data9$readID), ]
write.csv(data9,'查验2.csv')#其numMatches还是为1，物种采用species
dim(data9_2)
dim(Ndata5)+dim(data9_2)+23
22839+1827+23 #缺失值#原始reads=24690
#合并处理

#hitLength 和V2的taxonomy保留就可以
Ndata5<- data5[!duplicated(data5$readID), ]
Ndata5_2 <- select(Ndata5,hitLength,V2)

#直接分割处理
Ndata5_3 <- Ndata5_2   %>% 
  separate(V2,  sep = ";s__",    
           into = c("V2", "deleteVar"))  
write.csv(Ndata5_3,'属水平分割的检查.csv') #导出检查下
Ndata5_4 <- select(Ndata5_3,-deleteVar)

data9_2 <- data9[!duplicated(data9$readID), ]
data9_3 <- select(data9_2,hitLength,V2)

dim(data9_3)+dim(Ndata5_2)

dataall <- rbind(Ndata5_4,data9_3)
dim(dataall)

#hit过滤处理
dataall$V2 <- as.character(dataall$V2)
dataall2 <- mutate(dataall,V2=ifelse(hitLength<300,"lowhit",V2))
table(is.na(dataall2$V2))
dataall2$num <- 1
#分组统计合并taxonomy
dataall3 <- tapply(dataall2$num,dataall2$V2,sum)
dataall4 <- as.data.frame(dataall3)
Control105A_result <- rownames_to_column(dataall4,var='tax')
names(Control105A_result)[2] <- "Control105A"
write_csv(Control105A_result,'Control105A_result.csv')

⑨ 我拿自己一个样本比对上面代码的结果，接近是一致的，但是我此时设置的hitLength为0，就是没有过滤任何小于多少bp的，提示是否在nanofilt的时候就应该过滤bp到1000以上，但是queryLength就算是1500的也有hitLength是小于300bp的存在，还分配到species水平

⒑一个问题bug？：很好奇，我比对上都1000bp以上了，在种水平还是分不出来？

---

【一些参考1】对于分类分配，将FASTQ文件上载到EPI2ME平台的FASTQ 16S协议上，按质量对reads进行过滤，然后使用BLAST将分类分配给NCBI数据库，最小水平覆盖率为30％，最小精度为77％作为默认参数（ONT）。
–query_cov，覆盖度：满足相似度的情况下，同时要求查询序列的覆盖度达到多少；
time qiime feature-classifier classify-consensus-blast
【一些参考2】参考[github的MaestSi/MetONTIIME]，可以专门在conda上面安装(https://github.com/MaestSi/MetONTIIME)，我自己参考里面的代码，采用qiime2软件做blast处理序列，做出了的结果只能比对上genus 上面的参考1是设置的参考

#导入数据
qiime tools import \
--type 'SampleData[SequencesWithQuality]' \
--input-path /home/dengqr/dataset/nanopore_qiime/manifest.tsv \
--input-format 'SingleEndFastqManifestPhred33V2' \
--output-path nanoporeSequences.qza

#数据换成质控后的
qiime tools import \
--type 'SampleData[SequencesWithQuality]' \
--input-path /home/dengqr/dataset/nanopore_qiime/manifest2.tsv \
--input-format 'SingleEndFastqManifestPhred33V2' \
--output-path nanoporeSequencesclean.qza

/home/zhuangzc/dataset/nanopore/01nanofilt/

#【解释】https://docs.qiime2.org/2021.8/plugins/available/vsearch/dereplicate-sequences/?highlight=qiime%20vsearch%20dereplicate%20sequences
#Dereplicate sequence data and create a feature table and feature representative sequences.

time qiime vsearch dereplicate-sequences \
--i-sequences nanoporeSequencesclean.qza \
--o-dereplicated-table table_tmp.qza \
--o-dereplicated-sequences rep-seqs_tmp.qza

#无参/从头聚类 聚类是按序列相似度99%的水平执行的，以创建100%的OTU
#https://docs.qiime2.org/2021.8/plugins/available/vsearch/cluster-features-de-novo/?highlight=qiime%20vsearch%20cluster%20features%20de%20novo
time qiime vsearch cluster-features-de-novo \
--i-sequences rep-seqs_tmp.qza \
--i-table table_tmp.qza \
--p-perc-identity 1 \
--o-clustered-table table.qza \
--o-clustered-sequences rep-seqs.qza \
--p-threads 12

#半有参聚类半无https://docs.qiime2.org/2021.8/tutorials/otu-clustering/?highlight=vsearch
qiime vsearch cluster-features-open-reference \
--i-table table.qza \
--i-sequences rep-seqs.qza \
--i-reference-sequences /home/dengqr/dataset/shao/10db_HOMD/homd_otus.qza \
--p-perc-identity 0.77 \
--o-clustered-table table-or-77.qza \
--o-clustered-sequences rep-seqs-or-77.qza \
--o-new-reference-sequences new-ref-seqs-or-77.qza

time qiime feature-classifier classify-consensus-blast \
--i-query rep-seqs-or-77.qza \
--i-reference-reads /home/dengqr/dataset/shao/10db_HOMD/homd_otus.qza \
--i-reference-taxonomy /home/dengqr/dataset/shao/10db_HOMD/99-taxonomy.qza \
--o-classification taxonomy_77.qza \
--p-perc-identity 0.77 \
--p-query-cov 0.3

qiime taxa barplot \
--i-table table-or-77.qza \
--i-taxonomy taxonomy_77.qza \
--m-metadata-file /home/dengqr/dataset/nanopore_qiime/metadata.tsv \
--o-visualization taxa-bar-plots_77.qzv

#序列描述、查看
qiime demux summarize \
--i-data rep-sequences.qza \
--o-visualization demux_summary.qzv

#table描述
qiime feature-table summarize \
--i-table t/table.qza \
--o-visualization t/table.qzv \
--m-sample-metadata-file /home/dengqr/dataset/nanopore_qiime/metadata.tsv

qiime feature-table tabulate-seqs \
--i-data t/rep-seqs.qza \
--o-visualization t/rep-seqs.qzv

#数据库
/home/dengqr/dataset/shao/10db_HOMD/99-taxonomy.qza
/home/dengqr/dataset/shao/10db_HOMD/homd_otus.qza

#https://docs.qiime2.org/2021.8/plugins/available/feature-classifier/classify-consensus-blast/?highlight=qiime%20feature%20classifier%20classify%20consensus%20blast

if [ "$CLASSIFIER_UC" == "BLAST" ]; then

#对于分类分配，将FASTQ文件上载到EPI2ME平台的FASTQ 16S协议上，按质量对reads进行过滤，然后使用BLAST将分类分配给NCBI数据库，最小水平覆盖率为30％，最小精度为77％作为默认参数（ONT）。
--query_cov，覆盖度：满足相似度的情况下，同时要求查询序列的覆盖度达到多少；

time qiime feature-classifier classify-consensus-blast \
--i-query t/rep-seqs.qza  \
--i-reference-reads /home/dengqr/dataset/shao/10db_HOMD/homd_otus.qza \
--i-reference-taxonomy /home/dengqr/dataset/shao/10db_HOMD/99-taxonomy.qza \
--o-classification taxonomy2.qza
--p-perc-identity 0.77 \
--p-query-cov 0.3 

qiime taxa barplot \
--i-table table.qza \
--i-taxonomy taxonomy2.qza \
--m-metadata-file /home/dengqr/dataset/nanopore_qiime/metadata.tsv \
--o-visualization taxa-bar-plots2.qzv