纳米孔Nanopore-16S数据分析学习笔记

1.下载原始数据

本次学习分析的文章是这篇：https://academic.oup.com/gigascience/article/7/12/giy140/5202451
这篇文章的原始数据有点问题，使用sra和ena数据库直接下载都基本上会失败，sra只能下到一个10M左右的数据，转换格式成fastq后只能获得4.6M的数据。最后使用aspera connect下载可以成功。命令如下，我是黑果，其他系统格式应该类似，软件安装和使用参见我前面的学习记录。https://jiawen.zd200572.com/916.html

ascp -QT -l 3000m -P33001 -L- -k1 -i asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/run/ERR224/ERR2241542/RAW_R9.5_1_organism.tar.gz BioData/Nanopore/

下载的数据没有通过完整性测试，没办法进行下一步分析。后面重新下载了原始数据，见下面：

2.解压转换成fastq

#解压
tar zxvf RAW_R9.5_1_organism.tar.gz
#hd5转换成fastq

1.把下载的fastq格式转为fasta

其实这是一个很简单的过程，即使自己用个脚本或者使用命令行也能解决，介于想要重复作者结果，就按作者的原步骤进行。

#新建一个python2的conda环境
conda create -n python2 python=2
#安装依赖
conda install  -y  biopython=1.65 blast seqtk numpy vsearch mafft-ginsi
# seqtk转换fastq到fasta
seqtk seq  ERR2241540.fastq -a > ERR2241540.fasta

2.获得共识序列

这里走了点弯路，其实本文的参考文献里说明了是使用INC-seq这个流程进行前处理的，找到这个流程的github仓库，就可以使用了。

git clone https://github.com/CSB5/INC-Seq.git
############我是分界线，以下是我走的弯路#######
#安装三代纠错软件，卡在这好久，docker下载会卡住失败，源码编译会下载一个不存在的boost_1_58_0-headersonly.tbz2文件，最后用conda解决（友情提示，只有linux-64版本）
conda install pbdagcon -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
#########这些分界线包着的是不需要的步骤，都包含在了INC-seq的脚本里####
#开始获得共识序列#
./inc-seq.py -i  ERR2241540.fasta -o consensus.fa  --copy_num_thre 3  --iterative  -a poa
#INC-seq的流程包含了三个比对方法（blastn,  graphmap, poa），默认是采用blastn的，依照文章作者的参数，--copy_num_thre 3  --iterative ，只有poa方法能够获得共识序列，blastn提示：
    #Warning: [blastn] Examining 5 or more matches is recommended
    #Warning: [blastn] Examining 5 or more matches is recommended
    #Warning: [blastn] Examining 5 or more matches is recommended
    #Warning: [blastn] Examining 5 or more matches is recommended
    #PBDAGCON failed (trimmed more than 100 bases)!
    #Consensus construction failed!
#graphmap提示Candidate read found!
    #PBDAGCON failed (trimmed more than 100 bases)!
    #Consensus construction failed!
#完成后的文件列表在这里，当然，只有poa有结果。可以看出序列利用率是比较低的，好多序列由于长度不够，或者发现的片段不一致而过滤掉了。
tree -h 
.
├── [ 183]  bpipe.config
├── [   0]  consensus-blastn.fa
├── [   0]  consensus-gra.fa
├── [150K]  consensus-poa.fasta
├── [4.0K]  data
│   └── [ 17K]  inc_seq_test_read.fa
├── [2.3M]  ERR2241540.fasta
├── [4.6M]  ERR2241540.fastq
├── [7.4K]  inc-seq.py
├── [1.4K]  LICENSE
├── [1.2K]  pipeline.bpipe
├── [3.0K]  README.md
└── [4.0K]  utils

3.共识序列方向校正

awk -v k="Sample1" '/^>/{gsub(">","",$0);$0=">barcodelabel="k";"$0}1' chop.fasta > Sample1_barcode.fa

4.获得otu表

python nanoCLUST.py -i Sample1_barcode.fa -o ./nano -s 0,450,451,900,901,1300
比较遗憾，尽管没有报错，却没有得出结果，otu表为空，或许是原始数据有问题，或许前面参数出了问题，做了个探索。