单细胞数据库数据挖掘,单细胞数据处理

　　最近我和一个朋友交流了纳米孔的16S测序数据分析，发现我们真的没有从零开始做过一次这个数据分析，然后发现这方面的资料比较少，就研究了一下，分享给大家。坦白说，牛津纳米孔测序技术在16S多样性的研究上还是有一些不足，只能说勉强够用。主要应用场景是在一些现场快速检测，主要是致病菌。但相信随着测序精度和分析软件的提高，会越来越多地使用。得益于互联网的便利性和分享精神，今天我们可以轻松获取测序的原始数据，并自由分析。

　　最后的内容主要是指从它的一个兄弟的博客，https://ithelp.ithome.com.tw/articles/10200625之前提到过。

　　1.软件和数据准备处理牛津纳米孔的测序数据，首先当然要安装相关的专用软件，基本都是原厂的，原厂品质，放心！有minimap2和yacrd，用于嵌合体。数据来自一篇文章，文件不大。可以直接下载，也可以通过ascp下载。

　　#此过程请参考https://github.com/tetedange13/stage_M2BI_scripts#软件准备#首先使用NanoPlot检查测序质量，直接安装pip。要提高速度，可以使用清华pip install NanoPlot #的porechopgit clone，以及https://github.com/rrwick/Porechop.gitcd porechoppython 3 setup . py的nano filt install porechop-h # quality control。或者pip bioconda安装pip install nanofilt#minimap2，编译或下载预编译的二进制文件。在这里，编译git clone https://github.com/lh3/minimap2cd minimap 2 make # yacrd，一个install conda install yacrd#来获取数据。我是ascp下载/volumes/10.11/users/ZD 200572/applications/aspera \ Connect.app/Contents/Resources/ascp-I ~/asper aweb _ id _ DSA . OpenSSH-Tr-Q-L 100m-p 33001-L-fasp-g1k @ fasp . 1000 genomes . ebi . AC . uk:vol 1/fastq/err 277/07/err 2778177/err 2778177 . fastq . gz . #比较下载时请检查gi号。term term=33175[生物工程]#] #在这里，我已经把这个文件保存在网盘上了。请在微信官方账号回复“16S”获取下载地址链接：https://pan.baidu.com/s/1FOp6kU2dB_37_Wi_mH2F0Q抽取码：8gce 2。数据预处理基本上是一个质量控制过程。首先，检查测序质量。当然，纳米孔的质量还是低了点，尤其是从手里下载的数据是测序核心低配版的R9.4。未来，通过串联两个纳米孔，R10可以提高到95%。然后去除测序的接头，得到真正的测序序列。该不该去掉底漆？但估计数据库中的序列也应该不含引物，所以估计的引物影响不大。然后，进行过滤以去除明显不合格的序列。

　　#查质量，fastq最常用。这里试试官方的nano plot-T 4-fastq err 2778177 . fastq . gz-Plots hex dot。结果仍然是一个类似于FastQ质量控制报告的功能，但是缺少了一些统计指标。有两个很好的图表将读取的长度和质量分布放在一起。

　　# Unjoint，4线程，这一步需要几个小时。对我五年前的双核四线程笔记本电脑PoreChop-T4-I Input进行质量控制。FastQ-O调整。FastQ #，质量值为9，长度为200过滤。虽然有点低，但是有两个重定向nanofilt-q 9-l200修剪。fastq已清理。fastq #摆脱嵌合体。首先，minimap2进行全局比较，找到嵌合体。然后yacrd去掉minimap 2/minip2-x ava-ont-g 500 cleaned . fastq cleaned . fastq overlap . pafyacrd-I overlap . PAF-o report . yacrd-c 4-n 0.4 scrubb-I cleaned . fastq-o reads . SCR ubb . FQ #卡在这里卡的时间最长。需要准备几个主要的文件，需要整理的已经放在文件夹里了。剩下的可以用ncbi下载，NCBI很大，但是不放，或者直接用我建的。#看看文件的保留情况，大约一半的数据树-h grep q ERR2778177.fastq.gz-[224m]-[237m]被清理了。FastQ -[200m]读。斯克罗布。-[476m] Trimmed.fastq 3。最后比较了纳米孔16S的数据分析，之前翻译的综述总结了大概有两种。一种是像前面的16S数据分析一样对ASV/OTU进行聚类，但是因为90%左右的准确率，所以用85%的准确率聚类是有用的。另一种是不聚类直接比较，也就是我们这次学到的blast等比较工具的方法。根据综述作者的观点，这些工具无法完成物种注释的任务(不够准确)，因为它们不是专门为纳米孔测序数据设计的。笔者推荐minimap2和离心机。但是，很大一部分文章都是基于这些工具进行数据分析的。大家动手吧最后做个对比！

