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

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

软件权威文章推荐

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

① 拿到公司提供给你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 homdcd 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.qzatime 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.3qiime 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.tsvqiime 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%20blastif [ "$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

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

NBT：使用纳米孔测序从微生物组中得到完整闭环的细菌基因组
文章目录使用纳米孔测序从微生物组中得到完整成环的细菌基因组热心肠导读摘要前言结果图1 定义的12种细菌混合物中的序列分类学组成.每种细菌的读长分布和基因组组装图2:在两个健康的人类粪便微 ...
iMeta | 南科大夏雨组纳米孔测序揭示微生物可减轻高海拔冻土温室气体排放
点击蓝字关注我们基于纳米孔测序的宏基因组学揭示微生物作为生物过滤器减轻高海拔冻土的温室气体排放 https://doi.org/10.1002/imt2.24 Research Article V ...
使用纳米孔测序数据进行16S-DNA条形码研究的计算方法[综述]
摘要通过对16S核糖体RNA(16S rRNA)基因进行测序来评估细菌多样性已广泛用于环境微生物学中,特别是自从高通量测序技术问世以来.这些技术带来的另一项创新是需要开发新的策略来管理和研究生成的大 ...
MPB：微生物所王军组-人类肠道病毒粒子富集及纳米孔测序
为进一步提高<微生物组实验手册>稿件质量,本项目新增大众评审环节.文章在通过同行评审后,采用公众号推送方式分享全文,任何人均可在线提交修改意见.公众号格式显示略有问题,建议电脑端点击文末阅 ...
NBT：王运浩、区健辉等综述纳米孔测序技术
2021年11月8日,美国俄亥俄州立大学(Ohio State University)区健辉(Kin Fai Au)研究组在Nature Biotechnology在线发表综述论文Nanopore s ...
【不容错过】12月10日：纳米孔测序科研团队大会NCM 2020亚太区特别专场
一年一度的纳米孔测序科研团体大会(NCM 2020)主会场已于美国东部时间12月初在线上成功召开,汇集了全球超过50位领先的纳米孔测序学者,分享他们的纳米孔测序最新研究成果. 2020年12月10日, ...
中国首次纳米孔测序大会：不可错过的教学专场和技术诊断
纳米孔测序是由英国牛津Oxford Nanopore研发的最新一代高通量单分子测序技术,支持无扩增直接分析DNA/RNA,生成超长读长,实时测序能够即时实施数据分析.随着纳米孔技术在芯片.试剂.软件 ...
纳米孔测序揭示冻土冻融对土壤微生物群落变化的影响
抗生素耐药性是当今世界日趋严重的健康威胁,仅在美国每年就可导致2万人死亡.而鉴定特定的抗生素抗性微生物是迅速和恰当治疗的关键.就算是在缺乏暴露的情况下,环境中的细菌也能够扮演抗性基因储备库的角色.随着 ...
面对万亿级测序市场，纳米孔测序技术何去何从？
这是<肠道产业>第 482 篇文章 [直播预告]"Protein& Cell人类微生物组专刊线上论坛" 12月21日晚7点开播,敬请期待!(点击查看详情) 编者 ...
Cell二连发 | 广东CDC/耶鲁大学利用纳米孔测序揭示中/美新冠病毒基因组流行病学传播规律...
利用纳米孔测序技术实时测定病毒全基因组信息(Nanopore Real-time Sequencing),能够动态地分析病毒分子进化来研究病毒的变异及传播特征,这些信息对疫情发展不同阶段制定有效的防控 ...

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

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

EPI2ME在线网站处理

操作流程官方学习网站github

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

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

软件权威文章推荐

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

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

② 下载安装 github下载

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

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

结果文件

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

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

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

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

bug提示，需要变量需要character

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

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

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

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

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

最新文章

热门文章

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

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

EPI2ME在线网站处理

操作流程官方学习网站github

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

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

软件权威文章推荐

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

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

② 下载安装 github下载

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

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

结果文件

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

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

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

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

bug提示，需要变量需要character

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

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

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

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

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

最新文章

热门文章

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

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

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

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

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

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