Nanopore 16S测序数据分析流程之blast/last
最近有朋友和我交流纳米孔16S测序数据的分析,发现真的没有从头完成过一次这方面的数据分析,然后发现这方面的资料也比较少,于是学习一下,和大家分享。坦白说,牛津纳米孔测序技术在16S多样性研究方面还是有些不足的,只能说勉强够用,主要应用场景是在一些现场快速检测方面,主要是病原菌这种。但是,相信随着测序准确度的提高和分析软件的改进,相信它的应用会越来越多。感谢互联网的便利和分享精神,今天的我们可以方便地获得测序的原始数据,并可以自由进行分析。
last方面的内容主要参考自IT帮一个兄弟的博客,之前提到过https://ithelp.ithome.com.tw/articles/10200625
1.软件和数据准备
处理牛津纳米孔的测序数据,首先当然要安装相关的专用软件了,基本上原厂出的,原厂品质,放心!还有就是minimap2和yacrd,用于去嵌合。数据来自一篇文章,文件不大,直接下载或者ascp下载均可。
#这个流程可参考https://github.com/tetedange13/stage_M2BI_scripts
#软件准备
#首先是NanoPlot,用于检测测序质量,直接pip安装了,加速的话可使用清华源
pip install NanoPlot
#去接头的Porechop
git clone https://github.com/rrwick/Porechop.git
cd Porechop
python3 setup.py install
porechop -h
#质控的NanoFilt,pip或bioconda安装
pip install nanofilt
#minimap2,编译或者下载预编译的二进制文件均可,这里编译
git clone https://github.com/lh3/minimap2
cd minimap2 && make
#yacrd,conda安装
conda install yacrd
#获取数据,我是ascp下载
/Volumes/10.11/Users/zd200572/Applications/Aspera\ Connect.app/Contents/Resources/ascp -i ~/asperaweb_id_dsa.openssh -Tr -Q -l 100M -P33001 -L- fasp-g1k@fasp.1000genomes.ebi.ac.uk:vol1/fastq/ERR277/007/ERR2778177/ERR2778177.fastq.gz .
#比对工具last,一款来自日本的比对软件,需要编译
wget http://last.cbrc.jp/last-1060.zip
unzip last-1060.zip
cd last-1060
make
#blast diamond seqkit
conda install -c bioconda diamond blast seqkit
#如果wget或者浏览器下载地址是wget http://ftp.ebi.ac.uk/vol1/fastq/ERR277/007/ERR2778177/ERR2778177.fastq.gz
#ncbi的16S数据库,由于出国网络差,下了好久,一定要确保数据下载完整,注意下载时勾选gi号,以备后面有用
#https://www.ncbi.nlm.nih.gov/nuccore?term=33175%5BBioProject%5D
#这里我把这个文件存在了网盘,公众号内回复“16S”即可获得下载地址
链接: https://pan.baidu.com/s/1FOp6kU2dB_37_Wi_mH2F0Q 提取码: 8gce
2.数据预处理
基本上是质控的过程,先看下测序质量,当然纳米孔的质量还是有点低的,特别是手上下载的数据是低版本的测序芯处R9.4,未来的R10可以通过两个纳米孔串联提高到95%。接着就是去除测序的接头,获得真正的测序序列,是不是引物也应该切除?但是估计数据库里的序列也应该不包含引物,所以估计引物影响不大。然后,进行过滤,除去明显不符合要求的序列。
#查看质量,fastq最常用,这里用官方的试试
NanoPlot -t 4 --fastq ERR2778177.fastq.gz --plots hex dot
结果依然是类似于fastq质控报告的一个函数,不过统计指标少了几个,有两个把reads长度和质量分布放在一起的图不错。
#去接头,4线程,这一步耗时以小时计,对于我的五年前的双核心四线程笔记本
porechop -t 4 -i input.fastq -o trimmed.fastq
# 质控,质量值为9,长度200过滤,虽然有点低,这有两个重定向
NanoFilt -q 9 -l 200 < trimmed.fastq > cleaned.fastq
#去嵌合,先minimap2自全局比对,发现嵌合,然后yacrd去除
minimap2/minimap2 -x ava-ont -g 500 cleaned.fastq cleaned.fastq > overlap.paf
yacrd -i overlap.paf -o report.yacrd -c 4 -n 0.4 scrubb -i cleaned.fastq -o reads.scrubb.fq
#卡在这时间最长,主要几个文件要准备,需要整理的已经放在文件夹,其余的可以ncbi下载,较大,没放,或者直接用我建好的。
#看下文件的保留情况,大概去除了一半以上的数据
tree -h | grep q
├── [237M] ERR2778177.fastq.gz
├── [224M] cleaned.fastq
├── [200M] reads.scrubb.fq
└── [476M] trimmed.fastq
3. last比对
关于纳米孔16S的数据分析,之前翻译的那篇综述总结了大概有两种,一种是和之前的16S数据分析一样聚类ASV/OTU,但是由于90%左右的准确度,看有用85%的准确度聚类的。。。另一种就是,不聚类,直接进行比对,也就是我们这次学习的blast等比对工具的方法,根据综述作者的观点,这几种工具由于不是专门为纳米孔测序数据设计,不能比较好的完成物种注释的任务(不够准确),作者推荐的是minimap2和centrifuge这两个软件。但是,有很大一部分文章是以这几个工具进行数据分析的,我们也做一下,最后进行一个比较吧!
# 构建参考数据库
last-1060/src/lastdb -uYASS -R01 microbialdb ncbi-16S.fasta
# 比对https://gigascience.biomedcentral.com/articles/10.1186/s13742-016-0111-z#Sec12
last-1060/src/lastal -P 4 -q 1 -b 1 -Q 0 -a 1 -e 45 microbialdb cleaned.fastq > aligned.maf
# 转换成BLAST格式
last-1060/scripts/maf-convert blasttab aligned.maf > aligned.txt### 4.结果处理,获得物种丰度表由于自己水平还不够写脚本从比对结果中获得物种注释,脚本基本上来自台湾的那个同仁的。确切的说应该还是对比对结果的数据结构不够熟悉,争取后面多熟悉些。为了实现使用他的脚本,对数据库的信息进行了一些小的提取操作,用我暂时用得比较顺手的python完成的。我已经把处理好的数据上传了百度云,直接使用的话可以略过。--------------我是分界线的开始-----------`cat ncbi-16S.fasta | grep ">" > conv.txt #先获得序列名信息备用````python
dic_fa = {}
i = 0
with open('conv.txt') as f:for line in f:i += 1if line.strip() != '':ac_n = line.strip().split('|')[3]gi_n = line.strip().split('|')[1]dic_fa[ac_n] = ''dic_fa[ac_n] = [gi_n, line.strip()]
print(i)j = 0
fout = open('info.txt', 'w')
with open('sequence.gb') as f1:for line in f1:if line[:7] == "VERSION":seq_name = line.strip().split(' ')[1].split(' ')[0]if '/db_xref="taxon:' in line:tax_id = line.strip().split('/db_xref="taxon:')[1].split('"')[0]if seq_name in dic_fa.keys():fout.write(dic_fa[seq_name][1] + '\t' + tax_id + '\n') #+ seq_name + '\n')seq_name = ''tax_id = ''j += 1
print(j)fout.close()
--------------我是分界线的结束-----------
# 加载包,第一次看到这样的方式
suppressMessages(library(data.table))
suppressMessages(library(tidyverse))
# 设置文件位置
id_mapping_file <- "/Volumes/Data/Biodata/Nanopore-16S/info.txt"
blastOutputs <- "/Volumes/Data/Biodata/Nanopore-16S/aligned.txt"
#导入输入文件
id_mapping <- fread(id_mapping_file, header = FALSE, sep = "\t", stringsAsFactors = FALSE)
setkey(id_mapping, V1)blastoutput <- fread(blastOutputs, header = FALSE, sep = "\t", stringsAsFactors = FALSE)
blastoutput_filtered <- blastoutput %>% dplyr::filter(V3>=90 & V4>=700) %>% group_by(V1) %>% dplyr::filter(V3 == max(V3)) %>% dplyr::filter(V12 == max(V12))
allreadsnum <- length(unique(blastoutput$V1))
readStats <- table(blastoutput_filtered$V2)
refseqIds <- names(readStats)
df <- data.frame(RefseqID = refseqIds, Lineage = NA, ReadsNum = NA, OrgName = NA, ReadsPerc = NA)
#循环获得信息
for (i in 1:length(refseqIds)) {refseqid <- refseqIds[i]tmp <- id_mapping[.(refseqid)]lineage <- tmp$V3orgname <- tmp$V4readnum <- readStats[refseqid]perc <- readnum*100/allreadsnumdf$Lineage[i] <- lineagedf$ReadsNum[i] <- readnumdf$ReadsPerc[i] <- percdf$OrgName[i] <- orgname
}taxonomy <- paste(df$Lineage, df$OrgName, sep="; ")
output_df <- data.frame(Taxonomy = taxonomy, ReadsNumber = df$ReadsNum)
output_df <- aggregate(output_df$ReadsNumber, by = list(Taxonomy=output_df$Taxonomy), FUN = sum)
colnames(output_df)[2] <- "ReadsNumber"
ReadPercentage <- output_df$ReadsNumber*100/allreadsnum
output_df <- cbind(output_df, ReadPercentage)
output_df <- output_df[order(output_df$ReadPercentage, decreasing = TRUE),]
#output <- paste0("tax_", sub("\\.\\/Alignment\\/", "", file))
write.table(output_df, 'out.txt', sep = "\t", row.names = FALSE, quote = FALSE)
最后,如果顺利,就是结果了:
Taxonomy | ReadsNumber | ReadPercentage | |
---|---|---|---|
266 | Trichocoleus desertorum strain ATA4-8-CV2 | 1046 | 41.85674270 |
159 | Neosynechococcus sphagnicola strain sy1 | 248 | 9.92396959 |
267 | Tychonema bourrellyi strain CCAP 1459/11B | 110 | 4.40176070 |
111 | Kastovskya adunca strain ATA6-11-RM4 | 104 | 4.16166467 |
168 | Okeania plumata strain FK12-27 | 74 | 2.96118447 |
120 | Loriellopsis cavernicola strain LF-B5 | 53 | 2.12084834 |
40 | Cephalothrix komarekiana CCIBt 3277 | 42 | 1.68067227 |
33 | Caedimonas varicaedens strain 221 | 41 | 1.64065626 |
12 | Aliterella antarctica strain CENA408 | 33 | 1.32052821 |
118 | Limnoraphis robusta strain CCALA 966 | 31 | 1.24049620 |
Nanopore 16S测序数据分析流程之blast/last相关推荐
- fastq质量值_微生物组16S rRNA数据分析小结:从fastq测序数据到OTU table
推荐阅读 1.ggplot2绘制曼哈顿图示例2.phyloseq | 用 R 分析微生物组数据及可视化3.R语言PCA分析教程 | Principal Component Methods in R4. ...
- Medicine in Microecology:Nanopore三代测序人类肠道病毒组的方法
应用Nanopore三代测序技术解析人类肠道病毒组 Profiling of Human Gut Virome with Oxford Nanopore Technology Medicine in ...
- Medicine in Microecology:微生物所王军组发表Nanopore三代测序人类肠道病毒组的方法
文章目录 应用Nanopore三代测序技术解析人类肠道病毒组 摘要 引言 结果 1. 病毒组分离富集和测序 2. ONT测序揭示健康个体的病毒组成 3. 使用ONT序列进行病毒基因组组装 4. 使用O ...
- 学习全基因组测序数据分析2:FASTA和FASTQ
本文转载自微信公众号解螺旋的矿工,作者为黄树嘉,已获得授权.黄树嘉写了WGS系列的文章,堪称教科书级别的生物信息学习材料.虽然本平台只关注宏基因组领域,但此系列文章知识体系完善.干货满满,是值得每位专 ...
- 【重磅综述】长序列数据分析相关资源哪里找?一文读懂长序列测序数据分析的机遇与挑战!...
简介 标题:长序列测序数据分析的机遇与挑战 杂志:GenomeBiology 影响因子:10.806 发表时间:2020年05月08日 ...
- activiti自己定义流程之Spring整合activiti-modeler5.16实例(四):部署流程定义
注:(1)环境搭建:activiti自己定义流程之Spring整合activiti-modeler5.16实例(一):环境搭建 (2)创建流程模型:activiti自己定义流程之Spr ...
- 三代测序数据分析之文献推荐
三代测序数据分析之文献推荐 (2018-07-20 09:03:47) 转载▼ 分类: 文献推荐 1:下面这一篇是使用了目前可用的技术测序一个植物基因组并比较优劣 Paajanen P, Kett ...
- activiti自定义流程之Spring整合activiti-modeler5.16实例(四):部署流程定义
注:(1)环境搭建:activiti自定义流程之Spring整合activiti-modeler5.16实例(一):环境搭建 (2)创建流程模型:activiti自定义流程之Sprin ...
- [原创]FOCUS处理系统流程之:流程批量生成(个人专用懒人版)
根据标准流程生成整个工区流程文件 (1)可按按目录生成 (2)可替换其中的关键字 标准流程准确,则生成的流程将一点到底,所剩的操作变成鼠标的操作,懒之产物. 界面如下: ☆其它物探处理原创软件相关☆ ...
- [原创]FOCUS处理系统流程之:大文本文件极速合并(sps文件合并)
下载地址:csdn 软件界面1: (功能如题) : 1.普通文本文件的合并 2.大文件快速合并,空行及重复行删除等功能 3.合并物探 ...
最新文章
- linux基础 云,云计算之linux基础一
- Boost1.62.0 + VS2015 配置
- 【JavaSE04】Java中循环语句for,while,do···while-思维导图
- Neko and Aki's Prank
- 牛客多校6 - Binary Vector(组合数学+推公式)
- 美一8岁华裔男童体育课上头部重伤 家长吁调查
- ios 数字键盘左下角添加按钮_iOS8数字键盘加左下角完成button
- java的sql的like_[Java教程]SQL like 模糊查询, in
- jwt重放攻击_4个点搞懂JWT、JWS、JWE
- IIS6.0应用程序池回收和工作进程【转:http://www.cnblogs.com/freshman0216/archive/2008/06/02/1212460.html】...
- 只有单杀技能的飞鸽传书
- Principles of Reactive Programming 之Actors are Distributed (3)
- 金蝶KIS专业版13.0视频教程
- if语句 power query_判断(if)语句
- javaCRC8计算的坑
- 360网站域名拦截检测 非法网址检测系统原理
- day69_淘淘商城项目_02_dubbo介绍 + dubbo框架整合 + zookeeper + 商品列表查询实现 + 分页 + 逆向工程_匠心笔记
- 模型推荐丨政务大数据项目案例模型分享
- 分组统计group by
- 【软件工程1916|W(福州大学)_助教博客】团队第四次作业(第7次)成绩公示...
热门文章
- 06_JavaEE回顾笔记Ⅱ
- x509证书cer格式转pem格式
- Matlab的自相关函数corr
- 欧拉角与四元数互相转换
- bc547可以用8050代换吗_代换S8550 S8050三极管要特别注意放大倍数
- 小米盒子 smb Android,客厅里的多媒体 小米盒子SMB本地连接
- 红米k30 允许调用gpu调试层_高效渲染!RTX 3090卡皇打造NVIDIA STUDIO强力主机实战体验|nvidia|显卡|gpu|cpu|内存...
- [06]项目实战-移动端流体布局
- 十分钟搞清字符集和字符编码
- [雨林木风][番茄花园][电脑公司][深度论坛][龙帝国]系统光盘收录大全(精品)