　　#建立参考数据库last-1060/src/lastdb-uyass-R01微生物数据库NCBI-16s.fasta #比较https://giga science . biomed central . com/articles/10.1186/s 13742-016-0111-z # sec 12 last-1060/src/LASTAL-P4-Q1-B1-Q0-A1-e45微生物数据库已清理。FastQ已对齐。MAF #转换为blast格式last-1060/scripts/Maf-转换BLAST制表符对齐。mafaligned.txt # # # 4。结果处理，获得物种丰富度表。因为我的水平还不足以写脚本从比对结果中获取物种标注，所以脚本基本上来自台湾省的那个同事。确切的说，应该还是对结果的数据结构不太熟悉，争取以后再熟悉一些。为了实现他的脚本的使用，对数据库的信息进行了一些小的提取操作，这些操作都是用python来完成的，我暂时用的还挺方便的。我已经把处理过的数据上传到百度云了。直接用的话可以跳过。-- cat NCBI-16s . fasta grep conv . txt #先获取序列名信息，备用` ` ` ` ` Python DIC _ FA={} I=0 with open。=: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:在第行：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()-我是分割线。

　　# 加载包，第一次看到这样的方式隐藏消息(库(数据。表))抑制消息(库(整齐韵文))#设置文件位置id _ mapping _ file-/Volumes/Data/bio Data/nano pore-16S/info。txt b输出-/体积/数据/生物数据/纳米孔-16S/对齐。txt #导入输入文件id _ mapping-fread(id _ mapping _ file，header=FALSE，sep=\t ，string as factors=FALSE)setkey(id _ mapping，V1)b输出-弗雷德(b输出，header=FALSE，sep=\t ，string as factors=FALSE)b输出_已过滤-b输出% % DP lyr:filter(V3=90 V4=700)% % group _ by(V1)% % DP lyr:filter(V3==max(V3))% % DP lyr:filter(v1循环获得信息for(I in 1:length(refseqIds)){ refseqid-refseqIds[I]tmp-id _ mapping[.(refseqid)]Lineage-tmp $ V3 OrgName-tmp $ V4 readnum-readStats[refseqid]PERC-readnum * 100/allreadsnum df $ Lineage[I]-Lineage df $ ReadsNum[I]-readnum df $ reads PERC[I]-PERC df$OrgName[I]-OrgName } taxonomy-paste(df $ Lineage，df $ OrgName，sep=；)output _ df-data。frame(Taxonomy=Taxonomy，reads number=df $ reads num)output _ df-aggregate(output _ df $ reads number，by=list(Taxonomy=output _ df $ Taxonomy)，FUN=sum)colnames(output _ df)[2]--reads number read percentage-output _ df $ reads number * 100/all reads numput _ df _ df-cbind(output _ df，read percentage)output _ df[order(output _ df $ read percentage，降序=TRUE\\/Alignment\\/，，file))write.table(output_df， out.txt ，sep=\t ，row.names=FALSE，quote=FALSE)最后，如果顺利，就是结果了：

　　分类单元编号百分比266沙漠链霉菌菌株at4-8-cv 159新球聚球菌菌株sy 1459/11b 11布拉氏锥虫菌株中心菌株at6-11-rm4 1585864867

郑重声明：本文由网友发布，不代表盛行IT的观点，版权归原作者所有，仅为传播更多信息之目的，如有侵权请联系，我们将第一时间修改或删除，多谢。