一、獲取SRA號
登陸Genome Announcements網(wǎng)站（地址：https://mra.asm.org/），搜索感興趣的基因組，我搜索的是細菌基因組。
進入網(wǎng)站頁面后，點擊右上角的搜索符號。
我隨意選了一篇文章，是一篇菲律賓貝納姆河岸分離的六種細菌基因組序列草稿的文章。我在文章里找到了細菌信息。
我選取第四個細菌基因組為例。
SRR8491236
二、下載SRA序列
首先需要在linux下載SRAtoolkit，然后才能下載SRA數(shù)據(jù)。
可以參考我的簡書文章：
http://www.itdecent.cn/p/3df52aa857ad?v=1667697177169
然后再執(zhí)行命令將SRA數(shù)據(jù)下載下來：
prefetch SRR8491236
運行結(jié)果如下：
一般軟件會自動建立~/ncbi/public/sra文件夾，但我的不知什么緣由，經(jīng)常不會建立這個文件夾，所以我就手動建立了。
三、解壓SRA文件為fastq格式
命令如下：
fastq-dump --gzip --split-files SRR8491236
結(jié)果如下：
四、用fastqc進行數(shù)據(jù)質(zhì)量評價
這里需要安裝fastqc，可以參考我的簡書文章：
http://www.itdecent.cn/p/5d364dfc0f43
輸入以下命令：
fastqc SRR8491236_1.fastq.gz
fastqc SRR8491236_2.fastq.gz
結(jié)果如下：
每個測序樣本生成了兩個文件。
可以將html文件下載到本地來查看和分析。
FastQC Report解讀
FastqC有3種結(jié)果：綠色代表pass；黃色代表warn；紅色代表fail。
可以參考簡書文章：http://www.itdecent.cn/p/c81c7110eed4
五、Trimmomatic去接頭
Trimmomatic是一個廣受歡迎的Illumina 平臺數(shù)據(jù)過濾工具。
Trimmomatic支持多線程，處理數(shù)據(jù)速度快，主要用來去除Illumina 平臺的Fastq序列中的接頭，并根據(jù)堿基質(zhì)量值對Fastq進行修剪。
軟件有兩種過濾模式，分別對應(yīng)SE 和PE 測序數(shù)據(jù)，同時支持gzip和bzip2 壓縮文件。
Trimmomatic也支持phred-33 和phred-64 格式互相轉(zhuǎn)化，不過現(xiàn)在絕大部分Illumina 平臺的產(chǎn)出數(shù)據(jù)也都轉(zhuǎn)為使用phred-33 格式了。
Trimmomatic是基于Java開發(fā)的，因此需要提前安裝Java，才能使用Trimmomatic。
Trimmomatic安裝可以參考我的簡書文章：
http://www.itdecent.cn/p/8e854198db13
Trimmomatic過濾命令如下：
java -jar Trimmomatic-0.38/trimmomatic-0.38.jar PE -phred33 ~/ncbi/sra/SRR8491236_1.fastq.gz ./ncbi/sra/SRR8491236_2.fastq.gz ./output_forward_paired.fq.gz ./output_forward_unpaired.fq.gz ./output_reverse_paired.fq.gz ./output_reverse_unpaired.fq.gz ILLUMINACLIP:Trimmomatic-0.38/adapters/TruSeq2-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 MINLEN:75
命令執(zhí)行解釋：
PE -phred33 因為雙端測序需要輸入兩個文件。最終會輸出4個文件。其中兩個文件對應(yīng)于"paired"數(shù)據(jù)，即read1和read2都保留。另外兩個對應(yīng)于"unpaired"數(shù)據(jù)，在處理的過程中會過濾掉其中一端的reads。
~/ncbi/sra/SRR8491236_1.fastq.gz ./ncbi/sra/SRR8491236_2.fastq.gz 這里是輸入將要過濾的文件
./output_forward_paired.fq.gz ./output_forward_unpaired.fq.gz ./output_reverse_paired.fq.gz ./output_reverse_unpaired.fq.gz這里是輸出的四個結(jié)果文件，都保存在當前目錄下
ILLUMINACLIP:Trimmomatic-0.38/adapters/TruSeq2-PE.fa:2:30:10 切除TruSeq3-PE中提供的Illumina適配器。最初Trimmomatic將尋找種子匹配(16個堿基)，最多允許2個不匹配。如果在配對端讀取的情況下達到30分(約50個堿基，50*0.6)，或在單端讀取的情況下達到10分(約17個堿基, 17*0.6)，這些"seed"序列將被修剪。
SLIDINGWINDOW:5:20 掃描5個堿基寬滑動窗口，當每個堿基的平均質(zhì)量下降到20以下時進行剪切
LEADING:20 刪除前低質(zhì)量堿基或N堿基(低于質(zhì)量3)
TRAILING:20 刪除后低質(zhì)量堿基或N堿基(低于質(zhì)量3)
MINLEN:75 上述步驟完成后，刪除小于36個堿基的reads
參考文章：
https://www.cnblogs.com/Sunny-King/archive/2022/06/16/Bioinformatics-Trimmomatic.html
運行結(jié)果：
由結(jié)果可以知道，我的兩條雙端測序都保留的占88.40%，只有正向測序的序列保留的占8.39%，只有反向測序的序列保留的占1.09%，過濾掉2.13%。
六、Spades組裝基因組草圖
安裝Spades可以參考我的簡書文章：
http://www.itdecent.cn/p/f5b2aa3fda64?v=1667780854519
運行命令如下：
spades.py --careful --pe1-1 ~/ncbi/sra/SRR8491236_1.fastq.gz --pe1-2 ~/ncbi/sra/SRR8491236_2.fastq.gz -o ./SPAdesout_SRR8491236
命令解釋：
--careful 減少錯配和短插入失.讓程序運行MismatchCorrector----不建議用于中型和大型基因組。官方文檔中可以看到MismatchCorrector占用的disk space非常多。
注:默認情況下執(zhí)行read糾錯、組裝。
--pe1-1
--pe<#>-1 <file_name>:left reads, <#>代表第幾個文庫。
--pe<#>-2 <file_name>:right reads, <#>代表第幾個文庫。
note:只有一個文庫，指定#等于1即可。
-o 指定輸出文件目錄(必需)。
這里參考簡書文章：
http://www.itdecent.cn/p/f02475d61a5c
出現(xiàn)錯誤，內(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 -s60 output_forward_paired.fq 1000 > seqtksample1_1000.fq
seqtk sample -s60 output_reverse_paired.fq 1000 > seqtksample2_1000.fq
#用wc查看，可對比前后文件，判斷是否抽取成功
wc -l output_forward_paired.fq
wc -l seqtksample1_1000.fq
spades.py --careful --pe1-1 seqtksample1_1000.fq --pe1-2 seqtksample2_1000.fq -o ./SPAdesout_SRR8491236_1000
這里需要下載seqtk。關(guān)于下載seqtk，可以參考我的簡書文章：
http://www.itdecent.cn/p/316b706a8bc9?v=1667780899292
命令解釋：
gunzip -c #-c或--stdout或--to-stdout 把解壓后的文件輸出到標準輸出設(shè)備。
seqtk sample -s60 從pair end的原始fastq文件中抽取1000條reads。其中-s是seed，控制隨機抽取，但是要注意在抽R1和R2的時候，一定要用相同的seed，這樣才能保證抽出來的R1和R2仍然是配對的，否則有可能會錯位。后面1000表示抽取的reads數(shù)目。
wc -l 統(tǒng)計文件內(nèi)容的行數(shù)
參考文章：https://cloud.tencent.com/developer/article/1674827
我還是出錯了。
看到這個回答，我將SPAdes卸載了，安裝了最新版。我原來安裝的版本是3.12.0，現(xiàn)在安裝的最新版本是3.15.4。
然后我重試，又出現(xiàn)了問題：
根據(jù)這個回答，我輸入如下命令：
spades.py --sc --careful --pe1-1 seqtksample1_1000.fq --pe1-2 seqtksample2_1000.fq -o ./SPAdesout_SRR8491236_1000
結(jié)果如下：
感覺應(yīng)該是成功了吧。
七、Quast評價組裝的基因組效果
Quast安裝可以參考我的簡書：
http://www.itdecent.cn/p/b43d8d4ef14c?v=1667866727396
輸入命令如下：
quast.py ~/SPAdesout_SRR8491236_1000/contigs.fasta -o ~/SPAdesout_SRR8491236_1000/quast_out
結(jié)果如下：
將report.html下載到本地。
report.html的解讀：
（1）contig基本信息統(tǒng)計表
1.contigs 是組裝的contigs的總數(shù)
2.Largest contig 是組裝的contigs中長度最長的contig
3.Total length 是組裝的堿基的總數(shù)
4.Reference length 是參考基因組堿基的總數(shù)
5.GC (%) 是組裝基因組中G和C核苷酸的總數(shù)除以組裝基因組中總的長度
6.Reference GC (%) 是G和C核苷酸在參考基因組中的比例
7.N50 對于一個組裝出來的序列，不論是contig還是scaffold, 首先將各個序列根據(jù)長度從大到小排序，然后從第一個序列開始，將長度進行累加，直到累加的長度超過了總長度的50%，此時，最后一個累加的contig的長度就是N50的長度。N90同理。
8.NG50 對于參考基因組的序列，不論是contig還是scaffold, 首先將各個序列根據(jù)長度從大到小排序，然后從第一個序列開始，將長度進行累加，直到累加的長度超過了總長度的50%，此時，最后一個累加的contig的長度就是N50的長度。
只有在提供參考基因組的情況下，才計算該指標。
9.L50 (Lx, LG50, LGx)是等于或大于N50 (Nx, NG50, NGx)的contigs數(shù)。換句話說，例如，L50是覆蓋程序集一半的最小contigs數(shù)。
L50側(cè)重條數(shù)統(tǒng)計，N50偏重長度統(tǒng)計。
關(guān)于Contig、Scaffold和N50、L50詳細解讀：
Contig是reads拼成的連續(xù)的DNA片段，連續(xù)表達一個gene。通過雙端測序的contig可確定contig之間的關(guān)系得到scaffold，Scaffold是reads拼成的有g(shù)ap的DNA片段。理想情況下，一條染色體用同一個scaffold的表達。整個genome存在很多零碎片段，可舍棄。因為duplication產(chǎn)生很多overlap。
N50，L50和NG50是評價genome assembly的quality的標準，評價長度時使用N50，N50是一個contig的長度。不選用genome size的50%是因為1.這是估計的size值不一定準；2.sequence 僅覆蓋80%。評價數(shù)量使用L50，L50數(shù)量越小越好。
圖片是自己手繪，只能勉強看看。
詳細解讀可以看文章：
https://www.cnblogs.com/yuanjingnan/p/11725496.html
10.N's per 100 kbp即No. of mismatches per 100 kb: 每 100,000 個對齊堿基的平均錯配數(shù)。該指標不區(qū)分單核苷酸多態(tài)性（組裝基因組與參考基因組的真正差異）和單核苷酸錯誤（由于讀取錯誤或組裝算法錯誤）。
11.auN, auNG, auNA, auNGA 分別是Nx, NGx, NAx, NGAx曲線下的面積。
如果您想用一個數(shù)字來總結(jié)基因組組裝的連續(xù)性，auN (auNG等)是比N50 (NG50等)更好的選擇。它更穩(wěn)定，受contig長度大幅變化的影響較小，并考慮整個Nx (NGx等)曲線。
12.Nx (where 0<=x<=100): 最長長度的 contigs 占組裝堿基的百分比。
（2）contig長度累計曲線
橫坐標為contig個數(shù)，縱坐標為累加的長度，示意圖如下
對于所考慮的所有類型的累積圖，重疊群按堿基數(shù)從最大到最小排列。累積長度圖顯示前 x 個 contigs 中的堿基數(shù)，因為 x 從零變化到 contigs 的數(shù)量。類似地計算完整基因的累積數(shù)量和完整操縱子的累積數(shù)量。
（3）Nx 長度分布曲線
顯示隨著 x 變化的 Nx、NGx、NAx、NGAx 指標的趨勢，這比僅使用 N50 提供更多信息。
（4）GC含量分布圖
顯示重疊群中 GC 含量的分布。x 值顯示 Gc 的百分比。y 值顯示 GC 內(nèi)容為 x 的非重疊 100bp 窗口的數(shù)量。這種分布通常是 高斯分布；但是如果存在具有不同 GC 含量的污染物，通常會出現(xiàn)多個 高斯疊加。
這兩張圖作圖依據(jù)是Contigs被分解成不重疊的100 bp窗口。圖中顯示了每個GC百分比的窗口數(shù)。
關(guān)于結(jié)果解答有官方說明：
https://quast.sourceforge.net/docs/manual.html#sec4
也可以參考文章：
https://zerobio.github.io/archives/306037846.html
https://blog.51cto.com/u_10721944/5398083
本文總體參考的文章：
http://www.itdecent.cn/p/1e3fd96aac68
