一擅这、獲取SRA號
登錄Genome Announcements網(wǎng)站(https://mra.asm.org/)持偏,搜索關(guān)鍵詞“bacteria genome SRA”墅垮,在搜索到的細(xì)菌基因組文章中選擇一篇可免,找到文章記載的SRA號。
以下面文章中的SRR號 SRR6466501 為例:
二龙优、下載SRA文件
SRA (Sequence ReadArchive)數(shù)據(jù)庫羊异,是用于存儲二代測序的原始數(shù)據(jù),包括454, Illumina, SOLiD, lonTorrent, Helicos 和CompleteGenomics。除了原始序列數(shù)據(jù)外球化,SRA現(xiàn)在也存在raw reads在參考基因的比對信息秽晚。
根據(jù)SRA數(shù)據(jù)產(chǎn)生的特點,將SRA數(shù)據(jù)分為四類:Studies-- 研究課題筒愚、Experiments-- 實驗設(shè)計赴蝇、Runs-- 測序結(jié)果集、Samples-- 樣品信息巢掺。
SRA文件句伶,SRA數(shù)據(jù)庫用不同的前級加以區(qū)分:
●ERP或SRP表示Sudics;
●SRS表示Samplo;
●SRX表示Hxpcritmcrnt;
●SRR表示Runs;
從SRA數(shù)據(jù)庫上用prefetch下載sra文件,輸入命令如下:
prefetch SRR6466501 #下載
結(jié)果如下:
軟件自動建立~/ncbi/public/sra文件夾陆淀,查看:
三考余、 Fastq-dump解壓
cd ~/ncbi/public/sra #進(jìn)入到sra文件夾下
fastq-dump --gzip --split-files SRR6466501.sra #解壓
結(jié)果如下:
四、Fastqc質(zhì)控
*準(zhǔn)備:Fastqc安裝轧苫,可參考我的簡書http://www.reibang.com/p/e0659f09288c
使用以下命令進(jìn)行質(zhì)控:
fastqc SRR6466501_1.fastq.gz
fastqc SRR6466501_2.fastq.gz
結(jié)果:
FastQC報告
打開:SRR6466501_1_fastqc.html
SRR6466501_2_fastqc.html
五楚堤、Trimmomatic去接頭
切除接頭序列和低質(zhì)量堿基,此處使用數(shù)據(jù)過濾軟件Trimmomatic含懊。
Trimmomatic 是一個廣受歡迎的Ilumina平臺數(shù)據(jù)過濾工具身冬。
處理數(shù)據(jù)速度快,主要用來去除Illumina 平臺的Fastq序列中的接頭岔乔,并根據(jù)堿基質(zhì)量值對Fastq進(jìn)行修剪酥筝。
支持多線程,有兩種過濾模式雏门,分別對應(yīng)SE和PE測序數(shù)據(jù)嘿歌。
用法:
單末端測序模式
在SE模式下,只有一個輸入文件和一個過濾后的輸出文件
java -jar [Trimmomatic所在的絕對路徑] SE-phred33 [輸入文件] [輸出文件] [動作1] [動作2]......
雙末端測序模式
在PE模式下茁影,有兩個輸入文件,正向測序序列和反向測序序列宙帝,過濾之后,輸出文件有四個兩端序列都保留的為paried -端序列過濾后被遺棄另- -端保留為unparied
java jar [Trimmomatic的絕對路徑] PE -phred33
第一個輸入文件(forward端)正向端
第二個輸入文件(reverse端)反向端
輸出文件[forward_ paried forward_unpaired reverse_ paired reverse_ _unpaired ] [動作1] [動作2]......
另外也支持phred-33和phred-64格式互相轉(zhuǎn)化募闲,不過現(xiàn)在絕大部分Illumina平臺的產(chǎn)出數(shù)據(jù)都轉(zhuǎn)為使用phred-33格式了茄唐。
*準(zhǔn)備:Trimmomatic安裝:
unzip Trimmomatic-0.38.zip
cd Trimmomatic-0.38
java -jar trimmomatic-0.38.jar
Trimmomatic過濾命令如下:
java -jar /home/gxx/Trimmomatic-0.38/trimmomatic-0.38.jar PE -phred33 SRR6466501_1.fastq.gz SRR6466501_2.fastq.gz ./output_forward_paired.fq.gz ./output_forward_unpaired.fq.gz ./output_reverse_paired.fq.gz ./output_reverse_unpaired.fq.gz ILLUMINACLIP:/home/gxx/Trimmomatic-0.38/adapters/TruSeq2-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 MINLEN:75
結(jié)果如下:
動作說明
ILLUMINACLIP:
過濾reads中的Illumina測序街頭和引|物序列并決定是否去除反向互補(bǔ)的R1/R2中的R2
參數(shù)說明( PE測序要注意最后-一個參數(shù), SE最后兩個參數(shù)不用設(shè)置,參數(shù)之間用冒號連接):
<接頭和引物序列所在位置( Trimmomatic自帶在adapters里面,其中TruSeq2為2代測序, TruSeq3為三代測序)> :
<seed搜索允許多少個個堿基錯配)>:<alindrome比對分值閾值為多少>:<simple clip比對分值閾值為多少>:
<palindrome模式允許切除的最短接頭序列為多少默認(rèn)為8bp>:
<palindrome模式去除與R1完全反向的R2>
SLIDINGWINDOW:
從reads的5'端開始,進(jìn)行滑窗質(zhì)量過濾,切掉堿基質(zhì)量平均值低于閾值的滑窗
參數(shù)說明: <窗口大小>: <堿基平均質(zhì)量值閾值>
MAXINFO:
一個自動調(diào)整的過濾選項,在保證reads長度的情況下盡量降低測序錯誤率,最大化reads的使用價值
LEADING:
從reads的開頭切除質(zhì)量低于閾值的堿基
TRAILING:
從reads的末尾切除質(zhì)量低于閾值的堿基
CROP:
從reads的末尾切掉部分堿基是reads達(dá)到指定長度
HEADCROP:
從reads的開頭切掉指定數(shù)量的堿基
MINLEN:
如果經(jīng)過剪切后reads的長度低于閾值則丟棄這條reads
AVGQUAL:
如果reads的平均堿基質(zhì)量值低于閾值則丟棄這條reads
TOPHRED33:
將reads的堿基質(zhì)量值體系轉(zhuǎn)為phred-33
TOPHRED64:
將reads的堿基質(zhì)量值體系轉(zhuǎn)化為phred-64
六蝇更、SPAdes組裝基因組草圖
SPAdes:
?由俄羅斯科學(xué)院圣彼得堡理工大學(xué)計算生物學(xué)實驗室開發(fā),是目前評價最好的拼接工具之一呼盆。
?主要用于基因組拼接年扩,也可用于一、二访圃、三代測序的混合組裝;還可用于轉(zhuǎn)錄組從頭組裝(rnaSPAdes)和宏基因組拼接(metaSPAdes) 厨幻。
*準(zhǔn)備:SPAdes的安裝
python環(huán)境下,安裝命令如下:
wget http://cab.spbu.ru/files/release3.12.0/SPAdes-3.12.0-Linux.tar.gz
mkdir ~/Biosofts/spades
tar zvxf /home/gxx/SPAdes-3.12.0-Linux.tar.gz -C ~/Biosofts/spades/
~/Biosofts/spades/SPAdes-3.12.0-Linux/bin/spades.py -h
echo 'export PATH=~/Biosofts/spades/SPAdes-3.12.0-Linux/bin:$PATH' >> ~/.bashrc
source ~/.bashrc
spades.py -h
SPAdes運行:
spades.py --careful --pe1-1 SRR6466501_1.fastq.gz --pe1-2 SRR6466501_2.fastq.gz -o ./SPAdesout_SRR6466501new
遇到錯誤:out of memory,內(nèi)存不足况脆。
解決:將虛擬機(jī)設(shè)置中的系統(tǒng)內(nèi)存調(diào)至合理范圍內(nèi)最大饭宾,再次嘗試以上命令,但最終還是不行格了。看铆。
嘗試使用seqtk抽取1000條,命令如下:
#解壓
gunzip -c output_forward_paired.fq.gz >output_forward_paired.fq
gunzip -c output_reverse_paired.fq.gz >output_reverse_paired.fq
#抽取1000條
seqtk sample -s 60 output_forward_paired.fq 1000 >seqtksample1_1000.fq
seqtk sample -s 60 output_reverse_paired.fq 1000 >seqtksample2_1000.fq
#用wc查看盛末,可對比前后文件弹惦,判斷是否抽取成功
wc -l output_forward_paired.fq
wc -l seqtksample1_1000.fq
然后,再次嘗試SPAdes運行:
spades.py --careful --pe1-1 seqtksample1_1000.fq --pe1-2 seqtksample2_1000.fq -o ./SPAdesout_SRR6466501_1000new
結(jié)果如下:
七悄但、Quast評價組裝的基因組效果
*準(zhǔn)備:Quast的安裝棠隐,安裝命令如下
tar zvxf quast-5.0.0.tar.gz -C ~/Biosofts/
~/Biosofts/quast-5.0.0
python quast.py -h
運行代碼:
quast.py ~/ncbi/public/sra/SPAdesout_SRR6466501_1000new/contigs.fasta -o ~/ncbi/public/sra/SPAdesout_SRR6466501_1000new/quast_out
結(jié)果如下:
查看輸出的文件夾quast_out:
本地下載quast報告 report.html,并查看: