納米孔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表為空，或許是原始數據有問題，或許前面參數出了問題，做了個探索。
