纳米孔Nanopore-16S数据分析学习笔记
- 2020 年 3 月 3 日
- 筆記
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 [email protected]:vol1/run/ERR224/ERR2241542/RAW_R9.5_1_organism.tar.gz BioData/Nanopore/
下载的数据没有通过完整性测试，没办法进行下一步分析。后面重新下载了原始数据，见下面：
wget -c ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/ERR/ERR224/ERR2241540/ERR2241540.sra
2.把下载的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
3.获得共识序列
这里走了点弯路，其实本文的参考文献里说明了是使用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
4.共识序列方向校正
python chopSEQ.py -i ../../consensus-poa.fasta -l 1300 -m 1450 -f AGRGTTTGATCMTGGCTCAG -r GGGCGGWGTGTACAAG > chop.fasta
5.获得otu表
awk -v k="Sample1" '/^>/{gsub(">","",$0); $0=">barcodelabel="k";"$0}1' chop.fasta > Sample1_barcode.fa
python nanoCLUST.py -i Sample1_barcode.fa -o ./nano -s 0,450,451,900,901,1300
比较遗憾，尽管没有报错，却没有得出结果，otu表为空，或许是原始数据有问题，或许前面参数出了问题，做了个探索。
